Flow measurement method based on dynamic optimization of three pressure sensors

[0001] The present application is a continuation-in-part application of U.S. patent application Ser. No. 17/718,383 filed on Apr. 12, 2022, which claims the benefit of Chinese Patent Application No. 2021107141890 filed on Jun. 25, 2021, the disclosure of which is hereby incorporated herein by reference.

[0002] The present disclosure relates to the technical field of measuring flow rate, in particular to a device and method for measuring flow rate based on dynamic optimization of three pressure sensors.

[0003] The existing methods for measuring flow rate, such as differential pressure method, turbine method, vortex method, ultrasonic method and laser method, etc., have some defects. The methods for measuring flow rate, such as differential pressure method, turbine method and vortex method, etc., need an orifice, turbine or vortex generator set in the pipe. However, the orifice, turbine or vortex generator will bring throttling loss to the fluid in the pipe and influence the pattern of the fluid. At the same time, the accuracy of the method for measuring flow rate by turbine or vortex generator is limited by the accuracy of the number of turbine teeth, so it cannot measure the dynamic flow rate of the fluid accurately. Ultrasonic method and laser method are affected by the vibration in the environment, and there is a large error in the vibration environment.

[0004] In order to solve at least one technical problem in the above prior art, the embodiments of the present disclosure provides a device and method for measuring flow rate based on dynamic optimization of three pressure sensors. In the device and method, it can be optimized dynamically based on pressure sensor, the device has a simple and compact structure, and the device and method have high dynamic effect and strong anti-interference ability.

[0005] One aspect, the embodiment of the present disclosure provides a device for measuring flow rate based on dynamic optimization of three pressure sensors. The device comprises a test straight pipe, a flow rate calculation estimator and three pressure sensors;

• • [0006] the test straight pipe has an equal cross-section passage and is configured to provide a passage through which a measured fluid flows; a first pressure sensor and a second pressure sensor are respectively set on a boundary of both ends of the test straight pipe, and are configured to measure boundary pressures of two ends of the fluid in the test straight pipe; and a third pressure sensor is set on the test straight pipe between the first pressure sensor and the second pressure sensor to measure monitoring pressure of the fluid in the middle of the test straight pipe; and
  • [0007] the flow rate calculation estimator is respectively connected to the three pressure sensors, and is configured to acquire a flow rate of the fluid in the test straight pipe according to the boundary pressures and the monitoring pressure.

[0008] Another aspect, the embodiment of the present disclosure provides a method for measuring flow rate based on dynamic optimization of three pressure sensors, and the method is applied to the above device for measuring flow rate based on dynamic optimization of three pressure sensors. The method comprises:

• • [0009] acquiring boundary pressures of a fluid in both ends of a test straight pipe and a monitoring pressure of the fluid in the middle of the test straight pipe;
  • [0010] acquiring a predicted pressure of the fluid in the middle of the test straight pipe according to the boundary pressures;
  • [0011] optimizing a pipe hydrodynamics module according to the predicted pressure and the monitoring pressure to enable the predicted pressure constantly close to the monitoring pressure; and
  • [0012] acquiring a flow rate of the fluid in the test straight pipe according to the boundary pressures and the optimized pipe hydrodynamics module.

[0013] The device for measuring flow rate based on dynamic optimization of three pressure sensors in the embodiment of the present disclosure comprises a test straight pipe, a flow rate calculation estimator and three pressure sensors, has a simple and compact structure and is convenient to assemble. In the method for measuring flow rate based on dynamic optimization of three pressure sensors, pressure sensors are set on the fix positions of both ends of the test straight pipe and on the middle position of the test straight pipe, and are configured to measure the pressure of the fluid at the positions where the pressure sensors are set in real time, and the real-time dynamic flow rate of the fluid in the pipe can be estimated through a flow rate calculation estimator. This method has fast response speed and high measurement accuracy, can effectively eliminate the noise interference of vibration in the environment, and has strong anti-interference ability.

[0022] In the figures: 1—test straight pipe, 21—the first pressure sensor, 22—the second pressure sensor, 23—the third pressure sensor, 3—flow rate calculation estimator.

[0023] The present disclosure will be further described in detail below in conjunction with the drawings and embodiments. It can be understood that the specific embodiments described here are only used to explain related content, but not to limit the present disclosure. In addition, it should be noted that, for ease of description, only the parts related to the present disclosure are shown in the drawings.

[0024] It should be noted that the embodiments in the present disclosure and the features in the embodiments can be combined with each other if there is no conflict. Hereinafter, the present disclosure will be described in detail with reference to the drawings and in conjunction with the embodiments.

[0025] It should be noted that the step number in the present disclosure is only for the convenience of the explanation of the specific embodiment, and is not used to limit the sequence of step execution.

[0026] The method provided by the embodiment of the present disclosure can be executed by the relevant processor, and the following will take the processor as the execution subject as an example. It can be understood that the execution subject can be adjusted according to the actual demand. For example, the execution subject can also be servers, electronic devices, computers, etc.

[0027] One aspect of the embodiment of the present disclosure provides a device for measuring flow rate based on dynamic optimization of three pressure sensors. As shown in FIG. 1, a device for measuring flow rate based on dynamic optimization of three pressure sensors comprises a test straight pipe 1, pressure sensors and a flow rate calculation estimator 3.

[0028] The test straight pipe 1 has an equal cross-section passage and is configured to provide a passage through which the measured fluid flows. That is, the test straight pipe 1 is a straight pipe with fixed length and equal cross-section, which provides a passage through which the measured fluid flows and fixed points for setting the pressure sensors. In the test straight pipe, the fluid flows in from one end and flows out from the other end, and the cross-sectional area through which the fluid flows in the pipe does not change. The flow rate of fluid in the pipe is the object to be measured.

[0029] There are three pressure sensors. Wherein, two pressure sensors are respectively set on the boundary or edge of both ends of the test straight pipe 1, which can be called boundary pressure sensor. The other pressure sensor is set on a certain position on the test straight pipe between the two boundary pressure sensors, which can be called monitoring pressure sensor. The boundary pressure sensors are set on the boundary of the test straight pipe with a fixed distance and the monitoring pressure sensor is located between the boundary pressure sensors. The pressure sensors are configured to measure the fluid pressure at the position where the pressure sensor is set in real time. The fixed distance of the adjacent pressure sensors on the test straight pipe 1, that is, the relative distance between the two adjacent sensors, can be flexibly adjusted by the operator according to the actual measurement needs. Specifically, the first pressure sensor 21 and the second pressure sensor 22 are respectively set on the boundary of the two ends of the test straight pipe 1 to measure the boundary pressures of the two ends of the fluid in the test straight pipe 1, and the third pressure sensor 23 is set on the test straight pipe 1 between the first pressure sensor 21 and the second pressure sensor 22 to measure the monitoring pressure of the fluid in the middle of the test straight pipe 1.

