HEAT SUPPLY NETWORK HYDRAULIC CIRCUIT MODELING METHOD FOR COMPREHENSIVE ENERGY SYSTEM SCHEDULING

Abstract

A heat supply network hydraulic circuit modeling method for a comprehensive energy system scheduling is provided. The hydraulic analysis model is unified with the power network model, and the connection between the hydraulic dynamic state and the hydraulic steady state is established. Based on the characteristic equations of thermal pipelines, flow control valves and compressors, this method abstracts hydraulic circuit element models such as hydraulic resistance, hydraulic inductance, and hydraulic pressure source, establishes hydraulic branch characteristics of the heat supply network based on the above hydraulic circuit elements, establishes the hydraulic topology constraints of the heat supply network based on Kirchhoff-like voltage and current laws, and establishes the steady hydraulic network equation by combining the above hydraulic branch characteristics and hydraulic topology constraints.

Claims

1. A heat supply network hydraulic circuit modeling method for a comprehensive energy system scheduling, comprising: step 1 of establishing a device model for a heat supply network, comprising: step 1-1 of establishing a thermal pipeline model in the heat supply network, comprising: step 1-1-1 of establishing a mass conservation equation and a momentum conservation equation describing a one-dimensional flow process of water in the thermal 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 ρ, ν and p denote a density, a flow rate and a pressure of the water, respectively, λ, D and θ denote a friction coefficient, an inner diameter and an inclination angle of the thermal pipeline, respectively, g denotes an acceleration of gravity, and t and denote time and space, respectively; step 1-1-2 of establishing, based on a fact that the water is an incompressible fluid, a differential equation of the density of the water about the time and the space: $\frac{\partial \rho }{\partial t}=\frac{\partial \rho }{\partial x}=0;$ step 1-1-3 of ignoring the convection term $\frac{\partial \rho v^2}{\partial x}$ in the momentum conservation equation in the step 1-1-1, that is $\frac{\partial \rho v^2}{\partial x}\approx 0,$ and performing an incremental linearization approximation on the square term of the flow rate in the resistance term $\frac{\lambda \rho v^2}{2D},$ that is, letting ν.sup.2≈2ν.sub.baseν−ν.sub.base.sup.2, where ν.sub.base denotes a base value of the flow rate of the water in the thermal pipeline, taking a flow rate in a design condition; step 1-1-4 of substituting the steps 1-1-2 and 1-1-3 into the step 1-1-1 to get following equations: $\frac{\partial G}{\partial x}=0,\mathrm{and}\frac{1}{A}\frac{\partial G}{\partial t}+\frac{\partial p}{\partial x}+\frac{\lambda G_{\mathrm{base}}}{\rho A^{2}D}G-\frac{\lambda G_{\mathrm{base}}^2}{2\rho A^{2}D}+\rho g\sin \theta =0,$ where G denotes a mass flow of the water, G=ρνA denotes a cross-sectional area of the thermal pipeline, and G.sub.base denotes a base value of the mass flow corresponding to the base value of the flow rate, that is G.sub.base=ρν.sub.baseA; step 1-1-5 of establishing, based on the step 1-1-4, equations of a flow difference and a pressure drop at both ends of a micro-element dx of the thermal pipeline: $\mathrm{dG}=0,\mathrm{and}\mathrm{dp}=-\frac{\mathrm{dx}}{A}.Math.\frac{\mathrm{dG}}{\mathrm{dt}}-\frac{\lambda G_{\mathrm{base}}\mathrm{dx}}{\rho A^{2}D}G-\left(\rho g\sin \theta -\frac{\lambda G_{\mathrm{base}}^2}{2\rho A^{2}D}\right)\mathrm{dx},$ where dG denotes the flow difference at both ends of the micro-element of the thermal pipeline, and dp denotes the pressure drop at both ends of the micro-element of the thermal pipeline; step 1-1-6 of obtaining, based on the equations of the flow difference and the pressure drop at both ends of the micro-element of the thermal pipeline in the step 1-1-5, calculation equations of a hydraulic resistance R.sub.h, a hydraulic inductance L.sub.h and a hydraulic pressure source E.sub.h in the thermal pipeline as follows:
R.sub.h=λG.sub.base/(ρA.sup.2D),
L.sub.h=1/A, and
E.sub.h=ρg sinθ−λG.sub.base.sup.2/(2ρA.sup.2D), wherein the micro-element dx of the thermal pipeline denotes a hydraulic circuit comprising 3 elements, and the entire thermal pipeline denotes a distributed parameter hydraulic circuit; step 1-1-7 of establishing, based on element parameters of the distributed parameter hydraulic circuit of the thermal pipeline in the step 1-1-6, element parameters of a lumped parameter hydraulic circuit of the thermal pipeline:
R=R.sub.hl,
L=L.sub.hl, and
E=E.sub.hl, where R denotes a hydraulic resistance in the lumped parameter hydraulic circuit of the thermal pipeline, L denotes a hydraulic inductance in the lumped parameter hydraulic circuit of the thermal pipeline, E denotes a hydraulic pressure source in the lumped parameter hydraulic circuit of the thermal pipeline, and l denotes a length of the thermal pipeline; step 1-1-8 of performing Fourier transform on an excitation of the lumped parameter hydraulic circuit of the thermal pipeline, decomposing the excitation of the lumped parameter hydraulic circuit of the thermal pipeline into a plurality of sinusoidal steady state excitations at different frequencies, and establishing algebraic equations of a lumped parameter frequency domain hydraulic circuit corresponding to each frequency component co in the sinusoidal steady state excitation:
p.sub.1=p.sub.0−(R+jωL)G.sub.0−E, and
G.sub.1=G.sub.0, where p.sub.0 and G.sub.0 denote a pressure and a flow at a head end of the thermal pipeline, respectively, and p, and G, denote a pressure and a flow at a terminal end of the thermal pipeline; step 1-2 of establishing a flow control valve model in the heat supply network, comprising: step 1-2-1 of establishing an equation between a pressure difference p on both sides of the flow control valve and a mass flow G of the flow control valve:
p=k.sub.νG.sup.2, where k, denotes an opening coefficient of the flow control valve, and G denotes the mass flow of the water; step 1-2-2 of performing an incremental linearization approximation on the square term G.sup.2 of the mass flow in the step 1-2-1, that is G.sup.2=2G.sub.baseG−G.sub.base.sup.2, and converting the equation between the pressure difference p on both sides of the flow control valve and the mass flow G of the flow control valve in the step 1-2-1 into the following equation:
p=2k.sub.νG.sub.base.Math.G−k.sub.νG.sub.base.sup.2, step 1-2-3 of defining, based on the step 1-2-2, calculation equations of a hydraulic resistance R.sub.ν and a hydraulic pressure source E.sub.ν of the flow control valve as follows:
R.sub.ν=2k.sub.νG.sub.base, and
E.sub.ν=−k.sub.νG.sub.base.sup.2; step 1-3 of establishing a compressor model in the heat supply network, comprising: step 1-3-1 of establishing an equation between a pressure difference p on both sides of the compressor and a mass flow G of the water of the compressor at a given rotation speed:
p=−(k.sub.p1G.sup.2+k.sub.p2ω.sub.pG+k.sub.p3ω.sub.p.sup.2), where k.sub.p1, k.sub.p2 and k.sub.p3 denote inherent coefficients of the compressor, which are obtained from a factory nameplate of the compressor or obtained through an external characteristic testing and fitting, and ω.sub.p denotes a rotation frequency of the compressor; step 1-3-2 of performing an incremental linearization approximation on the square term G.sup.2 of the mass flow in the step 1-3-1, that is G.sup.2=2G.sub.baseG−G.sub.base.sup.2, and converting the step 1-3-1 into following equation:
p=−(2k.sub.p1G.sub.base+k.sub.p2ω.sub.p).Math.G−(k.sub.p3ω.sub.p.sup.2−k.sub.p1G.sub.base.sup.2); step 1-3-3 of defining, based on the step 1-3-2, calculation equations of a hydraulic resistance R.sub.p and a hydraulic pressure source E.sub.p of the compressor as follows:
R.sub.p=−(2k.sub.p1G.sub.base+k.sub.p2ω.sub.p), and
E.sub.p=k.sub.p1G.sub.base.sup.2−k.sub.p3ω.sub.p.sup.2; step 2 of establishing a hydraulic branch characteristic of the heat supply network, comprising: step 2-1 of establishing, based on the models of the thermal pipeline, the flow control valve and the compressor established in the step 1, a hydraulic branch characteristic equation of the heat supply network:
G.sub.b=y.sub.b(p.sub.b−E.sub.b), where G.sub.b denotes a base value of a mass flow corresponding to a base value of a flow rate in a hydraulic branch, p.sub.b denotes a hydraulic pressure difference at both ends of the hydraulic branch, y.sub.b denotes a branch admittance formed by a hydraulic resistance and a hydraulic inductance in the hydraulic branch, and E.sub.b denotes a sum of hydraulic pressure sources in the hydraulic branch; step 2-2 of writing hydraulic branch equations of all the hydraulic branches in the heat supply network into a matrix form as follows:
G.sub.b=y.sub.b(p.sub.b−E.sub.b), where G.sub.b, p.sub.b and E.sub.b denote a column vector composed of the mass flow of the water in all the hydraulic branches in the heat supply network, a column vector composed of the hydraulic pressure difference at both ends of each branch in all the hydraulic branches in the heat supply network, and a column vector composed of the hydraulic pressure sources in all the hydraulic branches in the heat supply network, respectively, and y.sub.b denotes a diagonal matrix composed of the admittances of all the branches of the heat supply network; step 3 of establishing hydraulic topology constraints of the heat supply network, comprising: step 3-1 of defining a node-branch correlation matrix A.sub.h in the heat supply network, which is a matrix of n rows and m columns, where n denotes a number of nodes in the heat supply network, and m denotes a number of branches in the heat supply network, (A.sub.h).sub.i,j denotes an element in an i.sup.th row and a j.sup.th column in denotes that the branch j is not connected to the node i, (A.sub.h).sub.i,j=1 denotes that the branch j flows out from the node i , and (A.sub.h).sub.i,j=−1 denotes that the branch flows into the node i; step 3-2 of establishing, based on Kirchhoff-like current law, a mass conservation equation of nodes of the heat supply network:
A.sub.hG.sub.b=G.sub.n, where G.sub.n denotes a column vector formed by water injection of each node, and when the heat supply network is a closed network, then G.sub.n=0; step 3-3 of establishing, based on Kirchhoff-like voltage law, a loop pressure drop equation of the heat supply network:
A.sub.h.sup.Tp.sub.n=p.sub.b, where p.sub.n denotes a column vector composed of a hydraulic pressure value of each node; step 4 of establishing a dynamic hydraulic network equation of the heat supply network, comprising: step 4-1 of substituting the hydraulic topology constraints established in the steps 3-2 and 3-3 into the branch characteristic equation established in the step 2-2 and obtaining an unreduced form of a hydraulic network equation of the heat supply network as follows:
A.sub.hy.sub.bA.sub.b.sup.Tp.sub.n=G.sub.n+A.sub.hy.sub.bE.sub.b; step 4-2 of defining a generalized node admittance matrix Y.sub.h and a generalized node injection vector G′.sub.n as follows:
Y.sub.h=A.sub.hy.sub.bA.sub.h.sup.T, and
G′.sub.n=G.sub.n+A.sub.hy.sub.bE.sub.b; step 4-3 of substituting the Y.sub.h and G′.sub.n defined in the step 4-2 into the unreduced form of the hydraulic network equation of the heat supply network in the step 4-1, and obtaining the hydraulic network equation in the heat supply network of following form:
Y.sub.hp.sub.n=G′.sub.n, wherein the above hydraulic network equation describes a hydraulic dynamic of the heat supply network; step 5 of deleting a hydraulic inductance element in the hydraulic circuit model of the heat supply network, recalculating, based on the step 4, the generalized node admittance matrix Y.sub.h, and only taking a zero-frequency component in a frequency domain to degenerate the dynamic hydraulic circuit model into a steady hydraulic circuit model, wherein the steady hydraulic circuit model is the heat supply network hydraulic circuit model for the comprehensive energy system control.

