Gradient descent iterative method for energy flow of integrated energy system considering multiple modes of compressors

Energy Conversion and Management 207 (2020) 112534

Contents lists available at ScienceDirect

Energy Conversion and Management journal homepage: www.elsevier.com/locate/enconman

Gradient descent iterative method for energy flow of integrated energy system considering multiple modes of compressors

T

⁎

Jinghua Li , Yujin Huang, Mengshu Zhu School of Electrical Engineering, Guangxi University, Nanning 530004, China

A R T I C LE I N FO

A B S T R A C T

Keywords: Integrated energy system Energy flow Compressor Operating modes Gradient descent

Energy flow plays an important role in optimal operation and planning analysis of integrated energy systems. However, in the energy flow calculation of integrated energy systems, the model of natural gas networks is inaccuracy. In traditional models, the gas compressor is modeled with a single compression ratio, without regard for the impact of different operating modes on energy flow. However, compressors often operate in different modes, so the constant compression ratio model will lead to inaccurate results of energy flow. To enhance the accuracy of energy flow, multiple compressor operating modes are considered in this paper, including constant compression ratio mode, constant outlet pressure mode, and constant natural gas flow mode. Also, for conciseness, the three modes are unified formulated by a general equation. Then, to ascertain the energy flow, a complete gradient descent iterative method is proposed, in which the partial derivatives of each state variable are calculated, and a new Jacobian matrix is derived. Compared with the non-gradient descent method, the computational speed of the proposed method is faster while the computational accuracy has not been reduced. The results tested on 14 node electricity-14 node heat-10 node gas and 1047 node electricity-32 node heat-48 node gas integrated systems show that the operating modes of compressors have great influences on the gas pressure distribution, which further affect the distributions of electricity and heat. Also, the superiority of the proposed gradient descent method is verified.

1. Introduction More and more attention has been paid to the technology for energy efficient utilization [1]. Integrated energy systems (IESs) can coordinate various energy carriers and improve energy efficiency [2], so the IES is becoming an important means to reverse the trends of fossil energy consumption and climate change [3]. Energy flow is the basis of optimal operation and planning analysis of the IES [4]. Based on the energy flow, the voltage, the temperature, and the gas pressure can be monitored and analyzed [5]. Hence, energy flow calculation has always been an active research area. At present, an IES energy flow model typically contains four parts: electricity network, heat network, natural gas network, and energy conversion of coupling components. The research on mathematical models of electricity networks [6] and heat networks [7] is relatively mature. The models of natural gas networks are still evolving. In particular, it has been shown that the gas compressor model is inadequate [8]. Liu and Mancarella [9] established steady-state models of electricity, heat and natural gas networks, and used a multi-energy efficiency matrix to model conversion components. However, compressors

⁎

of natural gas networks were not considered in [9]. Similarly, Sun and Chen [10] used the efficiency model to formulate the coupling of multiple energy networks, but still ignored compressors in the natural gas network. Further, the compressor was formulated as a constant compression ratio model in [11]. Using the constant compression ratio model, Martinez-Mares and Fuerte-Esquivel [11] considered the change of gas temperature to analyze the effect of temperature on the energy flow of integrated electricity-gas energy systems. In [12], based on the constant compression ratio model, the bi-directional energy conversion of the integrated electricity-gas energy system was considered to study the influence of wind power on the outputs of gas generators and power-to-gas units. In the integrated electricity-heat-gas energy system (IEHGES), Shabanpour-Haghighi and Seifi [13] considered compressors with constant compression ratio and proposed a model with nonlinear part-load efficiency performance for gas generators, combined heat and power (CHP) units and boilers. In order to further analyze the coupling among electricity, heat, and natural gas networks, Massrur et al. [14] used an energy hub to simulate the energy flow based on the constant compression ratio model. Compressors play an important role in the natural gas network,

Corresponding author at: School of Electrical Engineering, Guangxi University, No. 100 Daxue Road, Nanning 530004, China. E-mail address: [email protected] (J. Li).

https://doi.org/10.1016/j.enconman.2020.112534 Received 25 November 2019; Received in revised form 13 January 2020; Accepted 23 January 2020 0196-8904/ © 2020 Elsevier Ltd. All rights reserved.

Energy Conversion and Management 207 (2020) 112534

J. Li, et al.

Nomenclature

fCHP

Indices

fGB

i, j k

node i, j pipeline k

fG rc Bgu

Variables natural gas flow within pipeline k (m3/h) square of pressure at node i and j (bar2) vector of natural gas flow within each pipeline (m3/h) vector of natural gas flow through each compressor (m3/ h) Πc,in square of inlet pressure at compressor Πc,out square of outlet pressure at compressor ΔΠ vector of difference in square of pressure fu fu = [fc; fG], vector of natural gas flow through each nonpipeline element (m3/h) ΠI vector of square of pressure at non-outlet node ΠO vector of square of pressure at outlet node |Vi|, |Vj| voltage magnitude at node i and j (p.u.) θij voltage angle difference between node i and j (degree) PCP electrical power consumed by circulation pump (MW) Pc electrical power consumed by compressor (MW) m vector of mass flow within each pipe (kg/s) mq vector of injected mass flow at each node (kg/s) hf vector of head loss along each pipe (m) Ts,load vector of supply temperature at load node (°C) Ts,source vector of supply temperature at source node (°C) Tr,source vector of return temperature at source node (°C) Tstart temperature at starting node of heat pipe (°C) Tend temperature at ending node of heat pipe (°C) min mass flow coming into heat node (kg/s) mout mass flow leaving heat node (kg/s) Tin temperature of mass flow before mixing (°C) Tout mixture temperature (°C) mCP mass flow through circulation pump (kg/s) pc,in inlet pressure at compressor (bar) pc,out outlet pressure at compressor (bar) V vector of node voltage Tsn Tsn = Ts − Tam Trn Trn = Tr − Tam fk Πi, Πj f fc

fSL Gij Bij PS QS PG PGG PCHP PD PHP PEB Ah Kh Bh ΦD ΦG Cp To ΦCHP ΦHP ΦEB ΦGB Tam λ Lh ηGG Hgas ηHP ηEB ηGB ηCP ga hCP ηc α Y Ah,load Ah,source Cs bs Cr br

Parameters γ Kg,k Ag Bg fS fD fGG

exponent of natural gas flow resistance coefficient of gas pipeline k node-branch incidence matrix of natural gas network incidence matrix between gas node and compressor vector of natural gas flow at each node (m3/h) vector of natural gas flow consumed by load (m3/h) vector of natural gas flow consumed by gas-fired generator (m3/h)