[0030] The flow rate calculation estimator 3 is respectively connected to the three pressure sensors to receive the measured values of the three pressure sensors, and is configured to calculate or acquire the flow rate of the fluid in the test straight pipe 1 according to the boundary pressures and the monitoring pressure. The flow rate solution estimator 3 collects or acquires the fluid pressure values measured by the boundary pressure sensors as the boundary input, and collects or acquire the fluid pressure value measured by the monitoring pressure sensor as the correction basis. The flow rate solution estimator 3 is used for real-time dynamic optimization to estimate the accurate flow rate of the fluid in the pipe. The flow rate calculation estimator 3 can be formed by integrated circuit modules, which are integrated with modules to realize different functions.

[0031] Optionally, as shown in FIG. 2, the flow rate calculation estimator 3 comprises a pipe hydrodynamics module, a dynamic optimization module and a flow rate calculation module, and a pipe hydrodynamics model is integrated in the pipe hydrodynamics module.

[0032] The pipe hydrodynamics module calculates or acquires the predicted pressure of the fluid in the middle of the test straight pipe 1 based on the boundary pressures.

[0033] The dynamic optimization module optimizes the pipe hydrodynamics module based on the predicted pressure and monitoring pressure. That is, the pipe hydrodynamic model integrated in the pipe hydrodynamics module is optimized, so that the predicted pressure is constantly close to the monitoring pressure, that is, the difference between the predicted pressure and monitoring pressure is reduced and finally falls into a reference range. Wherein, the reference range can be adjusted or set according to the actual accuracy requirements. That is, the measured monitoring pressure information is used as the correction basis, and the pipe hydrodynamics module can be optimized by the dynamic optimization algorithm according to the correction basis, so that the estimated value of the monitoring pressure approximates the real value.

[0034] The flow rate calculation module calculate or acquire the flow rate of the fluid in the test straight pipe 1 based on the boundary pressures and the optimized pipe hydrodynamics module.

[0035] As shown in FIG. 3, another aspect of the embodiment of the present disclosure provides a method for measuring flow rate based on dynamic optimization of three pressure sensors. The method is applied to the device for measuring flow rate based on dynamic optimization of three pressure sensors disclosed above, and the method comprises the following steps.

[0036] The boundary pressures of the fluid in the both ends of the test straight pipe 1 and the monitoring pressure of the fluid in the middle of the test straight pipe 1 are acquired. Wherein, the boundary pressures of the fluid in the both ends of the test straight pipe 1 are measured or acquired by the first pressure sensor 21 and the second pressure sensor 22. The monitoring pressure of the fluid in the middle of the test straight pipe 1 is measured or acquired by the third pressure sensor 23.

[0037] The predicted pressure of the fluid in the middle of the test straight pipe 1 is acquired based on the boundary pressures. This step can be executed by the flow rate calculation estimator 3. In this step, the two boundary pressure values are inputted into the pipe hydrodynamics module to calculate the predicted pressure of the fluid in the middle of the test straight pipe 1. The predicted pressure is the pressure of the fluid at the position where the monitoring pressure sensor locates.

[0038] The pipe hydrodynamics module is optimized based on the predicted pressure and monitoring pressure to enable the predicted pressure constantly close to the monitoring pressure. This step can be executed by the flow rate calculation estimator 3. The pipe hydrodynamics module is optimized by the dynamic optimization module based on the predicted pressure and monitoring pressure to enable the predicted pressure constantly close to the monitoring pressure. That is, the measured monitoring pressure information is used as the correction basis, and the pipe hydrodynamics module is optimized dynamically by the dynamic optimization algorithm based on the correction basis to enable the estimated value of the monitoring pressure is close to the real value.

[0039] Based on the boundary pressures and the optimized pipe hydrodynamics module, the flow rate of the fluid in the test straight pipe 1 can be acquired. Through the optimization of the above step, the pipe hydrodynamics module is more accurate, so that the flow rate of the fluid in the test straight pipe 1 can be acquired accurately according to the real-time boundary pressures. In addition, optionally, the acquired flow rate data can be output to the human-computer interaction device for viewing and analysis.

[0040] It should be noted that, for the pipe hydrodynamics module and the flow rate calculation algorithm in the flow rate calculation estimator 3, it is need to set calculation time interval and spatial grid division.

• • [0041] 1) Determining the set positions of the two boundary pressure sensors and the monitoring pressure sensor, measuring the distances between adjacent pressure sensors, and save the set positions and the distances to the flow rate calculation estimator 3.
  • [0042] 2) Setting the parameters of the test straight pipe 1 and the fluid, dividing the spatial grid and setting the algorithm operation time step in the flow rate calculation estimator 3.
  • [0043] 3) Setting the parameters of the dynamic optimization algorithm.

[0044] At the beginning of measurement, the value of fluid pressure and flow rate need to be initialized. However, the optimal estimation algorithm adopted by the flow rate calculation estimator 3 does not depend on the initial value of the flow rate, and the final solution is not affected by the initial value of the flow rate, so that any initial value of flow rate can be accepted. Therefore, the assistance of additional flowmeter does not required.

[0045] Optionally, that the predicted pressure of the fluid in the middle of the test straight pipe 1 is acquired based on the boundary pressures comprises:

• • [0046] the boundary pressures are input into the pipe hydrodynamics module, and the predicted pressure is acquired by solving the pipe hydrodynamics equation according to the method of characteristics. Specifically, the boundary pressures are input into the pipe hydrodynamics module, which comprises the processing and solving of the hydrodynamics equation. For example, firstly, the equation is processed according to the method of characteristics to enable the equation transform to a solvable ordinary differential equation; secondly, the ordinary differential equation is solved to acquire the relationship equation between the predicted pressure and the boundary pressures; finally, the boundary pressures as the boundary condition are input into the pipe hydrodynamics module to calculate the predicted pressure.

[0047] The dynamical partial differential equation of fluid in a pipe includes fluid momentum equation and fluid continuity equation, that is, equation (1) and equation (2):

[0048] $\begin{array}{cc}\frac{\partial q}{\partial t}+\frac{A}{\rho }\frac{\partial p}{\partial x}+f\left(q\right)=0,& \left(1\right)\end{array}$ $\begin{array}{cc}\frac{\partial p}{\partial t}+\frac{\rho c^2}{A}\frac{\partial q}{\partial x}=0,& \left(2\right)\end{array}$

• • [0049] where p and q are respectively a pressure and a flow rate of the fluid in the pipe, ρ is a density of the fluid, A is a cross-sectional area of a flow in the pipe, c is a wave velocity the fluid, f(q) is a friction term of the fluid, which is related to the flow rate of the fluid; x is a coordinate along an axial direction of the pipe, and t is a time.