Description

BRIEF DESCRIPTION OF DRAWINGS

[0063] FIGS. 1a and 1b show schematic diagrams of a distributed parameter hydraulic circuit of a thermal pipeline involved in a method of the present disclosure, wherein FIG. 1a is the distributed parameter hydraulic circuit of the entire thermal pipeline, and FIG. 1b is a distributed parameter hydraulic circuit of a micro-element with a dx length in the thermal pipeline.

[0064] FIG. 2 is a schematic diagram of a lumped parameter equivalent hydraulic circuit of a thermal pipeline.

DESCRIPTION OF EMBODIMENTS

[0065] A heat supply network hydraulic circuit modeling method for a comprehensive energy system control proposed in the present disclosure includes the following steps:

[0066] step 1 of establishing a device model for a heat supply network, including:

[0067] step 1-1 of establishing a thermal pipeline model in the heat supply network, including:

[0068] step 1-1-1 of establishing a mass conservation equation and a momentum conservation equation describing a one-dimensional flow process of water in the thermal pipeline:

[00008]$\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,$

[0069] where ρ, ν and p denote a density, a flow rate and a pressure of the water, respectively, λ, D and θ denote a friction coefficient, an inner diameter and an inclination angle of the thermal pipeline, respectively, g denotes an acceleration of gravity, and t and x denote time and space, respectively;

