NATURAL GAS CIRCUIT MODELING METHOD FOR OPERATION CONTROL OF INTEGRATED ENERGY SYSTEM

Abstract

A natural gas circuit modeling method for operation control of an integrated energy system, including: establishing, based on the mass conservation and momentum conservation equations in the natural gas pipeline, as well as the state equation and the flow equation of the natural gas, a partial differential equation between the flow rate and the pressure in the natural gas pipeline; mapping the gas circuit to the frequency domain through Fourier transform and obtaining the lumped parameter model through the two-port equivalence; establishing, combined with the natural gas compressor equation, the general branch model of the natural gas circuit; defining the node-branch correlation matrix and the node-outflow branch correlation matrix to establish the topological constraint equation of the natural gas circuit; and establishing, combined with the general branch model of the natural gas circuit and the topological constraint equation of the natural gas circuit, the natural gas circuit equation.

Claims

1. A natural gas circuit modeling method for operation control of an integrated energy system, comprising: step 1 of establishing a pipeline model of a natural gas circuit, comprising: step 1-1 of establishing a mass conservation equation and a momentum conservation equation for a one-dimensional flow process of the natural gas in the pipeline: $\frac{\partial \rho }{\partial t}+\frac{\partial \rho v}{\partial x}=0,\mathrm{and}$ $\frac{\partial \rho v}{\partial t}+\frac{\partial \rho v^2}{\partial x}+\frac{\partial p}{\partial x}+\frac{\lambda \rho v^2}{2D}+\rho g\sin \theta =0,$ where ρ, v and p denote a density, a flow rate and a pressure of the natural gas, respectively, λ, D and θ denote a friction coefficient, an inner diameter and an inclination angle of the pipeline, respectively, provided by a management party of the natural gas circuit, g denotes an acceleration of gravity, and t and x denote time and space, respectively; step 1-2 of introducing two approximations into the momentum conservation equation in the step 1-1: one approximation is to ignore the convection term, that is $\frac{\partial \rho v^2}{\partial x}\approx 0;$ the other approximation is to perform incremental linearization approximation on the square term of the flow rate in the resistance term: that is v.sup.2≈2v.sub.bv−v.sub.b.sup.2, where v.sub.b denotes a base value of the flow rate of the natural gas in the natural gas pipeline, which is a flow rate in a design condition, to obtain the resistance term $\frac{\lambda \rho v^2}{2D}\approx \frac{\lambda \rho \left(2v_{b}v-v_b^2\right)}{2D}$ in the momentum conservation equation in the step 1-1, and simplifying the momentum conservation equation as: $\frac{\partial \rho v}{\partial t}+\frac{\lambda \rho 2v_{b}v-v_b^2}{2D}+\frac{\partial p}{\partial x}+\rho g\sin \theta =0;$ step 1-3 of substituting a state equation of the natural gas p=RTρ and a flow equation of the pipeline G=ρvA into the mass conservation equation and the simplified momentum conservation equation, and obtaining a space-time partial differential equation between a flow and the pressure of the natural gas in the pipeline: $A\frac{\partial p}{\partial t}+RT\frac{\partial G}{\partial x}=0,\mathrm{and}$ $\frac{1}{A}\frac{\partial G}{\partial t}+\frac{\partial p}{\partial x}+\frac{\lambda v_{b}G}{AD}-\frac{\lambda v_b^{2}p}{2RTD}+\frac{pg\sin \theta }{RT}=0,$ where R and T denote a gas constant and temperature of the natural gas, respectively, G denotes a mass flow of the natural gas, and A denotes a cross-sectional area of the natural gas pipeline; step 1-4 of establishing equations of a flow difference and a pressure drop at both ends of a micro-element on the natural gas pipeline: $\mathrm{dG}=-\frac{\mathrm{Adx}}{RT}.Math.\frac{\mathrm{dp}}{dt},\mathrm{and}$ $\mathrm{dp}=-\frac{dx}{A}.Math.\frac{\mathrm{dG}}{\mathrm{dt}}-\frac{\lambda v_b\mathrm{dx}}{AD}.Math.G-\frac{2\mathrm{gD}\sin \theta -\lambda v_b^2}{2RTD}\mathrm{dx}.Math.p,$ step 1-5 of defining, based on the equations of the flow difference and the pressure drop at both ends of the micro-element in the step 1-4, calculation equations of a gas resistance R.sub.g, a gas inductance L.sub.g, a gas capacitance C.sub.g and a controlled gas pressure source k.sub.g, in the natural gas pipeline as follows: $R_g=\lambda v_b/\left(\mathrm{AD}\right),$ $L_g=1/A,$ $C_g=A/\left(\mathrm{RT}\right),\mathrm{and}$ $k_g=\frac{2\mathrm{gD}\sin \theta -\lambda v_b^2}{2RTD},$ wherein a pipeline with a length of dx is represented as a gas circuit comprising four elements, and the entire pipeline is further represented as a distributed parameter gas circuit; step 1-6 of substituting R.sub.g, L.sub.g, C.sub.g and k.sub.g defined in the step 1-5 into the equations of the flow difference and the pressure drop at both ends of the micro-element in the step 1-4, mapping them to a frequency domain through Fourier transform, and obtaining an ordinary differential equation for each frequency component as follows: $\frac{\mathrm{dG}}{\mathrm{dx}}=-jw{C}_g.Math.p,\mathrm{and}$ $\frac{\mathrm{dp}}{\mathrm{dx}}=-\left(R_g+{\mathrm{wL}}_g\right)G-k_{g}p,$ and defining Z.sub.g=R.sub.g+jwL.sub.g, and Y.sub.g=jwC.sub.g; step 1-7 of using the two equations in the step 1-6 to solve a flow and a pressure at a terminal end of the natural gas pipeline as: $G_l=\left[-\frac{2Y_g{p}_0}{\sqrt{k_g^2+4Z_g{Y}_g}}\sinh \left(\frac{\sqrt{k_g^2+4Z_g{Y}_g}}{2}l\right)+G_0\cosh \left(\frac{\sqrt{k_g^2+4Z_g{Y}_g}}{2}l\right)+\frac{k_g{G}_o}{\sqrt{k_g^2+4Z_g{Y}_g}}\sinh \left(\frac{\sqrt{k_g^2+4Z_g{Y}_g}}{2}l\right)\right]e^{-\frac{k_{g}l}{2}},\mathrm{and}$ $p_l=\left[\text{⁠}{p}_0\cosh \left(\frac{\sqrt{k_g^2+4Z_g{Y}_g}}{2}l\right)-\frac{2k_g{G}_0}{\sqrt{k_g^2+4Z_g{Y}_g}}\sinh \left(\frac{\sqrt{k_g^2+4Z_g{Y}_g}}{2}l\right)-\frac{k_g{p}_0}{\sqrt{k_g^2+4Z_g{Y}_g}}\sinh \left(\frac{\sqrt{k_g^2+4Z_g{Y}_g}}{2}l\right)\right]e^{-\frac{k_{g}l}{2}},$ where G.sub.l and p.sub.l denote the flow of the natural gas and the pressure of the natural gas at the terminal end of the natural gas pipeline, respectively, G.sub.o and p.sub.o denote a flow of the natural gas and a pressure of the natural gas at a first end of the natural gas pipeline, respectively, and l denotes a length of the natural gas pipeline; step 1-8 of defining a propagation coefficient of the natural gas pipeline as γ.sub.gc=Z.sub.gY.sub.g, and defining a characteristic impedance of the natural gas pipeline as Z.sub.gc=Z.sub.g/Y.sub.g; step 1-9 of expressing, based on the two equations of the flow and the pressure at the terminal end of the natural gas pipeline in the step 1-7 and the two definitions in the step 1-8, the natural gas pipeline equation in a form of a linear two-port network: $\left[\begin{array}{c}{p}_l\\ G_l\end{array}\right]=\left[\begin{array}{cc}A& B\\ C& D\end{array}\right]\left[\begin{array}{c}{p}_0\\ G_0\end{array}\right],$ where A, B, C and D denote network parameters whose values are: $A=\left[\cosh \left(\frac{\sqrt{k_g^2+4{\gamma }_{\mathrm{gc}}}}{2}l\right)-\frac{k_g}{\sqrt{k_g^2+4{\gamma }_{\mathrm{gc}}}}\sinh (\frac{\sqrt{k_g^2+4{\gamma }_{\mathrm{gc}}}}{2}l\right].Math.e^{-\frac{k_{g}l}{2}},$ $B=-\frac{2}{\sqrt{k_g^2+4/Z_{gc}}}\sinh \left(\frac{\sqrt{k_g^2+4{\gamma }_{gc}}}{2}l\right).Math.e^{-\frac{k_{g}l}{2}},$ $C=-\frac{2}{\sqrt{k_g^2+4Z_{gc}}}\sinh \left(\frac{\sqrt{k_g^2+4{\gamma }_{gc}}}{2}l\right).Math.e^{-\frac{k_{g}l}{2}},\mathrm{and}$ $D=\left[\cosh \left(\frac{\sqrt{k_g^2+4{\gamma }_{gc}}}{2}l\right)+\frac{k_g}{\sqrt{k_g^2+4{\gamma }_{gc}}}\sinh (\frac{\sqrt{k_g^2+4{\gamma }_{gc}}}{2}l\right].Math.e^{-\frac{k_{g}l}{2}},$ step 1-10 of establishing, based on the linear two-port network equation in the step 1-9, a π_type equivalent gas circuit, and equivalent parameters are:
Z=−B,
K=1−AD−BC,
Y.sub.1=(AD−BC−A)/B, and
Y.sub.2=(1−D)/B; step 2 of establishing a general branch model of the natural gas circuit, comprising: step 2-1 of establishing a mathematical model of a natural gas compressor as follows:
p.sub.1=p.sub.2+E.sub.g, where p.sub.1 and p.sub.2 denote pressures on both sides of the natural gas compressor, and E.sub.g denotes a pressure increment provided by the natural gas compressor; step 2-2 of forming, based on the gas resistance, the gas inductance, the gas capacitance, the controlled gas pressure source in the step 1-5 and the natural gas compressor in the step 2-1, a general branch, a general branch equation being as follows:
G.sub.b=y.sub.b(p.sub.b+E.sub.b−k.sub.bp.sub.f), where G.sub.b denotes a flow in the branch, G.sub.b is an unknown quantity, p.sub.b is a gas pressure difference between two ends of the branch, p.sub.f is a pressure at a first end of the branch, when the first end of the branch is a gas source, p.sub.f is a known quantity, when the first end of the branch is not the gas source, p.sub.f is an unknown quantity, p.sub.t is a pressure at a terminal end of the branch, p.sub.t is an unknown quantity, y.sub.b is a branch admittance composed of the gas resistance, the gas inductance and the gas capacitance, k.sub.b and E.sub.b denote a component parameter of a controlled gas pressure source in the branch and a component parameter of the natural gas compressor, provided by the natural gas circuit management party; step 2-3 of writing branch equations of all branches in the natural gas circuit into a matrix form as follows:
G.sub.b=y.sub.b(p.sub.b+E.sub.b−k.sub.bp.sub.f), where G.sub.b denotes a vector formed by the flow of each branch, y.sub.b denotes a diagonal matrix formed by each branch admittance, p.sub.b denotes a vector formed by the pressure difference at both ends of each branch, E.sub.b is a vector formed by the pressure increment of each natural gas compressor, k.sub.b is a vector formed by the parameter of each controlled gas pressure source, and p.sub.f is a vector formed by the gas pressure at the first end of each branch; step 3 of establishing a topological constraint equation of the natural gas circuit, comprising: step 3-1 of defining a node-branch correlation matrix A.sub.g in the natural gas circuit, wherein the matrix is a matrix of n rows and m columns, where n denotes a number of nodes and m denotes a number of branches, and (A.sub.g).sub.i,j denotes an element of the i-th row and the j-th column, (A.sub.g).sub.i,j=0 denotes that a branch j is not connected to a node i, (A.sub.g).sub.i,j=1 denotes the branch j flowing out from the node i, and (A.sub.g).sub.i,j=−1 denotes the branch j flowing into the node i; step 3-2 of defining a node-outflow branch correlation matrix A.sub.g+ in the natural gas circuit, wherein the matrix retains an non-negative element in the matrix A.sub.g, that is, as for (A.sub.g+).sub.i,j, when the branch j flows out from the node i, the element is 1, and when the branch j does not flow out from the node i, the element is 0; step 3-3 of establishing a mass conservation equation of nodes of the natural gas circuit:
A.sub.gG.sub.g=G.sub.n, where G.sub.n denotes a column vector formed by a flow injection at each node, where a flow at a gas load node in the natural gas circuit is a known quantity, a flow at a gas source node is an unknown quantity, and a flow at a node of a non-gas load and a non-gas source is 0; step 3-4 of establishing a gas pressure equation of nodes of the natural gas circuit:
A.sub.g.sup.Tp.sub.n=p.sub.b, and
A.sup.t.sub.g+p.sub.n=p.sub.f, where p.sub.n denotes a column vector formed by pressures at each node, where a pressure at the gas source node in the natural gas circuit is a known quantity, a pressure at the gas load node is an unknown quantity, a pressure at the node of the non-gas load and the non-gas source is an unknown quantity; step 4 of establishing a natural gas circuit equation, comprising: step 4-1 of substituting the equations established in the steps 3-3 and 3-4 into the branch equation established in the step 2-3, and obtaining a natural gas circuit equation in an unreduced form as follows:
G.sub.n=A.sub.gy.sub.b(A.sub.g.sup.Tp.sub.n+E.sub.b−k.sub.bA.sub.g+.sup.Tp.sub.n), step 4-2 of defining a generalized node admittance matrix Y′.sub.g and a generalized node injection vector G′.sub.n as follows:
Y′.sub.g=A.sub.gy.sub.bA.sub.g.sup.T−A.sub.gy.sub.bk.sub.bA.sub.g+.sup.T, and
G′.sub.n=G.sub.n−A.sub.gy.sub.bE.sub.b; step 4-3 of substituting Y′.sub.g and G′.sub.n defined in the step 4-2 into the unreduced natural gas circuit equation in the step 4-1 and obtaining a following natural gas circuit model equation:
Y′.sub.gp.sub.n=G′.sub.n, solving the above natural gas circuit model to obtain an unknown node pressure in the natural gas circuit, and obtaining an unknown branch flow by using the branch equation to implement the operation control of the integrated energy system.

Description

BRIEF DESCRIPTION OF DRAWINGS

[0043] FIG. 1 is a diagram of a distributed parameter gas circuit of a natural gas pipeline, in which FIG. 1a is a schematic diagram of a distributed parameter gas circuit of the entire natural gas pipeline, and FIG. 1B is a schematic diagram of a distributed parameter gas circuit of a micro-element dx of the natural gas pipeline.

[0044] FIG. 2 is a schematic diagram of a lumped parameter equivalent gas circuit of a natural gas pipeline.

[0045] FIG. 3 is a schematic diagram of a general branch in a natural gas circuit.

DESCRIPTION OF EMBODIMENTS

[0046] The natural gas circuit modeling method for operation control of an integrated energy system proposed in the present disclosure includes the following steps.

[0047] Step 1 of establishing a pipeline model of a natural gas circuit, including:

[0048] step 1-1 of establishing a mass conservation equation and a momentum conservation equation for a one-dimensional flow process of the natural gas in the pipeline:

[00012]$\frac{\partial \rho }{\partial t}+\frac{\partial \rho v}{\partial x}=0,\mathrm{and}$$\frac{\partial \rho v}{\partial t}+\frac{\partial \rho v^2}{\partial x}+\frac{\partial p}{\partial x}+\frac{\lambda \rho v^2}{2D}+\rho g\sin \theta =0,$

[0049] where ρ, v and p denote a density, a flow rate and a pressure of the natural gas, respectively, λ, D and θ denote a friction coefficient, an inner diameter and an inclination angle of the pipeline, respectively, which are provided by a management party of the natural gas circuit, g denotes an acceleration of gravity, and t and x denote time and space, respectively;

[0050] step 1-2 of introducing two approximations into the momentum conservation equation in the step 1-1: one approximation is to ignore the convection term, that is

[00013]$\frac{\partial \rho v^2}{\partial x}\approx 0;$

the other approximation is to perform incremental linearization approximation on the square term of the flow rate in the resistance term: that is v.sup.2≈2v.sub.bv−v.sub.b.sup.2, where v.sub.b denotes a base value of the flow rate of the natural gas in the natural gas pipeline, which is a flow rate in a design condition, to obtain the resistance term

[00014]$\frac{\mathrm{\lambda \rho }{v}^2}{2D}\approx \frac{\mathrm{\lambda \rho }\left(2v_{b}v-v_b^2\right)}{2D}$

in the momentum conservation equation in the step 1-1, and simplifying the momentum conservation equation as:

[00015]$\frac{\partial \rho v}{\partial t}+\frac{\mathrm{\lambda \rho }\left(2v_{b}v-v_b^2\right)}{2D}+\frac{\partial p}{\partial x}+\rho g\sin \theta =0;$

[0051] step 1-3 of substituting a state equation of the natural gas p=RTρ and a flow equation of the pipeline G=ρvA into the mass conservation equation and the simplified momentum conservation equation, and obtaining a space-time partial differential equation between a flow and the pressure of the natural gas in the pipeline:

[00016]$A\frac{\partial p}{\partial t}+RT\frac{\partial G}{\partial x}=0,\mathrm{and}$$\frac{1}{A}\frac{\partial G}{\partial t}+\frac{\partial p}{\partial x}+\frac{\lambda v_{b}G}{AD}-\frac{\lambda v_b^{2}p}{2RTD}+\frac{pg\sin \theta }{RT}=0,$

[0052] where R and T denote a gas constant and temperature of the natural gas, respectively, G denotes a mass flow of the natural gas, and A denotes a cross-sectional area of the natural gas pipeline;

[0053] step 1-4 of establishing equations of a flow difference and a pressure drop at both ends of a micro-element on the natural gas pipeline:

[00017]$\mathrm{dG}=-\frac{Adx}{RT}.\frac{\mathrm{dp}}{\mathrm{dt}},\mathrm{and}$$\mathrm{dp}=-\frac{\mathrm{dx}}{A}.Math.\frac{\mathrm{dG}}{\mathrm{dt}}-\frac{\lambda v_b\mathrm{dx}}{AD}.Math.G-\frac{2gD\sin \theta -\lambda v_b^2}{2RTD}\mathrm{dx}.Math.p,$

[0054] step 1-5 of defining, based on the equations of the flow difference and the pressure drop at both ends of the micro-element in the step 1-4, calculation equations of a gas resistance R.sub.g, a gas inductance L.sub.s, a gas capacitance C.sub.g and a controlled gas pressure source k.sub.g, in the natural gas pipeline as follows:

R.sub.g=λv.sub.b/(AD),

L.sub.g=1/A,

C.sub.g=A/(RT), and

[00018]$k_g=\frac{2gD\sin \theta -\lambda v_b^2}{2\mathrm{RTD}},$

[0055] where a pipeline with a length of dx is represented as a gas circuit including four elements, and the entire pipeline is further represented as a distributed parameter gas circuit, and the distributed parameter gas circuit of the entire natural gas pipeline and the micro-element dx of the distributed parameter gas circuit of the natural gas pipeline are shown in FIG. 1;

[0056] step 1-6 of substituting R.sub.g, L.sub.g, C.sub.g and k.sub.g defined in the step 1-5 into the equations of the flow difference and the pressure drop at both ends of the micro-element in the step 1-4, mapping them to a frequency domain through Fourier transform, and obtaining an ordinary differential equation for each frequency component as follows:

[00019]$\frac{\mathrm{dG}}{\mathrm{dx}}=-jw{C}_g.Math.p,\mathrm{and}$$\frac{\mathrm{dp}}{\mathrm{dx}}=-\left(R_g+jw{L}_g\right)G-k_{g}p,$

[0057] and defining Z.sub.g=R.sub.g+jwL.sub.g, and Y.sub.g=jwC.sub.g;

[0058] step 1-7 of using the two equations in the step 1-6 to solve a flow and a pressure at a terminal end of the natural gas pipeline as:

[00020]$G_l=\left[-\frac{2Y_g{p}_0}{\sqrt{k_g^2+4Z_g{Y}_g}}\sinh \left(\frac{\sqrt{k_g^2+4Z_g{Y}_g}}{2}l\right)+G_0\cosh \left(\frac{\sqrt{k_g^2+4Z_g{Y}_g}}{2}l\right)+\frac{k_g{G}_o}{\sqrt{k_g^2+4Z_g{Y}_g}}\sinh \left(\frac{\sqrt{k_g^2+4Z_g{Y}_g}}{2}l\right)\right]e^{-\frac{k_{g}l}{2}},$$\mathrm{and}{p}_l=\left[\text{⁠}{p}_0\cosh \left(\frac{\sqrt{k_g^2+4Z_g{Y}_g}}{2}l\right)-\frac{2k_g{G}_0}{\sqrt{k_g^2+4Z_g{Y}_g}}\sinh \left(\frac{\sqrt{k_g^2+4Z_g{Y}_g}}{2}l\right)-\frac{k_g{p}_0}{\sqrt{k_g^2+4Z_g{Y}_g}}\sinh \left(\frac{\sqrt{k_g^2+4Z_g{Y}_g}}{2}l\right)\right]e^{-\frac{k_{g}l}{2}},$

[0059] where G.sub.l and p.sub.l denote the flow of the natural gas and the pressure of the natural gas at the terminal end of the natural gas pipeline, respectively, G.sub.o and p.sub.o denote a flow of the natural gas and a pressure of the natural gas at a first end of the natural gas pipeline, respectively, and l denotes a length of the natural gas pipeline;

[0060] step 1-8 of defining a propagation coefficient of the natural gas pipeline as γ.sub.gc=Z.sub.gY.sub.g, and defining a characteristic impedance of the natural gas pipeline as Z.sub.gc=Z.sub.g/Y.sub.g;

[0061] step 1-9 of expressing, based on the two equations of the flow and the pressure at the terminal end of the natural gas pipeline in the step 1-7 and the two definitions in the step 1-8, the natural gas pipeline equation in a form of a linear two-port network:

[00021]$\left[\begin{array}{c}{p}_l\\ G_l\end{array}\right]=\left[\begin{array}{cc}A& B\\ C& D\end{array}\right]\left[\begin{array}{c}{p}_0\\ G_0\end{array}\right],$

[0062] where A, B, C and D denote network parameters whose values are:

[00022]$A=\left[\cosh \left(\frac{\sqrt{k_g^2+4{\gamma }_{\mathrm{gc}}}}{2}l\right)-\frac{k_g}{\sqrt{k_g^2+4{\gamma }_{\mathrm{gc}}}}\sinh (\frac{\sqrt{k_g^2+4{\gamma }_{\mathrm{gc}}}}{2}l\right].Math.e^{-\frac{k_{g}l}{2}},$$B=-\frac{2}{\sqrt{k_g^2+4/Z_{gc}}}\sinh \left(\frac{\sqrt{k_g^2+4{\gamma }_{gc}}}{2}l\right).Math.e^{-\frac{k_{g}l}{2}},$$C=-\frac{2}{\sqrt{k_g^2+4Z_{gc}}}\sinh \left(\frac{\sqrt{k_g^2+4{\gamma }_{gc}}}{2}l\right).Math.e^{-\frac{k_{g}l}{2}},\mathrm{and}$$D=\left[\cosh \left(\frac{\sqrt{k_g^2+4{\gamma }_{gc}}}{2}l\right)+\frac{k_g}{\sqrt{k_g^2+4{\gamma }_{gc}}}\sinh (\frac{\sqrt{k_g^2+4{\gamma }_{gc}}}{2}l\right].Math.e^{-\frac{k_{g}l}{2}},$

[0063] step 1-10 of establishing, based on the linear two-port network equation in the step 1-9, a π-type equivalent gas circuit as shown in FIG. 2, and equivalent parameters are:

Z=−B,

K=1−AD+BC,

Y.sub.1=(AD−BC−A)/B, and

Y.sub.1=(1−D)/B.

[0064] Step 2 of establishing a general branch model of the natural gas circuit, includes:

[0065] step 2-1 of establishing a mathematical model of a natural gas compressor as follows:

p.sub.1=p.sub.2+E.sub.g,

[0066] where p.sub.1 and p.sub.2 denote pressures on both sides of the natural gas compressor, and E.sub.g denotes a pressure increment provided by the natural gas compressor;

[0067] step 2-2 of forming, based on the gas resistance, the gas inductance, the gas capacitance, the controlled gas pressure source in the step 1-5 and the natural gas compressor in the step 2-1, a general branch (as shown in FIG. 3), a general branch equation being as follows:

G.sub.b=y.sub.b(p.sub.b+E.sub.b−k.sub.bp.sub.f),

[0068] where G.sub.b denotes a flow in the branch, G.sub.b is an unknown quantity, which is obtained by solving the natural gas circuit equation to obtain the pressure of each node in the natural gas circuit. p.sub.b is a gas pressure difference between two ends of the branch, p.sub.f is a pressure at a first end of the branch, when the first end of the branch is a gas source, p.sub.f is a known quantity, when the first end of the branch is not the gas source, p.sub.f is an unknown quantity, p.sub.t is a pressure at a terminal end of the branch, p.sub.t is an unknown quantity, y.sub.b is a branch admittance composed of the gas resistance, the gas inductance and the gas capacitance, and k.sub.b and E.sub.b denote a component parameter of a controlled gas pressure source in the branch and a component parameter of the natural gas compressor, which are provided by the natural gas circuit management party;

[0069] step 2-3 of writing branch equations of all branches in the natural gas circuit into a matrix form as follows:

G.sub.b=y.sub.b(p.sub.b+E.sub.b−k.sub.bp.sub.f),

[0070] where G.sub.b denotes a vector formed by the flow of each branch, y.sub.b denotes a diagonal matrix formed by each branch admittance, p.sub.b denotes a vector formed by the pressure difference at both ends of each branch, E.sub.b is a vector formed by the pressure increment of each natural gas compressor, k.sub.b is a vector formed by the parameter of each controlled gas pressure source, and p.sub.f is a vector formed by the gas pressure at the first end of each branch.

[0071] Step 3 of establishing a topological constraint equation of the natural gas circuit, includes:

[0072] step 3-1 of defining a node-branch correlation matrix A.sub.g in the natural gas circuit, where the matrix is a matrix of n rows and m columns, where n denotes a number of nodes) and m denotes a number of branches, and (A.sub.g).sub.i,j denotes an element of the i-th row and the j-th column, (A.sub.g).sub.i,j=0 denotes that a branch j is not connected to a node i, (A.sub.g).sub.i,j=1 denotes the branch j flowing out from the node i, and (A.sub.g).sub.i,j=−1 denotes the branch j flowing into the node i;

[0073] step 3-2 of defining a node-outflow branch correlation matrix A.sub.g+ in the natural gas circuit, where the matrix retains an non-negative element in the matrix A.sub.g, that is, as for (A.sub.g+).sub.i,j, when the branch j flows out from the node i, the element is 1, and when the branch j does not flow out from the node i, the element is 0;

[0074] step 3-3 of establishing a mass conservation equation of nodes of the natural gas circuit:

A.sub.gG.sub.b=G.sub.n

[0075] where G.sub.n denotes a column vector formed by a flow injection at each node, where a flow at a gas load node in the natural gas circuit is a known quantity, a flow at a gas source node is an unknown quantity, and a flow at a node of a non-gas load and a non-gas source is 0;

[0076] step 3-4 of establishing a gas pressure equation of nodes of the natural gas circuit:

A.sub.g.sup.Tp.sub.n=p.sub.b, and

A.sub.g+.sup.Tp.sub.n=p.sub.f,

[0077] where p.sub.n denotes a column vector formed by pressures at each node, where a pressure at the gas source node in the natural gas circuit is a known quantity, a pressure at the gas load node is an unknown quantity, and a pressure at the node of the non-gas load and the non-gas source is an unknown quantity.

[0078] Step 4 of establishing a natural gas circuit equation, includes:

[0079] step 4-1 of substituting the equations established in the steps 3-3 and 3-4 into the branch equation established in the step 2-3, and obtaining a natural gas circuit equation in an unreduced form as follows:

G.sub.n=A.sub.gy.sub.b(A.sub.g.sup.Tp.sub.n+E.sub.b−k.sub.bA.sub.g′.sup.Tp.sub.n),

[0080] step 4-2 of defining a generalized node admittance matrix Y′.sub.g and a generalized node injection vector G′.sub.n as follows:

Y′.sub.g=A.sub.gy.sub.bA.sub.g.sup.T−A.sub.gy.sub.bk.sub.bA.sub.g+.sup.T, and

G′.sub.n=G.sub.n−A.sub.gy.sub.bE.sub.b;

[0081] step 4-3 of substituting Y′.sub.g and G′.sub.n defined in the step 4-2 into the unreduced natural gas circuit equation in the step 4-1 and obtaining a following natural gas circuit model equation:

Y′.sub.gp.sub.n=G′.sub.n,

[0082] solving the above natural gas circuit model to obtain an unknown node pressure in the natural gas circuit, and obtaining an unknown branch flow by using the branch equation to implement the operation control of the integrated energy system.