[0050] The hydrodynamic equation is a complex hyperbolic partial differential equation, which cannot be directly calculated. The method of characteristics is an effective solution for solving the hyperbolic partial differential equation. The process of solving the hyperbolic partial differential equation by using the method of characteristics is as follows.

[0051] Equation (3) and equation (4) can be obtained based on equation (1) and equation (2):

[0052] $\begin{array}{cc}\frac{\partial q}{\partial t}=-\frac{A}{\rho }\frac{\partial p}{\partial x}-f\left(q\right),& \left(3\right)\end{array}$ $\begin{array}{cc}\frac{\partial q}{\partial x}=-\frac{A}{\rho c^2}\frac{\partial p}{\partial t}.& \left(4\right)\end{array}$

[0053] When there is equation (5):

[0054] $\begin{array}{cc}\frac{dq}{\mathrm{dt}}=\frac{\partial q}{\partial t}+\frac{\partial q}{\partial x}\frac{dx}{\mathrm{dt}},& \left(5\right)\end{array}$

• • [0055] equation (6) can be obtained by substituting equation (3) and equation (4) into equation (5):

[0056] $\begin{array}{cc}\frac{dq}{\mathrm{dt}}=-\frac{A}{\rho c^2}\frac{dx}{\mathrm{dt}}\left(\frac{\partial p}{\partial t}+c^2\frac{\mathrm{dt}}{dx}\frac{\partial p}{\partial x}\right)-f\left(q\right).& \left(6\right)\end{array}$

[0057] Items in a bracket on the right side of equation (6) are compared with items on the right side of equation (5), if there is

[0058] $\begin{array}{cc}{c}^2\frac{\mathrm{dt}}{dx}=\frac{dx}{\mathrm{dt}},& \left(7\right)\end{array}$

• • [0059] it means that when

[0060] $\begin{array}{cc}\frac{dx}{\mathrm{dt}}=\pm c,& \left(8\right)\end{array}$

• • [0061] the right side of equation (6) can be converted into fully differential form, namely,

[0062] $\begin{array}{cc}\frac{dq}{\mathrm{dt}}=-\frac{A}{\rho c^2}\frac{dx}{\mathrm{dt}}\left(\frac{\partial p}{\partial t}+\frac{\partial p}{\partial x}\frac{dx}{\mathrm{dt}}\right)-f\left(q\right)=-\frac{A}{\rho c^2}\frac{dx}{\mathrm{dt}}\frac{dp}{\mathrm{dt}}-f\left(q\right).& \left(9\right)\end{array}$

[0063] Equation (8) is a condition under which equation (9) is valid, equation (8) is in a form of a straight line, and thus it is called a characteristic line. The above equation (9) is arranged to obtain the following equations, that is, on the characteristic line, there are

[0064] $\begin{array}{cc}\begin{array}{cc}\frac{dq}{\mathrm{dt}}+\frac{A}{\rho c}\frac{dp}{\mathrm{dt}}+f\left(q\right)=0,& \frac{dx}{\mathrm{dt}}=c,\end{array}& \left(10\right)\end{array}$ $\begin{array}{cc}\begin{array}{cc}\frac{dq}{\mathrm{dt}}-\frac{A}{\rho c}\frac{dp}{\mathrm{dt}}+f\left(q\right)=0,& \frac{dx}{\mathrm{dt}}=-c\end{array}.& \left(11\right)\end{array}$

[0065] The significance of the solution is that on the characteristic line (namely,

[0066] $\frac{dx}{\mathrm{dt}}=\pm c),$
the partial differential equation can be transformed into the ordinary differential equation, so that the hydrodynamic equation can be conveniently solved on the characteristic line. And its actual physical meaning can be expressed as that a pressure wave in the fluid propagates a distance of cΔt to left and right sides within a time interval Δt. The pipe is divided into spatial grids. As shown in FIG. 4, the state of a space point P at a time (t_k+Δt) is calculated from the state of space points A and B at a time t_k through two intersecting characteristic lines c⁺ and c⁻. That is, the state variables at all spatial points in the pipe can be obtained by the method of characteristics.

[0067] After the derivation of equation (3) to equation (11), the dynamic partial differential equation of the fluid in the pipe is solved by the method of characteristics and the computable ordinary differential equation is obtained. The physical meaning of the solution is shown in FIG. 4, that is, a flow rate q_P and a pressure p_P at a certain point P on the measured pipe at the current moment can be calculated based on a flow rate q_A and a pressure p_A at a point A and a flow rate q_B and a pressure p_B at a point B of the measured pipe at the previous moment through the characteristic lines. In order to obtain the calculation equation of the flow rate of the pipe, equation (10) and equation (11) can be written in the form of difference:

[0068] $\begin{array}{cc}{c}^{+}:q_P-q_A+\frac{A}{\rho c}\left(p_P-p_A\right)+\Delta tf\left(q_A\right)=0,& \left(12\right)\end{array}$ $\begin{array}{cc}{c}^{-}:q_P-q_B-\frac{A}{\rho c}\left(p_P-p_B\right)+\Delta tf\left(q_B\right)=0.& \left(13\right)\end{array}$

[0069] By arranging the above equation (12) and equation (13), equations for calculating the flow rate and the pressure at the point P at any position in the measured pipe (non-boundary) can be obtained as follows:

[0070] $\begin{array}{cc}{q}_P=\frac{1}{2}\left[q_A+q_B+\frac{A}{\rho c}\left(p_A-p_B\right)-\Delta t\left(f\left(q_A\right)+f\left(q_B\right)\right)\right],& \left(14\right)\end{array}$ $\begin{array}{cc}{p}_P=\frac{1}{2}\left[p_A+p_B+\frac{\rho c}{A}\left(q_A-q_B\right)-\Delta t\frac{\rho c}{A}\left(f\left(q_B\right)-f\left(q_A\right)\right)\right].& \left(15\right)\end{array}$

[0071] When the point P is located at the boundary of the measured pipe, the situation is respectively special. As shown in FIG. 5, there is only one characteristic line at the left and right boundaries. Therefore, the state of the point P on the boundary needs to be calculated by a characteristic line and external input boundary conditions, and the boundary conditions are pressure values of the left and right boundaries of the measured pipe, which are calculated by the pressure sensor.

[0072] When the point P is located at the left boundary, the value p_P collected by the pressure sensor is taken as the left boundary condition, and by combining with the c⁻ characteristic line, the flow rate at the point P can be obtained as follows:

[0073] $\begin{array}{cc}{q}_P=q_B+\frac{A}{\rho c}\left(p_p-p_B\right)-\Delta tf\left(q_B\right).& \left(16\right)\end{array}$