vector of natural gas flow consumed by combined heat and power unit (m3/h) vector of natural gas flow consumed by gas-fired boiler (m3/h) vector of natural gas flow supplied by natural gas station (m3/h) constant compression ratio incidence matrix between gas node and non-pipeline element vector of natural gas consumption at each node (m3/h) conductance between node i and j susceptance between node i and j injected active power (MW) injected reactive power (MW) electrical power generated by coal-fired unit (MW) electrical power generated by gas-fired generator (MW) electrical power generated by combined heat and power unit (MW) electrical power consumed by load (MW) electrical power consumed by heat pump (MW) electrical power consumed by electric boiler (MW) node-branch incidence matrix of heat network vector of pipe resistance coefficient of heat network loop-branch correlation matrix of heat network vector of heat power consumed by load (MW) vector of heat power supplied by source (MW) specific heat of water (MJ/(kg·°C)) vector of outlet temperature at load node (°C) vector of heat power generated by combined heat and power unit (MW) vector of heat power generated by heat pump (MW) vector of heat power generated by electric boiler (MW) vector of heat power generated by gas-fired boiler (MW) ambient temperature (°C) heat transfer coefficient of heat pipe (MW/(m·°C)) length of heat pipe (m) efficiency factor of gas-fired generator calorific value of natural gas (MJ/m3) efficiency of heat pump efficiency of electric boiler efficiency of gas-fired boiler efficiency of circulation pump gravity acceleration (m/s2) pump head of heat network (m) efficiency of compressor polytropic exponent of compressor admittance matrix the part corresponding to load node in matrix Ah the part corresponding to source node in matrix Ah coefficient matrix for supply temperature calculation solution vector for supply temperature calculation coefficient matrix for return temperature calculation solution vector for return temperature calculation

compressors to operate in different modes during dispatch periods. There are three operating modes: constant compression ratio mode, constant outlet pressure mode, and constant natural gas flow mode. These operating modes will directly change the pressure distribution, thus affecting the energy flow of the natural gas network [16]. Therefore, in order to obtain more accurate energy flow results, it is necessary to consider the diversity of compressor operating modes. To solve the energy flow model of the IES, two solution frameworks are currently available: the integrated framework and the decomposed

which compensate for the pressure loss of natural gas in the transmission process [15]. Although compressors have been considered in the above research, the different operating modes of compressors have been neglected. It is usually assumed that all compressors operate in the constant compression ratio mode. In this case, the conventional model with a constant compression ratio will lead to inaccurate energy flow results when compressors operate in different modes. Generally, the operating modes of compressors are set according to the system requirements to maintain the pressure level, so it is common for 2

Energy Conversion and Management 207 (2020) 112534

J. Li, et al.

framework. In the both frameworks, the gradient descent iterative method has been widely used. The gradient descent iterative method is effective for solving nonlinear equations, and the key of the method is the derivation of the Jacobian matrix [17]. In the integrated framework, multiple energy network models are solved as a whole. In the integrated electricity-gas energy system, Martinez-Mares and FuerteEsquivel [11] calculated the Jacobian matrix of the natural gas pipeline flow equation. Further, Zeng et al. [12] recalculated the Jacobian matrix considering power-to-gas units. In the integrated electricity-heat energy system, Liu et al. [18] derived the Jacobian matrix of the heat power balance equation, the loop pressure equation, the supply temperature equation, and the return temperature equation. In the IEHGES, Liu and Mancarella [9] calculated the Jacobian matrix considering gas generators, CHP units, heat pumps (HPs), and gas boilers. In addition, Shabanpour-Haghighi and Seifi [13] derived the Jacobian matrix considering the nonlinear part-load efficiency performance of units. In the decomposed framework, each energy network model is solved separately in each iteration. Based on the decomposed framework, Liu [19] solved both electricity and heat network models using the gradient descent iterative method, and calculated the Jacobian matrix of energy conversion equations when CHP units and HPs were considered as coupling components. Further, Masrur et al. [14] solved the electricity network model, the heat network model, and the natural gas network model using the holomorphic embedding method, the graph theory method, and the gradient descent iterative method, respectively. In the above research, the gradient descent iterative method has good applications in solving the energy flow problem of the IES. Although the Jacobian matrix of natural gas pipeline flow equations and energy conversion equations has been studied, these applications do not consider the equations for characterizing the diversity of compressor operating modes in natural gas networks. Osiadacz [16] proposed an iterative method based on the coefficient matrix of natural gas network equations for considering the diversity of compressor operating modes. However, because the elements of the iterative matrix of the method are not derivatives of the equations, the convergence speed of the method needs to be further improved. In the real-time monitoring of operation status of the IES, the slow computational speed will lead to the problem that the unsafe operation status cannot be alerted in time, such as the voltage, the temperature, and the pressure exceeding the limits. Hence, it is essential to enhance the computational speed of energy flow. In this paper, the diversity of compressor operating modes is first considered in the energy flow model of the IEHGES for more accurate results, and the effect of the compressor operating mode diversity on the energy flow results is analyzed. Subsequently, in order to improve the computational speed of energy flow, a complete gradient descent iterative method is proposed to solve the energy flow model based on the integrated framework. In the proposed method, the partial derivatives of each state variable are calculated and the Jacobian matrix is derived. The main contributions of this paper are as follows:

Fig. 1. Schematic diagram of the IEHGES.

2. Mathematical formulation of integrated electricity-heat-gas energy systems In the IEHGES, depicted in Fig. 1, electricity, heat, and natural gas networks are interconnected by coupling components such as gas-fired generators (GGs), CHP units, HPs, electric boilers (EBs), gas-fired boilers (GBs), circulation pumps (CPs), and electric compressors. The mathematical models of each network and coupling components are described below. 2.1. Natural gas network model In the natural gas network, pipelines, compressors, sources and loads are considered. The natural gas network model is composed of the pipeline flow equation, the nodal natural gas balance equation, and the compressor equation, as shown in Sections 2.1.1–2.1.4. 2.1.1. Pipeline flow equation There is a non-linear relationship between the natural gas flow through a pipeline and the pressure square difference at both ends of the pipeline, as follows: 1

(a) The diversity of compressor operating modes is considered in the energy flow model of the IEHGES, and the influence of multiple compressor operating modes on energy flow is analyzed. (b) The Jacobian matrix of energy flow equations considering the diversity of compressor operating modes is derived, and an improved iterative method with complete gradient descent is proposed to improve the computational speed of energy flow.

Sij (Πi − Πj ) ⎞ γ fk = φ (Πi − Πj ) = Sij ⎜⎛ ⎟ K g, k ⎝ ⎠

(1)

+ 1, Πi − Πj ⩾ 0 Sij = ⎧ − ⎨ ⎩ 1, Πi − Πj < 0

(2)

where φ is the function of pipeline flow; and Sij indicates the direction of natural gas flow. 2.1.2. Nodal natural gas balance equation The natural gas flow balance equation means that the sum of the flow through pipelines connected with the node, the flow through compressors connected with the node, and the flow consumption of the node is zero, as shown in Eq. (3). The flow consumption includes the consumption of the load, the GG, the CHP unit, and the GB, and the negative consumption of the natural gas station, as shown in Eq. (4).

The remainder of this paper is organized as follows. Section 2 introduces the mathematical formulation of multiple energy networks and coupling components. Section 3 presents the proposed integrated energy flow solution. Section 4 outlines case studies conducted. Section 5 presents concluding remarks.