[0070] step 1-1-2 of establishing, based on a fact that the water is an incompressible fluid, a differential equation of the density of the water about the time and the space:

[00009]$\frac{\partial \rho }{\partial t}=\frac{\partial \rho }{\partial x}=0;$

[0071] step 1-1-3 of ignoring the convection term

[00010]$\frac{\partial \rho v^2}{\partial x}$

in the momentum conservation equation in the step 1-1-1, that is

[00011]$\frac{\partial \rho v^2}{\partial x}\approx 0,$

and performing an incremental linearization approximation on the square term of the flow rate in the resistance term

[00012]$\frac{\lambda \rho v^2}{2D},$

that is, letting ν.sup.2≈2ν.sub.baseν−ν.sub.base.sup.2, where ν.sub.base denotes a base value of the flow rate of the water in the thermal pipeline, taking a flow rate in a design condition;

[0072] step 1-1-4 of substituting the steps 1-1-2 and 1-1-3 into the step 1-1-1 to get following equations:

[00013]$\frac{\partial G}{\partial x}=0,\mathrm{and}\frac{1}{A}\frac{\partial G}{\partial t}+\frac{\partial p}{\partial x}+\frac{\lambda G_{\mathrm{base}}}{\rho A^{2}D}G-\frac{\lambda G_{\mathrm{base}}^2}{2\rho A^{2}D}+\rho g\sin \theta =0,$