[0074] Similarly, when the point P is located at the right boundary, the pressure p_P obtained by the pressure sensor is taken as the right boundary condition, at this time, the flow rate at the point P can be obtained as follows:

[0075] $\begin{array}{cc}{q}_P=q_A-\frac{A}{\rho c}\left(p_p-p_A\right)-\Delta tf\left(q_A\right).& \left(17\right)\end{array}$

[0076] Considering the fluid friction term f(q) in equation (1), equations of the approximate calculation model of fluid friction in a laminar flow state are the following equations (18) and (19). The first term on the left side of equations (18) and (19) is fluid steady friction, and the second term is fluid dynamic friction.

[0077] $\begin{array}{cc}f\left(t+\Delta t\right)=\frac{8v}{r^2}q\left(t+\Delta t\right)+\frac{4v}{r^2}{\sum }_{j=1}^4{y}_j\left(t+\Delta t\right),& \left(18\right)\end{array}$ $\begin{array}{cc}{y}_j\left(t+\Delta t\right)=y_j\left(t\right)e^{-\frac{n_{j}v\Delta t}{r^2}}+m_j\left[q\left(t+\Delta t\right)-q\left(t\right)\right]e^{-\frac{n_{j}v\Delta t}{2r^2}},& \left(19\right)\end{array}$

• • [0078] where v is a kinematic viscosity of the fluid, r is a flow radius of the straight pipe, m_j and n_j are weight factors of the j^th term, and there are following equations:

[0079] $\begin{array}{cc}{m}_j=\mu m_{j-1},j\ge 3,& \left(20\right)\end{array}$ $\begin{array}{cc}{n}_j={\mu }^2{n}_{j-1},j\ge 2,& \left(21\right)\end{array}$

• • [0080] where the values of μ, m_j and n_j are shown in Table 1.

[0081]

TABLE 1
Weight factors of friction terms
	μ	m1	m2	n1
	3	2.2457	6.8400	42.849

[0082] Based on the above contents, the calculation of the state of the point P at any position on the pipe is related to a time-space grid division. As shown in FIG. 6, a time-space grid model is established for the pipe, where k represents a time of the calculation, k≥0; and i represents a spatial point of pipe division, 1≤i≤N. Δt represents a time interval, and Δx is a pipe grid space interval.

[0083] Equations (14) and (15) for calculating the internal points of the pipe can be written as expressions related to time k-space i, that is, the state located at the pipe i at the time k can be expressed by the state of adjacent grid points ((i−1) and (i+1)) at the time (k−1) as follows:

[0084] $\begin{array}{cc}{q}_i^k=\frac{1}{2}\left[q_{i-1}^{k-1}+q_{i+1}^{k-1}+\frac{A}{\rho c}\left(p_{i-1}^{k-1}-p_{i+1}^{k-1}\right)-\Delta t\left(f\left(q_{i-1}^{k-1}\right)+f\left(q_{i+1}^{k-1}\right)\right)\right],& \left(22\right)\end{array}$ $\begin{array}{cc}{p}_i^k=\frac{1}{2}\left[p_{i-1}^{k-1}+p_{i+1}^{k-1}+\frac{\rho c}{A}\left(q_{i-1}^{k-1}-q_{i+1}^{k-1}\right)+\Delta t\frac{\rho c}{A}\left(f\left(q_{i+1}^{k-1}\right)-f\left(q_{i-1}^{k-1}\right)\right)\right].& \left(23\right)\end{array}$

[0085] Similarly, when the point P is located at the left and right boundary points, that is, when i=1 and N, equations (16) and (17) can be written as:

[0086] $\begin{array}{cc}{q}_1^k=q_2^{k-1}+\frac{A}{\rho c}\left(p_1^k-p_2^{k-1}\right)-\Delta \mathrm{tf}\left(q_2^{k-1}\right),& \left(24\right)\end{array}$ $\begin{array}{cc}{q}_N^k=q_{N-1}^k-\frac{A}{\rho c}\left(p_N^k-p_{N-1}^{k-1}\right)-\Delta \mathrm{tf}\left(q_{N-1}^{k-1}\right),& \left(25\right)\end{array}$

• • [0087] where

[0088] $p_1^k$
and

[0089] $p_N^k$
are respectively the left boundary condition and the right boundary condition as follows:

[0090] $\begin{array}{cc}{p}_1^k=p_{\mathrm{left}}^k,& \left(26\right)\end{array}$ $\begin{array}{cc}{p}_N^k=p_{r\mathrm{ight}}^k.& \left(27\right)\end{array}$

[0091] Considering the friction term f(q) in the above equation, it can be seen from equations (18) and (19) that

[0092] $\sum {j=1}_4^{}{y}_j\left(t+\Delta t\right)$
(t+Δt) in the friction term also changes with time k-space i, that is, the friction intermediate term y_j(t+Δt) is also taken as

[0093] $y_{j_i}^k.$
Therefore, at the time k, the states at the grid point i of the pipe include

[0094] $q_i^k,p_i^k,y_{1_i}^k,y_{2_i}^k,y_{3_i}^k\mathrm{and}{y}_{4_i}^k.$
Equations (22) and (23) can be written as:

[0095] $\begin{array}{cc}{q}_i^k=\frac{1}{2}\left[q_{i-1}^{k-1}+q_{i+1}^{k-1}+\frac{A}{\rho c}\left(p_{i-1}^{k-1}-p_{i+1}^{k-1}\right)-\Delta t\frac{8\nu }{r^2}\left(q_{i-1}^{k+1}+q_{i+1}^{k-1}\right)-\Delta t\frac{4\nu }{r^2}\left({\sum }_{j=1}^4{y}_{j_{i-1}^{k-1}}+{\sum }_{j=1}^4{y}_{j_{i+1}^{k-1}}\right)\right],& \left(28\right)\end{array}$ $\begin{array}{cc}{p}_i^k=\frac{1}{2}\left[p_{i-1}^{k-1}+p_{i+1}^{k-1}+\frac{\rho c}{A}\left(q_{i-1}^{k-1}-q_{i+1}^{k-1}\right)+\Delta t\frac{\rho c}{A}\frac{8\nu }{r^2}\left(q_{i+1}^{k-1}-q_{i-1}^{k-1}\right)+\Delta t\frac{\rho c}{A}\frac{4\nu }{r^2}\left({\sum }_{j=1}^4{y}_{j_{i+1}^{k-1}}-{\sum }_{j=1}^4{y}_{j_{i-1}^{k-1}}\right)\right],& \left(29\right)\end{array}$

• • [0096] wherein the friction intermediate term

[0097] $y_{j_i}^k,j=1,2,3,4,$
is as follows:

[0098] $\begin{array}{cc}{y}_{j_i}^k=y_{j_i}^{k-1}{e}^{-\frac{n_j\nu \Delta t}{r^2}}+m_j\left(q_i^k-q_i^{k-1}\right)e^{-\frac{n_j\mathrm{\nu \Delta t}}{2r^2}}.& \left(30\right)\end{array}$

[0099] When i=1 and N, equations (32) to (35) of the boundary can be written as:

[0100] $\begin{array}{cc}{q}_1^k=q_2^{k-1}+\frac{A}{\rho c}\left(p_1^k-p_2^{k-1}\right)-\Delta t\frac{8\nu }{r^2}{q}_2^{k-1}-\Delta t\frac{4\nu }{r^2}\left({\sum }_{j=1}^4{y}_{j_2^{k-1}}\right),& \left(31\right)\end{array}$ $\begin{array}{cc}{q}_N^k=q_{N-1}^k-\frac{A}{\rho c}\left(p_N^k-p_{N-1}^{k-1}\right)-\Delta t\frac{8\nu }{r^2}{q}_{N-1}^{k-1}=-\Delta t\frac{4\nu }{r^2}\left({\sum }_{j=1}^4{y}_{j_{N-1}^{k-1}}\right),& \left(32\right)\end{array}$ $\begin{array}{cc}{p}_1^k=p_{\mathrm{left}}^k,& \left(33\right)\end{array}$ $\begin{array}{cc}{p}_N^k=p_{r\mathrm{ight}}^k.& \left(34\right)\end{array}$
Let

[0101] $\alpha =\frac{A}{\rho c},\beta =\Delta t\frac{4\nu }{r^2}\mathrm{and}\gamma =-\frac{\mathrm{\nu \Delta }t}{2r^2},$
arrange equations (36) to (42) to obtain the following equations for calculating the grid state of the pipe.

[0102] Equations for calculating the internal points of the measured pipe are as follows:

[0103] $\begin{array}{cc}{q}_i^k=\frac{1-2\beta }{2}\left(q_{i-1}^{k-1}+q_{i+1}^{k-1}\right)+\frac{\alpha }{2}\left(p_{i-1}^{k-1}-p_{i+1}^{k-1}\right)-\frac{\beta }{2}\left({\sum }_{j=1}^4{y}_{j_{i-1}^{k-1}}+{\sum }_{j=1}^4{y}_{j_{i+1}^{k-1}}\right),& \left(35\right)\end{array}$ $\begin{array}{cc}{p}_i^k=\frac{1-2\beta }{2\alpha }\left(q_{i-1}^{k-1}-q_{i+1}^{k-1}\right)+\frac{1}{2}\left(p_{i-1}^{k-1}+p_{i+1}^{k-1}\right)+\frac{\beta }{2\alpha }\left({\sum }_{j=1}^4{y}_{j_{i+1}^{k-1}}-{\sum }_{j=1}^4{y}_{j_{i-1}^{k-1}}\right).& \left(36\right)\end{array}$

[0104] Equation for calculating the friction intermediate term is as follows:

[0105] $\begin{array}{cc}{y}_{j_i}^k=m_j{e}^{n_j\gamma }{q}_i^k-m_j{e}^{n_j\gamma }{q}_i^{k-1}+e^{2n_j\gamma }{y}_{j_i}^{k-1},j=1,2,3,4,& \left(37\right)\end{array}$

[0106] Equations for calculating the boundary points of the measured pipe are as follows:

[0107] $\begin{array}{cc}{q}_1^k=\left(1-2\beta \right)q_2^{k-1}-\alpha p_2^{k-1}-\beta {\sum }_{j=1}^4{y}_{j_2}^{k-1}+\alpha p_1^k,& \left(38\right)\end{array}$ $\begin{array}{cc}{q}_N^k=\left(1-2\beta \right)q_{N-1}^k+\alpha p_{N-1}^{k-1}-\beta {\sum }_{j=1}^4{y}_{j_{N-1}}^{k-1}+\alpha p_N^k,& \left(39\right)\end{array}$ $\begin{array}{cc}{p}_1^k=p_{\mathrm{left}}^k,& \left(40\right)\end{array}$ $\begin{array}{cc}{p}_N^k=p_{\mathrm{right}}^k.& \left(41\right)\end{array}$

[0108] In practical engineering applications, the time interval Δt of the grid division is generally determined by a hardware system, the space interval Δx is generally set artificially by an algorithm, and the wave velocity c of the fluid is determined by the fluid itself. There are various error sources in the calculation process, it is almost impossible for the characteristic line to pass through the grid points exactly, and thus grid interpolation calculation needs to be performed, as shown in FIG. 7. In general, the space interval Δx of the grid is set in the algorithm, so that the characteristic line falls between two grid points. At this time, the point A at the time (k−1) can be calculated by interpolation based on the adjacent grid points ((i−1) and i), and the point B at the time (k−1) can be calculated by interpolation based on the adjacent grid points (i and (i+1)), thereby compensating the arithmetic error caused by the characteristic line not passing through the grid, and the boundary points are the same.

[0109] Due to the relatively small grid division and considering the simplicity of research, linear interpolation method is adopted. At the time k, the states at the grid point i of the pipe include the flow rate

[0110] $q_i^k,$
the pressure

[0111] $p_i^k,$
the friction intermediate terms

[0112] $y_{1_i}^k,y_{2_i}^k,y_{3_i}^k\mathrm{and}{y}_{4_i}^k.$
Here, for the convenience of description, it is uniformly represented by

[0113] $x_i^k.$
Therefore, the linear interpolation of the points A and B on the characteristic lines that fall between the grid points is calculated as follows:

[0114] $\begin{array}{cc}{x}_{A_i}^{k-1}=\eta x_{i-1}^{k-1}+\left(1-\eta \right)x_i^{k-1},& \left(42\right)\end{array}$ $\begin{array}{cc}{x}_{B_i}^{k-1}=\left(1-\eta \right)x_i^{k-1}+\eta x_{i+1}^{k-1},& \left(43\right)\end{array}$

• • [0115] where

[0116] $\eta =\frac{c\Delta t}{\Delta x}$
represents an interpolation coefficient. When η∈0, 1, the characteristic line must fall between the adjacent grid points, as shown in FIG. 7.

[0117] $q_{i-1}^{k-1},q_{i+1}^{k-1},p_{i-1}^{k-1},p_{i+1}^{k-1},y_{j_{i-1}}^{k-1},y_{j_{i+1}}^{k-1},q_2^{k-1},p_2^{k-1},y_{j_2}^{k-1},q_{N-1}^k,p_{N-1}^k\mathrm{and}{y}_{j_{N-1}}^{k-1}$
in equations (43), (44), (46) and (47) all need to be calculated by interpolation, and are replaced with