3

Ag f + Bg fc + fS = 0

(3)

fS = fD + fGG + fCHP + fGB − fG

(4)

Energy Conversion and Management 207 (2020) 112534

J. Li, et al.

Ah m = m q

2.1.3. Compressor equation Traditionally, the compressor model with constant compression ratio is defined by Eq. (5):

Πc,out = rc2 Πc,in

2.3.2. Loop pressure equation The pipe friction causes head losses, and the loop pressure equation indicates that the sum of the head losses in a closed loop is zero.

(5)

On the other hand, a general equation is used to formulate the diversity of compressor operating modes:

w1 Πc,in + w2 Πc,out + w3 fc = d

Ag φ (ΔΠ ) + Bgu fu + fSL = 0

(7)

fSL = fD + fGG + fCHP + fGB

(8)

W1 ΠI + W2 ΠO + W3 fu − d = 0

(9)

Bh h f = 0

(15)

Φ D = Cp m q (Ts,load − To)

(16)

ΦG = Cp m q (Ts,source − Tr,source)

(17)

ΦG = ΦCHP + Φ HP + Φ EB + ΦGB

(18)

2.3.4. Temperature drop equation Due to energy loss, the temperature of the mass flow will decrease after flowing through the pipe, and the ratio of the temperature change is exponential with the length of the pipe, as shown in Eq. (19).

λLh ⎞ Tend = (Tstart − Tam ) exp ⎜⎛− ⎟ + Tam ⎝ Cp m ⎠

(19)

2.3.5. Mixture temperature equation The ideal mixing process for water is presented in Eq. (20). In this process, the total heat energy of water remains constant, that is, the sum of the product of mass flow and temperature is constant.

2.2. Electricity network model The electricity network model is composed of the nodal active power balance equation and the nodal reactive power balance equation. In Eq. (10), the active power balance equation shows that the sum of the active power of lines connected with the node is equal to the injected active power of the node. As shown in Eq. (11), the injected active power includes the active power provided by the coal-fired unit, the GG and the CHP unit, as well as the active power consumed by the load, the HP, the EB, the CP and the compressor. In Eq. (12), the reactive power balance equation shows that the sum of the reactive power of lines connected with the node is equal to the injected reactive power of the node.

( ∑ mout ) Tout =

∑ (min Tin)

(20)

2.4. Coupling component model GGs consume natural gas and generate electricity [20]:

fGG =

3600PGG ηGG Hgas

(21)

CHP units can produce both heat and electricity [21]. There are two types of turbines: back-pressure steam turbine (CHPI) and extraction steam turbine (CHPII) [22]. The mathematical formula defining the CHPI is given by Eqs. (22) and (23), and that for the CHPII by Eq. (24):

|Vi | ∑ |Vj |(Gij cos θij + Bij sin θij ) − PS, i = 0 (10) (11)

|Vi | ∑ |Vj |(Gij sin θij − Bij cos θij ) − QS, i = 0 j∈i

(14)

2.3.3. Nodal heat power balance equation The heat power consumed by the load is proportional to the product of mass flow and the difference between the load supply temperature and the load outlet temperature, as shown in Eq. (16). Similarly, the heat power supplied by the source is proportional to the product of mass flow and the difference between the source supply temperature and the source return temperature, as shown in Eq. (17). As shown in Eq. (18), heat sources include the CHP unit, the HP, the EB and the GB.

2.1.4. Natural gas network equation After extending the compressor to non-pipeline elements and using Eq. (6) to model the non-pipeline elements including compressors and gas sources, the natural gas network model can be obtained as follows:

PS, i = (PG, i + PGG, i + PCHP, i ) − (PD, i + PHP, i + PEB, i + PCP, i + Pc, i )

h f = Kh m|m|

(6)

where w1, w2, w3, and d are coefficients that characterize the operating modes of a compressor. There are usually three operating modes for a compressor. Mode I: the compressor operates at a constant compression ratio rc; thus, w1 = rc2, w2 = −1, w3 = d = 0. Mode II: the compressor operates at a preset outlet pressure whose square is Πc,set; thus, w1 = w3 = 0, w2 = 1, d = Πc,set. Mode III: the compressor operates at a preset natural gas flow fc,set; thus, w1 = w2 = 0, w3 = 1, d = fc,set. In this model, when the operating mode of the compressor changes, only the parameters of the general equation need to be updated, instead of changing the compressor equation and recalculating the Jacobian matrix.

j∈i

(13)

z CHPI =

ΦCHPI PCHPI

(22)

fCHPI =

3600PCHPI ηCHPI Hgas

(23)

(12)

2.3. Heat network model

ZCHPII = In the heat network model, the mass flow balance equation, the loop pressure equation, the nodal heat power balance equation, the temperature drop equation and the mixture temperature equation are considered.

ΦCHPII ηCHPII fCHPII Hgas/3600 − PCHPII

(24)

where ΦCHPI, PCHPI, zCHPI, and ηCHPI are the heat power output (MW), the electrical power output (MW), the heat-to-power ratio, and the electrical efficiency of the CHPI, respectively; fCHPI is the natural gas flow (m3/h) consumed by the CHPI; ΦCHPII and PCHPII are the heat power output (MW) and the electrical power output (MW) of the CHPII, respectively; fCHPII is the natural gas flow (m3/h) consumed by the CHPII; ZCHPII is the ratio describing the trade-off between electrical power and heat power of the CHPII; and ηCHPII is the electrical

2.3.1. Mass flow balance equation The mass flow balance equation describes the continuity of mass flow, that is, the mass flow entering a node is equal to the flow consumption at the node plus the mass flow leaving the node. 4

Energy Conversion and Management 207 (2020) 112534

J. Li, et al.

efficiency of the CHPII in full condensing mode. HPs [23] and EBs [24] convert electricity into heat:

Φ HP = PHP ηHP

(25)

Φ EB = PEB ηEB

(26)

ΠI ⎤ f ⎡ R1 R2 Bgu,I ⎤ ⎡ ΠO ⎥ + ⎡ SL,I ⎤ = ⎡ 0 ⎤ ⎢ R T R3 Bgu,O ⎥ ⎢ ⎢ f ⎥ ⎣0⎦ ⎢ ⎥ 2 ⎣ ⎦ ⎣ fu ⎦ ⎣ SL,O ⎦

where R1 represents the connection between pipelines; R2 represents the connection between pipelines and non-pipeline elements; R3 represents the connection between non-pipeline elements; Bgu,I is the non-outlet part of Bgu; Bgu,O is the outlet part of Bgu; fSL,I is the nonoutlet part of fSL; and fSL,O is the outlet part of fSL. Combining Eq. (9), Eq. (36) can be rewritten as follows:

GBs burn natural gas to generate heat [25]:

fGB =

3600ΦGB ηGB Hgas

(27)

In the heat network, CPs consume electricity to provide hydraulic pressure to ensure that the water returns to the heat source [26]:

PCP =