[0073] where G denotes a mass flow of the water, G=ρνA, A denotes a cross-sectional area of the thermal pipeline, and G.sub.base denotes a base value of the mass flow corresponding to the base value of the flow rate, that is G.sub.base=ρν.sub.base.sup.A;

[0074] step 1-1-5 of establishing, based on the step 1-1-4, equations of a flow difference and a pressure drop at both ends of a micro-element dx of the thermal pipeline:

[00014]$\mathrm{dG}=0,\mathrm{and}\mathrm{dp}=-\frac{\mathrm{dx}}{A}.Math.\frac{\mathrm{dG}}{\mathrm{dt}}-\frac{\lambda G_{\mathrm{base}}\mathrm{dx}}{\rho A^{2}D}G-\left(\rho g\sin \theta -\frac{\lambda G_{\mathrm{base}}^2}{2\rho A^{2}D}\right)\mathrm{dx},$

[0075] where dG denotes the flow difference at both ends of the micro-element of the thermal pipeline, and dp denotes the pressure drop at both ends of the micro-element of the thermal pipeline;

[0076] step 1-1-6 of obtaining, based on the equations of the flow difference and the pressure drop at both ends of the micro-element of the thermal pipeline in the step 1-1-5, calculation equations of a hydraulic resistance R.sub.h, a hydraulic inductance L.sub.h and a hydraulic pressure source E.sub.h in the thermal pipeline as follows:

[0077] R.sub.h=λG.sub.base/(ρA.sup.2D)

[0078] L.sub.h=1/A, and

[0079] E.sub.h=ρg sinθ−λG.sub.base.sup.2/(2ρA.sup.2D),

[0080] where the micro-element of the thermal pipeline denotes a hydraulic circuit including 3 elements, and the entire thermal pipeline denotes a distributed parameter hydraulic circuit, and the distributed parameter hydraulic circuit of the entire thermal pipeline and the distributed parameter hydraulic circuit of the micro-element dx of the thermal pipeline are shown in FIGS. la and lb;