[0118] $x_{A_i}^{k-1}\mathrm{and}{x}_{B_i}^{k-1}$
in equations (42) and (43). Arrange and rewrite equations (43) to (49) to obtain the following equations.

[0119] Equations for calculating the internal points of the measured pipe are as follows:

[0120] $\begin{array}{cc}{q}_i^k=\left[\frac{\eta \left(1-2\beta \right)}{2}{q}_{i-1}^{k-1}+\left(1-2\beta \right)\left(1-\eta \right)q_i^{k-1}+\frac{\eta \left(1-2\beta \right)}{2}{q}_{i+1}^{k-1}\right]& \left(44\right)\end{array}$ $\frac{\mathrm{\alpha \eta }}{2}\left[p_{i-1}^{k-1}-p_{i+1}^{k-1}\right]-\left[{\sum }_{j=1}^4\left(\frac{\beta \eta }{2}{y}_{1_{i-1}^{k-1}}+\beta \left(1-\eta \right)y_{1_i^{k-1}}+\frac{\beta \eta }{2}{y}_{1_{i+1}^{k-1}}\right)\right],$ $\begin{array}{cc}{p}_i^k=\left[\frac{\eta \left(1-2\beta \right)}{2\alpha }\left(q_{i-1}^{k-1}-q_{i+1}^{k-1}\right)\right]+\left[\frac{\eta }{2}{p}_{i-1}^{k-1}+\left(1-\eta \right)p_i^{k-1}+\frac{\eta }{2}{p}_{i+1}^{k-1}\right]-\left[\frac{\beta \eta }{2\alpha }{\sum }_{j=1}^4\left(y_{j_{i-1}^{k-1}}-y_{j_{i-1}^{k-1}}\right)\right].& \left(45\right)\end{array}$

[0121] Equation for calculating the friction intermediate term is as follows:

[0122] $\begin{array}{cc}{y}_{j_i}^k=m_j{e}^{n_{j\gamma }}{q}_i^k-m_j{e}^{n_{j\gamma }}{q}_i^{k-1}+e^{2n_{j\gamma }}{y}_{j_i}^{k-1},j=1,2,3,4.& \left(46\right)\end{array}$

[0123] Equations for calculating the boundary points of the measured pipe are as follows:

[0124] $\begin{array}{cc}{q}_1^k=\left(1-2\beta \right)\left[\left(1-\eta \right)q_1^{k-1}+\eta q_2^{k-1}\right]-\alpha \left[\left(1-\eta \right)p_1^{k-1}+\eta p_2^{k-1}\right]-\beta \left[{\sum }_{j=1}^4\left(\left(1-\eta \right)y_{j_1}^{k-1}+\eta y_{j_2}^{k-1}\right)\right]+\alpha p_1^k,& \left(47\right)\end{array}$ $\begin{array}{cc}{q}_N^k=\left(1-2\beta \right)\left[\eta q_{N-1}^k+\left(1-\eta \right)q_N^k\right]+\alpha \left[\eta p_{N-1}^{k-1}+\left(1-\eta \right)q_N^k\right]-\beta \left[{\sum }_{j=1}^4\left(\eta y_{j_{N-1}}^{k-1}+\left(1-\eta \right)y_{j_N}^{k-1}\right)\right]-\alpha p_N^k,& \left(48\right)\end{array}$ $\begin{array}{cc}{p}_1^k=p_{\mathrm{left}}^k,& \left(49\right)\end{array}$ $\begin{array}{cc}{p}_N^k=p_{\mathrm{right}}^k,& \left(50\right)\end{array}$

[0125] Further, it can be seen from the above contents that by taking the boundary pressure value of i=1 or N as the boundary condition, the states of boundary points and internal points are continuously calculated, where the states include the flow rate

[0126] $q_i^k,$
the pressure

[0127] $p_i^k,$
the friction intermediate terms

[0128] $y_{1_i}^k,y_{2_i}^k,y_{3_i}^k\mathrm{and}{y}_{4_i}^k.$
In order to use the optimal estimation method of Kalman filter, it is necessary to establish a state space model of a stochastic system, including state equation and measurement equation.

[0129] First, a (6N×1)-dimensional state vector X^k is constructed as follows:

[0130] $\begin{array}{cc}{X}^k={\left[q^k,p^k,y_1^k,y_2^k,y_3^k,y_4^k\right]}^T.& \left(51\right)\end{array}$

[0131] Equation (51) includes the following N-dimensional flow rate vector q^k, N-dimensional pressure vector p^k and (4N×1)-dimensional friction intermediate term vector:

[0132] $q^k={\left[q_1^k,q_2^k,\dots q_i^k,\dots q_{N-1}^k,q_N^k\right]}^T,$ $p^k={\left[p_1^k,p_2^k,\dots p_i^k,\dots p_{N-1}^k,p_N^k\right]}^T,$ $y_j^k=\left[y_1^k,y_2^k,y_3^k,y_4^k\right]={\left[y_{j_1}^k,y_{j_2}^k,\dots y_{j_i}^k,\dots ,y_{j_{N-1}}^k,y_{j_N}^k\right]}^T,j=1,2,3,4,$

• • [0133] where k represents the time, k≥0; i represents the number of the grid point, 1≤i≤N, that is, it means that the measured pipe is divided into N grid points; j represents the serial number of the friction intermediate term, j=1, 2, 3, 4.

[0134] As shown in FIG. 1, the value collected by the pressure sensor at the middle position of the measured pipe is taken as the measured value, and a (1×1)-dimensional measurement vector Z^k is constructed as follows:

[0135] $\begin{array}{cc}{Z}^k={\left[p_{mid}\right]}^T.& \left(52\right)\end{array}$

[0136] According to the state vector X^k and the measurement vector Z^k, the state space model of the stochastic system is given as follows:

[0137] $\begin{array}{cc}\{\begin{array}{c}{X}^k={\Phi }^{\frac{k}{k-1}}{X}^{k-1}+B^k{u}^k+{\Gamma }^{\frac{k}{k-1}}{W}^{k-1}\\ Z^k=H^k{X}^k+V^k\end{array},& \left(53\right)\end{array}$

• • [0138] where

[0139] ${\Phi }^{\frac{k}{k-1}}$
is a known structural parameter of the system, and is called (6N×6N)-dimensional state one-step transition matrix; u^k is a known deterministic input sequence, and in this model, u^k is a (2×1)-dimensional boundary pressure value; B^k is an input coefficient matrix;

[0140] ${\Gamma }^{\frac{k}{k-1}}$
is a (6N×2)-dimensional noise distribution matrix of the system; W^k−1 is a (2×1)-dimensional noise vector of the system; H^k is a (1×6N)-dimensional measurement matrix, and V^k is a (1×1)-dimensional measurement noise vector.

[0141] For the convenience of calculation, first consider the first term of the state equation of equation (53), and let

[0142] $q^k={\mathrm{MX}}^{k-1},$ $p^k={\mathrm{NX}}^{k-1},$ $y_j^k={\mathrm{Aq}}^k+{\mathrm{DX}}^{k-1}={\mathrm{AMX}}^{k-1}+{\mathrm{DX}}^{k-1}=\left(\mathrm{AM}+D\right)X^{k-1},$

• • [0143] where M, N, A and D are (N×6N)-dimensional matrix, (N×6N)-dimensional matrix, (4N×N)-dimensional matrix and (4N×6N)-dimensional matrix, respectively. The first term of the state equation of equation (53) can be written as:

[0144] $\begin{array}{cc}{X}^k=\left[\begin{array}{c}{q}^k\\ p^k\\ y_j^k\end{array}\right]={\Phi }^{\frac{k}{k-1}}{X}^{k-1}=\left[\begin{array}{c}M\\ N\\ AM+D\end{array}\right]X^{k-1}.& \left(54\right)\end{array}$

[0145] According to the calculation equations (44) to (50) of the fluid of the pipe, the above matrix can be written as:

[0146] $\begin{array}{cc}M=\left[\begin{array}{cccccc}\left(1-2\beta \right)\left(1-\eta \right)& \eta \left(1-2\beta \right)& 0& \dots & 0& 0\\ \frac{\eta \left(1-2\beta \right)}{2}& \left(1-2\beta \right)\left(1-\eta \right)& \frac{\eta \left(1-2\beta \right)}{2}& \dots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& \dots & \left(1-2\beta \right)\left(1-\eta \right)& \frac{\eta \left(1-2\beta \right)}{2}\\ 0& 0& 0& \dots & \eta \left(1-2\beta \right)& \left(1-2\beta \right)\left(1-\eta \right)\end{array}\begin{array}{cccccc}-\alpha \left(1-\eta \right)& -\mathrm{\alpha \eta }& 0& \dots & 0& 0\\ \frac{\mathrm{\alpha \eta }}{2}& 0& -\frac{\mathrm{\alpha \eta }}{2}& \dots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& \dots & 0& -\frac{\mathrm{\alpha \eta }}{2}\\ 0& 0& 0& \dots & \mathrm{\alpha \eta }& \alpha \left(1-\eta \right)\end{array}\underset{\underset{\times 4}{︸}}{\begin{array}{cccccc}-\beta \left(1-\eta \right)& -\mathrm{\beta \eta }& 0& \dots & 0& 0\\ \frac{-\mathrm{\beta \eta }}{2}& -\beta \left(1-\eta \right)& \frac{-\mathrm{\beta \eta }}{2}& \dots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& \dots & -\beta \left(1-\eta \right)& \frac{-\mathrm{\beta \eta }}{2}\\ 0& 0& 0& \dots & -\mathrm{\beta \eta }& -\beta \left(1-\eta \right)\end{array}}\right],& \left(55\right)\end{array}$ $\begin{array}{cc}M=\left[\begin{array}{cccccc}0& 0& 0& \dots & 0& 0\\ \frac{\eta \left(1-2\beta \right)}{2\alpha }& 0& \frac{\eta \left(2\beta -1\right)}{2\alpha }& \dots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& \dots & 0& \frac{\eta \left(2\beta -1\right)}{2\alpha }\\ 0& 0& 0& \dots & 0& 0\end{array}\begin{array}{cccccc}0& 0& 0& \dots & 0& 0\\ \frac{\eta }{2}& \left(1-\eta \right)& \frac{\eta }{2}& \dots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& \dots & \left(1-\eta \right)& \frac{\eta }{2}\\ 0& 0& 0& \dots & 0& 0\end{array}\underset{\underset{\times 4}{︸}}{\begin{array}{cccccc}0& 0& 0& \dots & 0& 0\\ \frac{-\mathrm{\beta \eta }}{2\alpha }& 0& \frac{\mathrm{\beta \eta }}{2\alpha }& \dots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& \dots & 0& \frac{\mathrm{\beta \eta }}{2\alpha }\\ 0& 0& 0& \dots & 0& 0\end{array}}\right],& \left(56\right)\end{array}$ $\begin{array}{cc}A=\left[\begin{array}{ccc}{m}_1{e}^{n_1\gamma }& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & m_1{e}^{n_1\gamma }\\ m_2{e}^{n_2\gamma }& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & m_2{e}^{n_2\gamma }\\ m_3{e}^{n_3\gamma }& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & m_3{e}^{n_3\gamma }\\ m_4{e}^{n_4\gamma }& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & m_4{e}^{n_4\gamma }\end{array}\right],& \left(57\right)\end{array}$ $\begin{array}{cc}D=\left[\begin{array}{cccccc}\begin{array}{ccc}-m_1{e}^{n_1\gamma }& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & -m_1{e}^{n_1\gamma }\end{array}& 0_{N\times N}& \begin{array}{ccc}{e}^{2n_1\gamma }& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & e^{2n_1\gamma }\end{array}& 0_{N\times N}& 0_{N\times N}& 0_{N\times N}\\ \begin{array}{ccc}-m_2{e}^{n_2\gamma }& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & -m_2{e}^{n_2\gamma }\end{array}& 0_{N\times N}& 0_{N\times N}& \begin{array}{ccc}{e}^{2n_2\gamma }& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & e^{2n_2\gamma }\end{array}& 0_{N\times N}& 0_{N\times N}\\ ⋮& 0_{N\times N}& 0_{N\times N}& 0_{N\times N}& \ddots & 0_{N\times N}\\ \begin{array}{ccc}-m_4{e}^{n_4\gamma }& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & -m_4{e}^{n_4\gamma }\end{array}& 0_{N\times N}& 0_{N\times N}& 0_{N\times N}& 0_{N\times N}& \begin{array}{ccc}{e}^{2n_4\gamma }& \dots & 0\\ ⋮& \ddots & ⋮\\ 0& \dots & e^{2n_4\gamma }\end{array}\end{array}\right],& \left(58\right)\end{array}$

• • [0147] where α, β, γ have been given in the above contents; η is an interpolation coefficient; m_j, n_j, j=1, 2, 3, 4 are weight factors, and the equations are (20) and (21).

[0148] Other matrices of the state space model (53) of the stochastic system are as follows:

[0149] $\begin{array}{cc}{u}^k=\left[\begin{array}{c}{p}_{\mathrm{left}}\\ p_{ri\mathrm{ght}}\end{array}\right],B^k=\left[\begin{array}{cc}\alpha & 0\\ ⋮& ⋮\\ 0& -\alpha \\ 1& 0\\ ⋮& ⋮\\ 0& 1\\ 0_{4N\times 1}& 0_{4N\times 1}\end{array}\right],& \left(59\right)\end{array}$ $\begin{array}{cc}{W}^{k-1}=\left[\begin{array}{c}{w}_{\mathrm{left}}^{k-1}\\ w_{ri\mathrm{ght}}^{k-1}\end{array}\right],{\Gamma }^{\frac{k}{k-1}}=\left[\begin{array}{cc}1& 0\\ ⋮& ⋮\\ 0& 1\\ 1& 0\\ ⋮& ⋮\\ 0& 1\\ 0_{4N\times 1}& 0_{4N\times 1}\end{array}\right],& \left(60\right)\end{array}$ $\begin{array}{cc}{H}^k=\left[\begin{array}{ccc}{0}_{1\times N}& \stackrel{\stackrel{N}{︷}}{\begin{array}{ccc}\underset{\underset{\frac{N}{2}}{︸}}{\begin{array}{ccc}0& \dots & 1\end{array}}& \dots & 0\end{array}}& 0_{4\times 1}\end{array}\right]V^k=\left[w_{mid}^{k-1}\right],& \left(61\right)\end{array}$

• • [0150] where

[0151] $w_{\mathrm{left}}^{k-1}{w}_{\mathrm{right}}^{k-1}{w}_{\mathrm{mid}}^{k-1}$
are measured white noises measured by the pressure sensors disposed at the left boundary, the right boundary and the middle position of the pipe at the time (k−1), that is, satisfy the following equation:

[0152] $\begin{array}{cc}\{\begin{array}{cc}E\left[W_k\right]=0,& E\left[W_k{W}_j^T\right]=Q_k\\ E\left[V_k\right]=0,& E\left[V_k{V}_j^T\right]=R_k\\ E\left[W_k{V}_j^T\right]=0& \end{array}.& \left(62\right)\end{array}$

[0153] It can be seen from the above contents, it can be seen that the state space model of the stochastic system is constructed by equations (51) to (62), mainly the state equation and measurement equation of equation (53). The state equation contains the calculation equation of the grid fluid state of the pipe, and describes that the states (the flow rate, the pressure and the friction intermediate term) at respective positions of the measured pipe are calculated in real time based on the boundary input conditions (the pressures of left and right boundaries). The measurement equation (that is, information of the pressure sensor in the middle position of the measured pipe is obtained as measurement values to construct the equation) continuously optimizes the states in the state equation in real time.

[0154] Optionally, that the pipe hydrodynamics module is optimized based on the predicted pressure and monitoring pressure to enable the predicted pressure constantly close to the monitoring pressure comprises:

• • [0155] based on the predicted pressure and monitoring pressure, by optimizing the algorithm module, calculating the uncertainty of the algorithm model and the observation model under the experimental conditions, acquiring the weight value of the monitoring pressure and observation pressure based on the uncertainty, performing the weighted average between the predicted pressure and the monitoring pressure, the predicted pressure value with better accuracy can be acquired and the pipe hydrodynamics module can be corrected, so that the predicted pressure is constantly close to the accurate pipe pressure which is the monitoring pressure.

[0156] In other words, substituting the boundary pressure values acquired by the first and second pressure sensors into the pipe hydrodynamics model to acquire the predicted pressure at the position of the third pressure sensor on the test straight pipe; combining with the optimization algorithm model, correcting the parameters of the model continuously in the process of solving, so that the predicted pressure is continuously close to the monitoring pressure acquired by the third pressure sensor, and the optimization process will come to an end until the error is within the acceptable range.

[0157] Optionally, that based on the boundary pressures and the pipe hydrodynamics module, the flow rate of the fluid in the test straight pipe is acquired comprises:

• • [0158] based on the boundary pressures and the optimized pipe hydrodynamics module, the predicted pressure and the initial flow rate of the fluid in the straight pipe are calculated firstly; then, while continuously correcting the predicted pressure and optimizing the pipe hydrodynamics module, the flow rate of the fluid in the straight pipe calculated from the pipe hydrodynamics equation is constantly close to the real value.

[0159] It can be seen from the above description that the principle of the device and method for measuring flow rate based on dynamic optimization of three pressure sensors in the embodiments of the present disclosure is to optimally estimate the current flow rate value of the fluid in the test straight pipe 1 according to the pressure values of the boundary pressure sensors and the monitoring pressure sensor at the current time. Specifically, the principle is as follows.

• • [0160] 1) The fluid pressure value measured by the boundary pressure sensors at the current time is acquired or collected;
  • [0161] 2) the boundary pressures acquired or collected by the flow rate calculation estimator 3 are substituted into the pipe hydrodynamics module as an external condition, and the pressure at the position of the monitoring pressure sensor at the current time is calculated according to the pipe hydrodynamics module;
  • [0162] 3) the real fluid pressure measured by the monitoring pressure sensor at the current time is acquired or collected;
  • [0163] 4) in the flow rate calculation estimator 3, the real monitoring pressure is taken as a reference, and the real monitoring pressure is updated using the dynamic optimization algorithm model, until the pressure calculated in step 2) is close to that calculated in step 3); and
  • [0164] 5) based on the boundary pressure measurement values at both ends, the flow rate of the test straight pipe 1 is accurately calculated by using the dynamic optimized pipe hydrodynamics module.