⎡ ΔfSL,I ⎤ ⎡ R1 R2 Bgu,I ⎤ ⎡ ΠI ⎤ ⎡ fSL,I ⎤ ⎡ 0 ⎤ ⎥= 0 ⎢ Δf ⎥ = ⎢ R2T R3 Bgu,O ⎥ ⎢ ΠO ⎥ + ⎢ f ⎥ ⎢ ⎥ ⎢ SL,O ⎥ ⎢ 0 ⎥ ⎢ SL,O ⎥ ⎢ f ⎣ Δd ⎦ ⎣ W1 W2 W3 ⎦ ⎣ u ⎦ ⎢ ⎦ ⎣ ⎦ ⎣ −d ⎥

mCP ga hCP 106ηCP

(28)

⎡⎛ pc,out ⎞ ⎢⎜ ⎟ 3600ηc (α − 1) ⎢⎝ pc,in ⎠ ⎣ 0.1pc,in fc α

α−1 α

⎤ − 1⎥ ⎥ ⎦

3.2. Basic framework of iterative solution (29) Eqs. (30), (31), and (37) are solved simultaneously as a whole by iteration. The iterative form is as follows:

3. Gradient descent iterative method for energy flow calculation

x (t + 1) = x (t ) − (J (t ) )−1ΔF (t )

Based on Eqs. (1)–(29), the integrated model can be expressed by Eqs. (30)–(32).

⎡ Ag φ (ΔΠ ) + Bgu fu + fSL = 0 ⎤ Fg (x e , xh , x g ) = ⎢ fSL = fD + fGG + fCHP + fGB ⎥ ⎢ ⎥ ⎢ ⎣W1 ΠI + W2 ΠO + W3 fu − d = 0 ⎥ ⎦

⎤ 0⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

3.3. Calculation of iterative matrix Two iterative matrices are developed to deal with the energy flow model. 3.3.1. Iterative matrix with non-gradient descent direction Eqs. (39) and (40) define the iterative matrix with the non-gradient descent direction (NGDD). ∂F

⎡ ∂xee ⎢ ∂F J1 = ⎢ ∂x eh ⎢ ⎢ ∂Fg ⎢ ∂x e ⎣

(31)

(32)

In order to establish a unified matrix model of natural gas networks considering the diversity of compressor operating modes, Eq. (7) can be rewritten as (33)

R = Ag Ψ −1AgT

(34)

∂F h ∂x h ∂Fg ∂x h

∂Fe ∂xg

⎤ ⎥ ∂F h ⎥ ∂xg ⎥ ⎥ J1gg ⎥ ⎦

(39)

(40)

In Eq. (39), the formulations of electricity network elements and heat network elements are as presented by Liu [19]. J1gg is equal to the coefficient matrix of the natural gas network equation shown in Eq. (37). However, because the iterative direction based on J1 is not the direction of gradient descent, the convergence speed needs to be further improved when using the NGDD method to solve the energy flow of the IEHGES.

3.1. Model rewriting of natural gas network

RΠ + Bgu fu + fSL = 0

∂Fe ∂x h

⎡ R1 R2 Bgu,I ⎤ J1gg = ⎢ R2T R3 Bgu,O ⎥ ⎢ ⎥ ⎣ W1 W2 W3 ⎦

where Fe, Fh, and Fg represent equations of the electricity network, the heat network and the natural gas network, respectively; xe represents the variables of electricity network equations, including voltage angle vector θ and voltage magnitude vector |V|; xh represents the variables of heat network equations, including mass flow vector m, supply temperature vector Ts, and return temperature vector Tr; and xg represents the variables of natural gas network equations, including pressure variable vector Π and natural gas flow vector fu.

3.3.2. Iterative matrix with gradient descent direction In order to improve the convergence speed and reduce the computational time, the partial derivatives of each variable in Eq. (37) are obtained, and an iterative matrix with the gradient descent direction (GDD) is derived. ∂F

⎡ ∂xee ⎢ ∂F h J2 = ⎢ ∂x e ⎢ ⎢ ∂Fg ⎢ ∂x e ⎣

where Ψ is the diagonal matrix

Ψk = K g, k |fk |γ − 1

(38)

where x is the state variables, comprising xe, xh, and xg; ΔF is the mismatches of Fe, Fh, and Fg; t is the number of iterations; and J is the iterative matrix.

Real(V (YV )∗) − PS = 0 ⎤ ⎡ Fe (x e , xh , x g ) = ⎢ PS = PG + PGG + PCHP − PD − PHP − PEB − PCP − Pc ⎥ ⎥ ⎢ Imag(V (YV )∗) − QS = 0 ⎦ ⎣ (30)

Cp Ah,load m (Ts,load − To) − Φ D = 0 ⎡ ⎢Cp Ah,source m (Ts,source − Tr,source) − ΦG = ⎢ ΦG = ΦCHP + Φ HP + Φ EB + ΦGB Fh (x e , xh , x g ) = ⎢ ⎢ Bh Kh m|m| = 0 ⎢ Cs Tsn − bs = 0 ⎢ ⎢ C r Trn − br = 0 ⎣

(37)

where ΔfSL,I, ΔfSL,O, and Δd are the mismatches of natural gas equations; and [fSL,I; fSL,O] = fSL, and fSL is explained in Eq. (8). Thus, the model of natural gas networks is converted from Eq. (32) to Eq. (37).

The electrical power consumed by a compressor is as follows [27]:

Pc =

(36)

(35)

Eq. (33) can be written as 5

∂Fe ∂x h ∂F h ∂x h ∂Fg ∂x h

∂Fe ∂xg

∂F

⎤ ⎡ ∂xe ⎥ ⎢ e ∂F h ⎥ ⎢ ∂F h ∂xg ⎥ = ⎢ ∂x e ∂Fg ⎥ ⎢ ∂Fg ⎢ ∂xg ⎥ ⎦ ⎣ ∂x e

∂Fe ∂x h ∂F h ∂x h ∂Fg ∂x h

∂Fe ∂xg

⎤ ⎥ ⎥ ⎥ ⎥ J2gg ⎥ ⎦ ∂F h ∂xg

(41)

Energy Conversion and Management 207 (2020) 112534

J. Li, et al.

∂Δf

J2gg

⎡ SL,I ⎢ ∂ΠI ∂ΔfSL,O =⎢ ⎢ ∂ΠI ⎢ ∂Δd ⎢ ∂ΠI ⎣

∂ΔfSL,I

∂ΔfSL,I

∂ΠO ∂ΔfSL,O

∂fu

⎤ ⎥

⎡ =⎢ ⎥ ⎢ ⎥ ⎣ ⎥ ⎦

∂ΔfSL,O ⎥

∂ΠO

∂fu

∂Δd ∂ΠO

∂Δd ∂fu

∂ΔfSL

∂ΔfSL

∂Π

∂fu

∂Δd ∂Π

∂Δd ∂fu

⎤ ⎥ ⎥ ⎦ (42)

where ΔfSL = [ΔfSL,I; ΔfSL,O], Π = [ΠI; ΠO]. Firstly, ∂ΔfSL/∂Π can be written in the following form:

∂ΔfSL ∂Π

= Ag ΓAgT

(43)

where Γ is the diagonal matrix

Γk =

fk 1 γ Πi − Πj

(44)

Then, ∂ΔfSL/∂fu is equal to

∂ΔfSL ∂fu

Bgu,I ⎤ =⎡ = Bgu,IO ⎢ Bgu,O ⎥ ⎣ ⎦

(45) Fig. 3. Structure diagram in Case 1.

Finally, ∂Δd/∂Π and ∂Δd/∂fu are obtained

∂Δd = [W1 W2 ] = W4 ∂Π

(46)

∂Δd = W3 ∂fu

(47)

Table 1 Parameters of coupling components in Case 1.

Therefore, J2gg can be written as

Ag ΓAgT Bgu,IO ⎤ J2gg = ⎡ ⎢ W W3 ⎥ 4 ⎣ ⎦

(48)

Components

Parameters

GG: calorific value (Hgas), efficiency factor (ηGG) CHPI: heat-to-power ratio (zCHPI), electrical efficiency (ηCHPI) CHPII: ratio (ZCHPII), electrical efficiency (ηCHPII) CP: gravity acceleration (ga), efficiency (ηCP) Compressor: polytropic exponent (α), efficiency (ηc)

39 MJ/m3, 0.55 1.25, 0.35 4, 0.25 9.80 m/s2, 0.65 1.27, 0.8

3.4. Solution steps for energy flow of integrated electricity-heat-gas energy systems

conducted via MATLAB R2017a on an Intel Core i5, 2.90 GHz CPU personal computer with Windows 10.

The flow chart of the energy flow solution is shown in Fig. 2, and the main steps are as follows:

4.1. Case 1: the 14 node electricity-14 node heat-10 node gas integrated system In this section, the system structure and simulation parameters of the 14 node electricity-14 node heat-10 node gas integrated system are introduced, and the results of computational accuracy, compressor operating mode diversity, and computational speed are analyzed.

Step 1: Initialization of state variables. In the electricity network, voltage magnitude and voltage angle are set to one per unit and 0 degree, respectively. In the heat network, the initial values for supply temperature and return temperature are selected with the temperature of heat sources and the outlet temperature of heat loads, respectively. In the natural gas network, the nodal pressures at the ends of each pipeline are considered to have a difference of 5–10% between the receiving node and the sending node [13]. Step 2: Calculate the mismatches of Fe, Fh, and Fg according to Eqs. (30), (31) and (37). If the maximum mismatch is greater than or equal to the mismatch tolerance, go to Step 3; if the maximum mismatch is less than the mismatch tolerance, go to Step 4. Step 3: Calculate the Jacobian matrix according to Eq. (41) and update the state variables according to Eq. (38). Then, go to Step 2. Step 4: Output the energy flow results. 4. Case studies

4.1.1. Description of Case 1 The first IES consists of a 14-node IEEE test system, a 14-node heat network [26], and a 10-node natural gas network [16]. The heat network contains 2 heat sources and 13 pipes. Moreover, the natural gas network includes 1 gas source and 10 pipelines. Fig. 3 shows the connection of networks and coupling components. The system contains 1 GG, 1 CHPI unit, 1 CHPII unit, 2 CPs, and 2 compressors. The parameters of the coupling components are summarized in Table 1. The variables to be solved include 14 voltage magnitudes |V|, 14 voltage angles θ, mass flow m of 13 pipes, 14 supply temperatures Ts, 14 return temperatures Tr, 10 pressure variables Π, and natural gas flow fu of 3 non-pipeline elements.

Two simulated IESs are adopted as test systems to analyze the influence of compressor operating modes on energy flow results and verify the effectiveness of the GDD method. The simulations are

4.1.2. Comparisons of computational accuracy To verify the computational accuracy of the GDD method in solving the energy flow of the IEHGES, the energy flow results of GDD and

Fig. 2. Flow chart of energy flow solution. 6

Energy Conversion and Management 207 (2020) 112534

J. Li, et al.

Table 2 Comparisons of energy flow results in Case 1. Voltage results

Temperature results

Method

NGDD

GDD

Node

|V| (p.u.)

θ (°)

|V| (p.u.)

E1 E2 E3 E4 E5 E6 E7 E8 E9 E10 E11 E12 E13 E14

1.060 1.045 1.010 1.021 1.028 1.070 1.055 1.090 1.040 1.036 1.047 1.056 1.046 1.023

0 −0.984 −2.124 −3.897 −2.642 −0.344 −6.015 −6.015 −7.133 −6.218 −3.430 −1.267 −2.018 −5.922

1.060 1.045 1.010 1.021 1.028 1.070 1.055 1.090 1.040 1.036 1.047 1.056 1.046 1.023

Pressure results

Method

NGDD

GDD

θ (°)

Node

Ts (°C)

Tr (°C)

Ts (°C)

0 −0.984 −2.124 −3.897 −2.642 −0.344 −6.015 −6.015 −7.133 −6.218 −3.430 −1.267 −2.018 −5.922

H1 H2 H3 H4 H5 H6 H7 H8 H9 H10 H11 H12 H13 H14

99.651 99.375 99.358 98.959 98.744 98.564 98.475 96.387 98.909 99.493 99.218 99.905 100 100

49.565 50 49.612 49.648 49.903 50 50 50 49.689 49.564 50 49.644 49.412 49.602

99.651 99.375 99.358 98.959 98.744 98.564 98.475 96.387 98.909 99.493 99.218 99.905 100 100

Method

NGDD

GDD

Tr (°C)

Node

p (bar)

p (bar)

49.565 50 49.612 49.648 49.903 50 50 50 49.689 49.564 50 49.644 49.412 49.602

G1 G2 G3 G4 G5 G6 G7 G8 G9 G10

50 58.687 49.621 49.631 58.785 58.786 49.301 49.345 59.161 59.214

50 58.687 49.621 49.631 58.785 58.786 49.301 49.345 59.161 59.214

Fig. 4. Energy transmission results in Case 1.

return temperature Tr, and natural gas pressure p, and the energy flow results of the GDD method and the NGDD method are the same. Therefore, the GDD method and the NGDD method have the same computational accuracy. Fig. 4 shows the energy transmission results of electricity, heat and natural gas networks. In the electricity network, the voltage magnitudes are between 0.9 and 1.1 p.u., satisfying the operating requirements. The power flow is transported from the GG at Node E6, the CHPI unit at Node E1, and the CHPII unit at Node E2 to the load. In the heat network, the mass flow is transported from the CHPI unit at Node H13 and the CHPII unit at Node H14 to the load to provide heat power for the load. During the transmission, due to energy loss, the supply temperature decreases continuously from 100 °C at Nodes H13 and H14 to 96.387 °C at Node H8, and the supply temperature decreases with the increase of transmission distance. In the natural gas network, the

Table 3 Scenarios with different compressor operating modes in Case 1. Scenario