[0081] step 1-1-7 of establishing, based on element parameters of the distribution parameter hydraulic circuit of the thermal pipeline in the step 1-1-6, element parameters of a lumped parameter hydraulic circuit of the thermal pipeline, where the lumped parameter hydraulic circuit of the thermal pipeline is shown in FIG. 2:

[0082] R=R.sub.hl,

[0083] L=L.sub.hl, and

[0084] E=E.sub.hl,

[0085] where R denotes a hydraulic resistance in the lumped parameter hydraulic circuit of the thermal pipeline, L denotes a hydraulic inductance in the lumped parameter hydraulic circuit of the thermal pipeline, E denotes a hydraulic pressure source in the lumped parameter hydraulic circuit of the thermal pipeline, and l denotes a length of the thermal pipeline;

[0086] step 1-1-8 of performing Fourier transform on an excitation of the lumped parameter hydraulic circuit of the thermal pipeline, decomposing the excitation of the lumped parameter hydraulic circuit of the thermal pipeline into a plurality of sinusoidal steady state excitations at different frequencies, and establishing algebraic equations of a lumped parameter frequency domain hydraulic circuit corresponding to each frequency component ω in the sinusoidal steady state excitation:

[0087] p.sub.l=p.sub.0−(R+jωL)G.sub.0−E, and

[0088] G.sub.l=G.sub.0,

[0089] where p.sub.0 and G.sub.0 denote a pressure and a flow at a head end of the thermal pipeline, respectively, and p.sub.l and G.sub.l denote a pressure and a flow at a terminal end of the thermal pipeline;

[0090] step 1-2 of establishing a flow control valve model in the heat supply network, including:

[0091] step 1-2-1 of establishing an equation between a pressure difference P on both sides of the flow control valve and a mass flow G of the flow control valve:

[0092] p=k.sub.νG.sup.2,

[0093] where k.sub.ν denotes an opening coefficient of the flow control valve, and G denotes the mass flow of the water;

[0094] step 1-2-2 of performing an incremental linearization approximation on the square term G.sup.2 of the mass flow in the step 1-2-1, that is G.sup.2=2G.sub.baseG−G.sub.base.sup.2, and converting the equation between the pressure difference p on both sides of the flow control valve and the mass flow G of the flow control valve in the step 1-2-1 into the following equation:

p=2k.sub.νG.sub.base.Math.G−k.sub.νG.sub.base.sup.2,

[0095] step 1-2-3 of defining, based on the step 1-2-2, calculation equations of a hydraulic resistance R.sub.ν and a hydraulic pressure source E.sub.ν of the flow control valve as follows:

[0096] R.sub.ν=2k.sub.νG.sub.base, and

[0097] E.sub.ν=−k.sub.νG.sub.base.sup.2;

[0098] step 1-3 of establishing a compressor model in the heat supply network, including:

[0099] step 1-3-1 of establishing an equation between a pressure difference p on both sides of the compressor and a mass flow G of the water of the compressor at a given rotation speed:

p.sub.l=−(k.sub.p1G.sup.2+k.sub.p2ω.sub.pG+k.sub.p3ω.sub.p.sup.2),

[0100] where .sup.kP1 .sup.kP2 and .sup.kP3 denote inherent coefficients of the compressor, which are obtained from a factory nameplate of the compressor or obtained through an external characteristic testing and fitting, and ω.sub.p denotes a rotation frequency of the compressor;

[0101] step 1-3-2 of performing an incremental linearization approximation on the square term G.sup.2 of the mass flow in the step 1-3-1, that is G.sup.2=2G.sub.baseG−G.sub.base.sup.2, and converting the step 1-3-1 into following equation:

p=−(2k.sub.p1G.sub.base+k.sub.p2ω.sub.p).Math.G−(k.sub.p3ω.sub.p.sup.2−k.sub.plG.sub.base.sup.2);

[0102] step 1-3-3 of defining, based on the step 1-3-2, calculation equations of a hydraulic resistance .sup.Rp and a hydraulic pressure source E.sub.p of the compressor as follows:

R.sub.p=−(2k.sub.p1G.sub.base+k.sub.p2ω.sub.p), and

E.sub.p=k.sub.p1G.sub.base.sup.2−k.sub.p3ω.sub.p.sup.2;