[0165] In the description of this specification, the descriptions of the terms “one embodiment/mode”, “some embodiments/modes”, “examples”, “specific examples”, or “some examples” means that the specific features, structures, materials, or characteristics described in connection with the embodiment/mode or example are included in at least one embodiment/mode or example of the present application. In this specification, the schematic expression of the above terms does not necessarily refer to the same embodiment/mode or example. Moreover, the specific features, structures, materials, or characteristics described can be combined in any suitable manner in any one or more embodiments/modes or examples. In addition, without contradicting each other, those skilled in the art may combine different embodiments/modes or examples and features of the different embodiments/modes or examples described in this specification.

[0166] In addition, the terms “first” and “second” are used for description purposes only, and cannot be understood as indicating or implying relative importance or implicitly indicating the number of technical features indicated. Thus, the features defined as “first” and “second” may include at least one of the features either explicitly or implicitly. In the description of the present application, the meaning of “plurality” is at least two, such as two, three, etc., unless specifically defined otherwise.

[0167] Those skilled in the art should understand that the above-mentioned embodiments are only for clearly illustrating the present disclosure, rather than limiting the scope of the present disclosure. For those skilled in the art, other changes or modifications can be made on the basis of the above disclosure, and these changes or modifications are still within the scope of the present disclosure.