Compressor connection position Electricity node

Natural gas node

1

E4, E9

2

E4, E9

3

E4, E9

Inlet: G7, G8 Outlet: G9, G10 Inlet: G7, G8 Outlet: G9, G10 Inlet: G7, G8 Outlet: G9, G10

Operating mode

I, I I, II I, III

NGDD methods are compared. Table 2 shows the comparisons of energy flow results. As can be seen from Table 2, the energy flow results include voltage magnitude |V|, voltage angle θ, supply temperature Ts, 7

Energy Conversion and Management 207 (2020) 112534

J. Li, et al.

natural gas flow is transported from the natural gas station at Node G1 to the load, and the natural gas flow consumed by the GG at Node G3, the CHPI unit at Node G6 and the CHPII unit at Node G5 is also provided. As the natural gas flow from Node G1 to Nodes G7 and G8, the pressure gradually decreases. After compensated by the two compressors, the pressures of Nodes G9 and G10 are increased. These results are reasonable and verify the accuracy of the energy flow.

Table 4 Feasibility of energy flow models in different scenarios in Case 1. Scenario

Constant model

Diversified model

1 2 3

Feasible Infeasible Infeasible

Feasible Feasible Feasible

4.1.3. Analysis of compressor operating mode diversity For the convenience of expression, the energy flow model considering only the constant compression ratio of the compressor is called the constant model, and the energy flow model considering the diversity of compressor operating modes is called the diversified model. To analyze the influence of compressor operating modes on the energy flow of the natural gas network, the energy flow model feasibility and the natural gas pressure are compared in different compressor operating modes. There are three scenarios in Table 3. Modes I, II, and III are explained in Section 2.1.3. Table 4 shows the feasibility of the constant model and the diversified model in these three scenarios. The constant model is not feasible in Scenarios 2 and 3 because it cannot cope with different compressor operating modes. For example, in Scenario 2, the compressor connected to electricity network Node E4 operates in Mode I, while the compressor connected to electricity network Node E9 operates in Mode II. The constant model cannot simulate the energy flow in Scenario 2 because different compressor operating modes are considered simultaneously. On the contrary, the diversified model can effectively deal with the energy flow problem in different scenarios. Fig. 5 demonstrates the pressure results in Scenarios 1, 2, and 3. The pressures of Nodes 2, 5, 6, 9, and 10 differ significantly in different scenarios. For example, the pressures at Node 2 in Scenarios 1, 2 and 3 are 58.687, 59.476 and 57.430 bar, respectively. The result shows that the operating modes of compressors affect the energy flow results of the natural gas network. Concomitantly, as the compressors are connected to the electricity network, the power flow results of the electricity network are also affected, and the results of the electricity network also have an impact on the results of the heat network. Therefore, the operating modes of compressors will affect the energy flow results of the IES, and the diversified model can more accurately reflect the steadystate operation of the IES.

Fig. 5. Pressure results for different scenarios in Case 1.

4.1.4. Comparisons of computational speed In order to show the efficiency of the GDD method in solving the energy flow of the IEHGES, convergence curves and computational time of GDD and NGDD methods are compared. The convergence curves and the computational time are given in Fig. 6. As can be seen from Fig. 6, with the mismatch tolerance of 10-5, the number of iterations of the GDD method is only 7, while the NGDD method requires 19 iterations. Hence, the maximum mismatch of the GDD method decreases more rapidly. The reason why the convergence speed of the GDD method is faster is that the model is solved along the iterative direction of gradient descent and the solution of the model can be found faster by the GDD method. The computational time of NGDD

Fig. 6. Convergence curves and computational time with the mismatch tolerance of 10−5 in Case 1.

Table 5 Number of iterations and computational time with different mismatch tolerances in Case 1. Method

NGDD

Mismatch tolerance

Number of iterations

Computational time (s)

Number of iterations

Computational time (s)

13 16 19 22 25 28

0.027 0.029 0.031 0.033 0.035 0.038

6 7 7 8 9 9

0.018 0.020 0.020 0.021 0.023 0.023

−3

10 10−4 10−5 10−6 10−7 10−8

GDD

8

Percentage of time decrease

33.333% 31.034% 35.484% 36.364% 34.286% 39.474%

Energy Conversion and Management 207 (2020) 112534

J. Li, et al.

Table 6 Connection position and parameters of coupling components in Case 2. Components

Connection position

GG CHPI HP EB GB CP Compressor

Parameters

Electricity node

Heat node

Natural gas node

E339, E332, E155, E550 E1 E1030 E102 — E1, E1030, E102, E379 E29, E180, E414, E265, E945, E602

— H30 H32 H29 H31 H30, H32, H29, H31 —

G12, G16, G27, G14 G47 — — G42 — Inlet: G12, G20, G20, G21, G48, G24 Outlet: G13, G21, G48, G22, G25, G46

Hgas: 39 MJ/m3, ηGG: 0.55 zCHPI: 1.25, ηCHPI: 0.35 ηHP: 3 ηEB: 0.95 ηGB: 0.92 ga: 9.80 m/s2, ηCP: 0.65 α: 1.27, ηc: 0.8

Table 7 Comparisons of energy flow results in Case 2. Voltage results

Temperature results

Method

NGDD

GDD

Node

|V| (p.u.)

θ (°)

|V| (p.u.)

E1 E117 E233 E349 E465 E581 E697 E813 E929 E1045

1.060 0.987 1.000 1.040 1.011 1.041 0.963 1.003 1.010 1.033

1.640 −6.487 −9.589 1.184 33.955 5.513 −20.637 13.595 1.340 9.410

1.060 0.987 1.000 1.040 1.011 1.041 0.963 1.003 1.010 1.033

Pressure results

Method

NGDD

GDD

θ (°)

Node

Ts (°C)

Tr (°C)

Ts (°C)

1.640 −6.487 −9.589 1.184 33.955 5.513 −20.637 13.595 1.340 9.410

H1 H4 H7 H10 H13 H16 H19 H22 H25 H28

69.850 69.800 68.559 69.620 69.542 68.765 68.730 68.842 69.030 69.451

29.716 29.647 30 29.742 29.762 30 30 30 30 30

69.850 69.800 68.559 69.620 69.542 68.765 68.730 68.842 69.030 69.451

Method

NGDD

GDD

Tr (°C)

Node

p (bar)

p (bar)

29.716 29.647 30 29.742 29.762 30 30 30 30 30

G1 G6 G11 G16 G21 G26 G31 G36 G41 G46

50 49.903 48.236 57.085 58.695 63.552 63.292 63.272 63.269 63.405

50 49.903 48.236 57.085 58.695 63.552 63.292 63.272 63.269 63.405

Table 8 Scenarios with different compressor operating modes in Case 2. Scenario

Compressor connection position

Operating mode

Electricity node

Natural gas node

1

E29, E180, E414, E265, E945, E602

2

E29, E180, E414, E265, E945, E602

3

E29, E180, E414, E265, E945, E602