[0103] step 2 of establishing a hydraulic branch characteristic of the heat supply network, including:

[0104] step 2-1 of establishing, based on the models of the thermal pipeline, the flow control valve and the compressor established in the step 1, a hydraulic branch characteristic equation of the heat supply network:

G.sub.b=y.sub.b(p.sub.b−E.sub.b),

[0105] where G.sub.b denotes a base value of a mass flow corresponding to a base value of a flow rate in a hydraulic branch, p.sub.b denotes a hydraulic pressure difference at both ends of the hydraulic branch, y.sub.b denotes a branch admittance formed by a hydraulic resistance and a hydraulic inductance in the hydraulic branch, and E.sub.b denotes a sum of hydraulic pressure sources in the hydraulic branch;

[0106] step 2-2 of writing hydraulic branch equations of all the hydraulic branches in the heat supply network into a matrix form as follows:

G.sub.b=y.sub.b(p.sub.b−E.sub.b),

[0107] where G.sub.b, p.sub.b and E.sub.b denote a column vector composed of the mass flow of the water in all the hydraulic branches in the heat supply network, a column vector composed of the hydraulic pressure difference at both ends of each branch in all the hydraulic branches in the heat supply network, and a column vector composed of the hydraulic pressure sources in all the hydraulic branches in the heat supply network, respectively, and y.sub.b denotes a diagonal matrix composed of the admittances of all the branches of the heat supply network;

[0108] step 3 of establishing hydraulic topology constraints of the heat supply network, including:

[0109] step 3-1 of defining a node-branch correlation matrix A.sub.h in the heat supply network, which is a matrix of n rows and m columns, where n denotes a number of nodes in the heat supply network, and m denotes a number of branches in the heat supply network, (A.sub.h).sub.i,j denotes an element in an i.sup.th row and a j.sup.th column in (A.sub.h).sub.i,j, (A.sub.h).sub.i,j=0 denotes that the branch j is not connected to the node i, (A.sub.h).sub.i,j=1 denotes that the branch flows out from the node i, and (A.sub.h).sub.i,j=−1 denotes that the branch flows into the node i;

[0110] step 3-2 of establishing, based on Kirchhoff-like current law, a mass conservation equation of nodes of the heat supply network:

A.sub.hG.sub.b=G.sub.n,

[0111] where G.sub.n denotes a column vector formed by water injection of each node, and when the heat supply network is a closed network, then G.sub.n=0;

[0112] step 3-3 of establishing, based on Kirchhoff-like voltage law, a loop pressure drop equation of the heat supply network:

A.sub.h.sup.Tp.sub.n=p.sub.b,

[0113] where p.sub.n denotes a column vector composed of a hydraulic pressure value of each node;

[0114] step 4 of establishing a dynamic hydraulic network equation of the heat supply network, including:

[0115] step 4-1 of substituting the hydraulic topology constraints established in the steps 3-2 and 3-3 into the branch characteristic equation established in the step 2-2 and obtaining an unreduced form of a hydraulic network equation of the heat supply network as follows:

A.sub.hy.sub.bA.sub.h.sup.Tp.sub.n=G.sub.n+A.sub.hy.sub.bE.sub.b;

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

[0117] Y.sub.h=A.sub.hy.sub.bA.sub.h.sup.T, and

[0118] G′.sub.n=G.sub.n+A.sub.hy.sub.bE.sub.b;

[0119] step 4-3 of substituting the Y.sub.h and G′.sub.n defined in the step 4-2 into the unreduced form of the hydraulic network equation of the heat supply network in the step 4-1, and obtaining the hydraulic network equation in the heat supply network of following form:

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

[0120] where the above hydraulic network equation describes a hydraulic dynamic of the heat supply network, and has a unified form with the network equation of the power network;

[0121] step 5 of deleting a hydraulic inductance element in the hydraulic circuit model of the heat supply network, recalculating, based on the step 4, the generalized node admittance matrix Y.sub.h , and only taking a zero-frequency component in a frequency domain to degenerate the dynamic hydraulic circuit model into a steady hydraulic circuit model. When the transient hydraulic dynamic process does not need to be considered, the steady hydraulic circuit model can be used as the heat supply network hydraulic circuit model for the comprehensive energy system control.