Inlet: G12, G20, G20, G21, G48, G24 Outlet: G13, G21, G48, G22, G25, G46 Inlet: G12, G20, G20, G21, G48, G24 Outlet: G13, G21, G48, G22, G25, G46 Inlet: G12, G20, G20, G21, G48, G24 Outlet: G13, G21, G48, G22, G25, G46

Constant model

Diversified model

1 2 3

Feasible Infeasible Infeasible

Feasible Feasible Feasible

I, II, II, I, I, I I, II, II, I, I, III

The percentages of time decrease are more than 30%, and the average percentage of time decrease is 34.996%.

Table 9 Feasibility of energy flow models in different scenarios in Case 2. Scenario

I, I, I, I, I, I

4.2. Case 2: the 1047 node electricity-32 node heat-48 node gas integrated system In this section, the system structure and simulation parameters of the 1047 node electricity-32 node heat-48 node gas integrated system are introduced, and the results of computational accuracy, compressor operating mode diversity and computational speed are analyzed.

and GDD methods is 0.031 and 0.020 s, respectively. The computational time of the GDD method is reduced by 35.484% than that of the NGDD method. Hence, the computational speed is enhanced greatly by the GDD method. In order to further analyze the convergence characteristics and the computational speed of NGDD and GDD methods, the number of iterations and computational time with different mismatch tolerances are compared. Table 5 shows the number of iterations and computational time with different mismatch tolerances. When the mismatch tolerance changes from 10−3 to 10−8, the number of iterations of the NGDD method changes from 13 to 28, while that of the GDD method changes from 6 to 9. Therefore, the GDD method has faster convergence speed. Moreover, with different mismatch tolerances, the computational time of the GDD method is less than that of the NGDD method.

4.2.1. Description of Case 2 To further analyze the influence of compressor operating modes on energy flow results and verify the superiority of the GDD method in solving the energy flow problem, the large-scale IES in Case 2 is simulated. The system consists of a 1047-node IEEE test system, a 32-node heat network [19], and a 48-node natural gas network [28]. There are 4 heat sources and 32 pipes in the heat network. In addition, 8 sources and 45 pipelines are contained in the natural gas network. It can be seen from the connection position and parameters of the coupling components shown in Table 6 that there are 4 GGs, 1 CHPI unit, 1 HP, 1 EB, 1 GB, 4 9

Energy Conversion and Management 207 (2020) 112534

J. Li, et al.

preset reference voltage magnitude and angle are 1.05 p.u. and 0 degree respectively, the voltage magnitude and angle change around the preset value and meet the operating requirements. In the heat system, the supply temperature of heat sources is 70 °C, so the supply temperatures of load Nodes H1–H28 are less than 70 °C. In addition, the return temperatures of the non-mixing Nodes H7, H16, H19, H22, H25, and H28 are equal to the outlet temperatures of the nodes. In the natural gas system, when the pressure of gas sources is 50 bar, the pressures of Nodes G6 and G11 are less than 50 bar because of pressure loss, and the pressures of Nodes G16, G21, G26, G31, G36, G41, and G46 are more than 50 bar due to the supplement of compressors. From the energy flow results presented in Table 7, it is clear that the results of the GDD method are the same as those of the NGDD method; thus, the computational accuracy of the GDD method is the same as that of the NGDD method. 4.2.3. Analysis of compressor operating mode diversity The preset scenarios are presented in Table 8. There are three scenarios. Modes I, II, and III are explained in Section 2.1.3. In Scenario 2, the compressors connected to electricity network Nodes E29, E265, E945, and E602 operate in Mode I, and the compressors connected to electricity network Nodes E180 and E414 operate in Mode II. Table 9 provides the feasibility of the constant model and the diversified model in these three scenarios. The constant model cannot simulate the energy flow in Scenarios 2 and 3 with different compressor operating modes, but the diversified model can effectively deal with the energy flow problem in all scenarios. Fig. 7 shows the pressure results in Scenarios 1, 2, and 3. There are significant differences in the pressure results of Nodes 21, 26, 31, 36, 41, and 46. For example, the pressures at Node 46 in Scenarios 1, 2 and 3 are 63.405, 61.351 and 62.896 bar, respectively. In the meantime, the energy flow results of the electricity network and the heat network also change in different scenarios. Consequently, it is necessary to consider the diversity of compressor operating modes.

Fig. 7. Pressure results for different scenarios in Case 2.

4.2.4. Comparisons of computational speed To demonstrate the superiority of the GDD method in solving the large-scale energy flow problem, the convergence curves and the computational time of the IES in Case 2 are analyzed. Fig. 8 shows the convergence curves and the computational time of the GDD method and the NGDD method. With the mismatch tolerance of 10−5, the number of iterations of the GDD method is only 6, whereas that of the NGDD method is 20. Therefore, the maximum mismatch of the GDD method decreases more rapidly. Moreover, the computational time of NGDD and GDD methods is 0.356 and 0.123 s, respectively. The computational time of the GDD method is reduced by 65.449% than that of the NGDD method. Therefore, the GDD method is more efficient in solving the energy flow problem of the large-scale IES. The number of iterations and the computational time with different mismatch tolerances are provided in Table 10. When the mismatch tolerance changes from 10−3 to 10−8, the number of iterations of the NGDD method increased rapidly from 14 to 29, while that of the GDD method only increased from 6 to 8. In addition, with different mismatch

Fig. 8. Convergence curves and computational time with the mismatch tolerance of 10−5 in Case 2.

CPs, and 6 compressors. The variables to be solved include 1047 voltage magnitudes |V|, 1047 voltage angles θ, mass flow m of 32 pipes, 32 supply temperatures Ts, 32 return temperatures Tr, 48 pressure variables Π, and natural gas flow fu of 14 non-pipeline elements.

4.2.2. Comparisons of computational accuracy The energy flow results including voltage magnitude |V|, voltage angle θ, supply temperature Ts, return temperature Tr, and natural gas pressure p are shown in Table 7. In the electric power system, since the

Table 10 Number of iterations and computational time with different mismatch tolerances in Case 2. Method

NGDD

Mismatch tolerance

Number of iterations

Computational time (s)

Number of iterations

Computational time (s)

14 17 20 23 26 29

0.272 0.313 0.356 0.404 0.447 0.486

6 6 6 7 8 8

0.123 0.123 0.123 0.147 0.162 0.162

−3

10 10−4 10−5 10−6 10−7 10−8

GDD

10

Percentage of time decrease

54.779% 60.703% 65.449% 63.614% 63.758% 66.667%

Energy Conversion and Management 207 (2020) 112534

J. Li, et al.

tolerances, the computational time of the GDD method is less than that of the NGDD method. The percentages of time decrease are more than 50%, and the computational time of the GDD method is reduced by 62.495% than that of the NGDD method on average.

[2]

[3]

5. Conclusion [4]

In this study, an energy flow model for the IEHGES that considers the diversity of compressor operating modes is established. Further, by calculating the partial derivatives of each state variable and deriving the Jacobian matrix, a gradient descent iterative method for solving the energy flow model is developed. The case studies presented illustrate the necessity of considering the diversity of compressor operating modes. The numerical results show that compared with the NGDD method, the number of iterations and the computational time of the GDD method are significantly reduced under the same mismatch tolerance. The results also verify the efficiency of the GDD method in solving the large-scale energy flow problem of the IEHGES. Using the model and method presented in this paper, the energy flow of the IEHGES containing compressors with different operating modes can be solved conveniently and effectively to provide more accurate and efficient energy flow results. In the IES, different energy systems have different dynamic processes. In the power system, the transmission speed of electrical energy is very fast and the inertia is very small. On the other hand, the transmission speed of heat energy and natural gas is slow and the inertia is large. Therefore, more research needs to be done in the future. Future work includes: the dynamic modelling of heat networks, natural gas networks, and coupling components; the dynamic coupling among electricity, heat, and natural gas networks; the energy flow considering the dynamic characteristics; and the integrated optimal energy flow.

[5] [6]

[7] [8] [9]

[10]

[11]

[12]

[13]

[14]

[15] [16]

CRediT authorship contribution statement

[17]

Jinghua Li: Conceptualization, Writing - review & editing, Supervision, Project administration, Funding acquisition. Yujin Huang: Methodology, Software, Validation, Formal analysis, Writing - original draft. Mengshu Zhu: Investigation, Resources, Data curation, Visualization.

[18] [19] [20]

[21]

Declaration of Competing Interest

[22]

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

[23]

Acknowledgment

[24]

This work was supported in part by the National Natural Science Foundation of China under Grant 51977042 and 51377027.

[25]

[26]

Appendix A. Supplementary data [27]

Supplementary data to this article can be found online at https:// doi.org/10.1016/j.enconman.2020.112534.

[28]

References [1] Zhou J, Wu Y, Wu C, Deng Z, Xu C, Hu Y. A hybrid fuzzy multi-criteria decision-

11

making approach for performance analysis and evaluation of park-level integrated energy system. Energy Convers Manage 2019;201:112134. Huang W, Zhang N, Kang C, Li M, Huo M. From demand response to integrated demand response: review and prospect of research and application. Protection Control Modern Power Syst 2019;4(4):148–50. Koirala BP, Koliou E, Friege J, Hakvoort RA, Herder PM. Energetic communities for community energy: a review of key issues and trends shaping integrated community energy systems. Renew Sustain Energy Rev 2016;56:722–44. Massrur HR, Niknam T, Fotuhi-Firuzabad M. Investigation of carrier demand response uncertainty on energy flow of renewable-based integrated electricity–gas–heat systems. IEEE Trans Ind Inf 2018;14(11):5133–42. Huang Z, Yu H, Chu X, Peng Z. Energetic and exergetic analysis of integrated energy system based on parametric method. Energy Convers Manage 2017;150:588–98. Li J, Fang J, Zeng Q, Chen Z. Optimal operation of the integrated electrical and heating systems to accommodate the intermittent renewable sources. Appl Energy 2016;167:244–54. Fu X, Sun H, Guo Q, Pan Z, Xiong W, Wang L. Uncertainty analysis of an integrated energy system based on information theory. Energy 2017;122:649–62. Zhou Y, Gu C, Wu H, Song Y. An equivalent model of gas networks for dynamic analysis of gas-electricity systems. IEEE Trans Power Syst 2017;32(6):4255–64. Liu X, Mancarella P. Modelling, assessment and Sankey diagrams of integrated electricity-heat-gas networks in multi-vector district energy systems. Appl Energy 2016;167:336–52. Sun Q, Chen Y. Multi-energy flow calculation method for we-energy based energy internet. 1st IEEE International Conference on Energy Internet (ICEI 2017). Beijing, China. 2017. p. 30–5. Martinez-Mares A, Fuerte-Esquivel CR. A unified gas and power flow analysis in natural gas and electricity coupled networks. IEEE Trans Power Syst 2012;27(4):2156–66. Zeng Q, Fang J, Li J, Chen Z. Steady-state analysis of the integrated natural gas and electric power system with bi-directional energy conversion. Appl Energy 2016;184:1483–92. Shabanpour-Haghighi A, Seifi AR. An integrated steady-state operation assessment of electrical, natural gas, and district heating networks. IEEE Trans Power Syst 2016;31(5):3636–47. Massrur HR, Niknam T, Aghaei J, Shafie-khah M, Catalão JPS. Fast decomposed energy flow in large-scale integrated electricity-gas-heat energy systems. IEEE Trans Sustain Energy 2018;9(4):1565–77. Vidović D, Sutlović E, Majstrović M. Steady state analysis and modeling of the gas compressor station using the electrical analogy. Energy 2019;166:307–17. Osiadacz A. Simulation and analysis of gas networks. USA: Gulf Publish Company; 1987. Meng Q, Guan Q, Jia N, Zhang L, Wang Y. An improved sequential energy flow analysis method based on multiple balance nodes in gas-electricity interconnection systems. IEEE Access 2019;7:95487–95. Liu X, Wu J, Jenkins N, Bagdanavicius A. Combined analysis of electricity and heat networks. Appl Energy 2016;162:1238–50. Liu X. Combined analysis of electricity and heat networks. Cardiff: Cardiff University; 2013. Yang Z, Gao C, Zhao M. Coordination of integrated natural gas and electrical systems in day-ahead scheduling considering a novel flexible energy-use mechanism. Energy Convers Manage 2019;196:117–26. Li H, Kang S, Lu L, Liu L, Zhang X, Zhang G. Optimal design and analysis of a new CHP-HP integrated system. Energy Convers Manage 2017;146:217–27. Anand H, Narang N, Dhillon JS. Unit commitment considering dual-mode combined heat and power generating units using integrated optimization technique. Energy Convers Manage 2018;171:984–1001. Cao Y, Wei W, Wu L, Mei S, Shahidehpour M, Li Z. Decentralized operation of interdependent power distribution network and district heating network: a marketdriven approach. IEEE Trans Smart Grid 2019;10(5):5374–85. Wang H, Zhang R, Peng J, Wang G, Liu Y, Jiang H. GPNBI-inspired MOSFA for Pareto operation optimization of integrated energy system. Energy Convers Manage 2017;151:524–37. Ma T, Wu J, Hao L. Energy flow modeling and optimal operation analysis of the micro energy grid based on energy hub. Energy Convers Manage 2017;133:292–306. Shabanpour-Haghighi A, Seifi AR. Simultaneous integrated optimal energy flow of electricity, gas, and heat. Energy Convers Manage 2015;101:579–91. Qiao Z, Guo Q, Sun H, Pan Z, Liu Y, Xiong W. An interval gas flow analysis in natural gas and electricity coupled networks considering the uncertainty of wind power. Appl Energy 2017;201:343–53. Chen S, Wei Z, Sun G, Cheung KW, Sun Y. Multi-linear probabilistic energy flow analysis of integrated electrical and natural-gas systems. IEEE Trans Power Syst 2017;32(3):1970–9.