Optimal operation of integrated energy system including power thermal and gas subsystems

Abstract

As a form of hybrid multi-energy systems, the integrated energy system contains different forms of energy such as power, thermal, and gas which meet the load of various energy forms. Focusing mainly on model building and optimal operation of the integrated energy system, in this paper, the dist-flow method is applied to quickly calculate the power flow and the gas system model is built by the analogy of the power system model. In addition, the piecewise linearization method is applied to solve the quadratic Weymouth gas flow equation, and the alternating direction method of multipliers (ADMM) method is applied to narrow the optimal results of each subsystem at the coupling point. The entire system reaches its optimal operation through multiple iterations. The power-thermal-gas integrated energy system used in the case study includes an IEEE-33 bus power system, a Belgian 20 node natural gas system, and a six node thermal system. Simulation-based calculations and comparison of the results under different scenarios prove that the power-thermal-gas integrated energy system enhances the flexibility and stability of the system as well as reducing system operating costs to some extent.

Graphical abstract

Keywords

integrated energy system / power-to-gas / dist-flow / piecewise linearization / alternating direction method of multipliers (ADMM)

1 Introduction

In recent years, with an increasing number of high-pressure natural gas pipelines and gas-fired units installed in the power system, power and natural gas systems have become highly coupled with each other. In addition, with the improvement of the efficiency of the power to gas (PtG) technology and the broad application of the combined heat and power (CHP) unit, the three different forms of energy, i.e., electricity, heat, and gas are closely coupled together. In the power thermal and gas integrated energy system, various forms of energy are monitored and transformed into each other at the coupling point. Load forecasting directly affects the economy of system operation [1], and different forms of energy conversion in the integrated energy system can better meet customer demand for different types of energy demands [2].

Although the power system and the natural gas system have different network structures and physical properties, the coupling and economy of the two system infrastructures are interdependent, and the two systems can still be analyzed separately [3]. There are many methods to calculate the power flow of an independent power system. Considering the natural gas system model can be formed by the analogy to the power system, the dist-flow method is selected to calculate the power flow in this article. This method can quickly calculate the voltage amplitude of each bus, the active and reactive power values of each branch to meet the optimal operating conditions in the integrated energy system.

In the power-thermal-gas integrated energy system, it not only includes multiple energy subsystems but also considers the coordinated operation of each subsystem after being coupled by an energy conversion device [4]. Thus, the production and operation cost of the entire system is minimized through the coordinated dispatch and unified management of the integrated energy system that meets the load of each subsystem [5]. After that, renewable energy can be added to improve the utilization rate of renewable energy and reduce carbon emissions [6]. To achieve the optimal operation condition, it is important to modify the security constraints and coordinated scheduling of coupled energy systems to reduce the potential safety hazards of the energy system, prevent operational failures, and save operation costs [7]. The optimal operating condition can be achieved by optimizing the entire system with coupling constraints within each system [8], and the operation costs of the system that satisfies the load of each subsystem can be reduced.

Although the thermal system has its network structure and load demand, it plays a role in coupling the power and gas subsystems in the entire optimization calculation [9]. However, the power subsystem and the gas subsystem are independent in a certain degree and the optimal operating conditions of their respective subsystems can be calculated separately [10]. To achieve the optimal operating conditions of the entire integrated energy system, it is necessary to use the alternating direction method of multipliers(ADMM) algorithm to calculate the optimal operating conditions of the power and gas subsystems at the coupling points [11]. The only information that needs to be exchanged when the ADMM algorithm used is the different optimal values at the coupling point energy conversion device, including gas-fired units and CHP units [12]. At the same time, these two subsystems can still keep their privacy and adopt corresponding models and algorithms to solve their system optimization calculations.

This paper is contributive because, for the overall optimal operation of the integrated energy system, the dist-flow method is used to quickly solve the bus voltage and the active power of each branch in the power subsystem. The natural gas model is established by analogy with the power system model, and the piecewise linearization method is used to address the quadratic gas flow function, which facilitates the calculation of the optimal operation of each subsystem. In addition, the power subsystem and the gas subsystem calculate the optimal situation of the subsystem according to their conditions separately. The only exchanged information is the calculation results of the two subsystems at the coupling point. Moreover, the ADMM algorithm and multiple iteration methods are used to narrow the deviation between the different optimal results to the tolerance. This method provides a basis for future research on the privacy protection settings of the subsystems and multi-agent collaborative control in the integrated energy system.

2 Power, gas, and thermal subsystem model

2.1 Power subsystem model

In the power-thermal-gas integrated energy system, the regular AC power flow model is adopted in the power system, whose model can be expressed by Eq. (1),

P=\mathrm{Re}\left\{U{\left(YU\right)}^{\mathrm{*}}\right\},

(1)Q=\mathrm{Im}\left\{U{\left(YU\right)}^{\mathrm{*}}\right\},

where P, Q, Y, and U represent active power, reactive power, system admittance matrix, and bus voltage vector, respectively.

2.1.1 Dist-Flow method

The Dist-Flow method is proposed by Felix F. Wu and Mesut E. Baran to calculate the power flow of the power system that especially focused on the distribution network [13]. The Dist-Flow model has several variables in a single radial network structure, including node voltage, active and reactive power of each branch, as well as active and reactive power of load demand (e.g., 33-bus power network). The structure of a single radial power system network is shown in Fig. 1.

(2a)\{\begin{array}{l}{P}_{i+\mathrm{1}}=P_i+P_{i+\mathrm{1}}^{\mathrm{G}}-P_{i+\mathrm{1}}^{\mathrm{D}}-P_{i+\mathrm{1}}^l-P_i^{\mathrm{loss}},\hfill \\ Q_{i+\mathrm{1}}=Q_i+Q_{i+\mathrm{1}}^{\mathrm{G}}-Q_{i+\mathrm{1}}^{\mathrm{D}}-Q_{i+\mathrm{1}}^l-Q_i^{\mathrm{loss}},\hfill \\ V_{i+\mathrm{1}}^{\mathrm{2}}=V_i^{\mathrm{2}}-\mathrm{2}\left(r_i{P}_i+x_i{Q}_i\right)+\left(r_i^{\mathrm{2}}+x_i^{\mathrm{2}}\right){\displaystyle \frac{P_i^{\mathrm{2}}+Q_i^{\mathrm{2}}}{V_i^{\mathrm{2}}}},\hfill \end{array}

(2b)P_i^{\mathrm{loss}}=r_i{\displaystyle \frac{P_i^{\mathrm{2}}+Q_i^{\mathrm{2}}}{V_i^{\mathrm{2}}}},

(2c)Q_i^{\mathrm{loss}}=x_i{\displaystyle \frac{P_i^{\mathrm{2}}+Q_i^{\mathrm{2}}}{V_i^{\mathrm{2}}}},

(2d)S_{i+\mathrm{1}}^{\mathrm{E}}=S_i^{\mathrm{E}}+S_i^{\mathrm{E},\mathrm{DG}}-S_i^{\mathrm{E},\mathrm{D}}-S_i^{E,l^{\prime }}-S_i^{\mathrm{E},\mathrm{loss}},

where S_i^{\mathrm{E}}=P_i+j{Q}_i is the power flow of the head of branch i while S_i^{\mathrm{E},\mathrm{DG}} is the generated power from distributed generator; S_i^{\mathrm{E},\mathrm{D}}=P_i^{\mathrm{D}}+j{Q}_i^{\mathrm{D}} is the load demand at branch, i while \mathrm{}{S}_i^{E,l^{\prime }} is the lateral branch power flow of branch i; and P_i^{\mathrm{loss}} and Q_i^{\mathrm{loss}} are lost active and reactive power at branch i respectively.

The Dist-Flow method is broadly used in the distribution system network because the magnitude of node voltage is the point needing to be considered as the angle of nodal voltage could be ignored. In addition, it is suitable for the single radial network structure. Moreover, the active and reactive power of the branch have been introduced.

2.1.2 Power subsystem constraints

The constraints of the power system contain load balance constraints (Eq. (3)), branch power flow limits (Eq. (4)), bus voltage magnitude limits (Eq. (5)), gas-fired generator, and PtG system capacity (Eqs. (6)–(8)). When there are multiple traditional generators in the power subsystem, the minimum ON/OFF time and ramping UP/DOWN limits (Eqs. (9)–(12)), startup and shutdown costs (Eqs. (13)–(14), system spinning reserve demand (Eq. (15)) and system stability requirements (Eq. (16)) should also be considered.

{\sum }_{i\in N_{\mathrm{E}}}{P}^{\mathrm{E}}+{\sum }_{i\in N_{\mathrm{DG}}}{P}^{\mathrm{DG}}+{\sum }_{i\in N_{\mathrm{GE}}}{P}^{\mathrm{GE}}

(3)-{\sum }_{i\in N_{\mathrm{PtG}}}{P}^{\mathrm{PtG}}-{\sum }_{i\in N_{\mathrm{EB}}}{P}^{\mathrm{EB}}={\sum }_{i\in N_{\mathrm{E}}}{P}^{\mathrm{load}},

(4)S_i\le S_i^{\mathrm{Max}},

(5)V_i^{\min }\le V_i\le V_i^{\mathrm{Max}},

(6)\{\begin{array}{l}{P}_i^{\mathrm{GE},\mathrm{Min}}\le P_i^{\mathrm{G}}\le P_i^{\mathrm{GE},\mathrm{Max}},\hfill \\ Q_i^{\mathrm{GE},\mathrm{Min}}\le Q_i^{\mathrm{G}}\le Q_i^{\mathrm{GE},\mathrm{Max}},\hfill \end{array}

(7)P_i^{\mathrm{GE},\mathrm{Min}}\le P_i\le {\lambda }^{\mathrm{ge}}{P}_i^{\mathrm{Gas}},(\mathrm{if}\mathrm{\hspace{0.17em}}\mathrm{\hspace{0.17em}}i\in N_{\mathrm{GE}})

(8)\mathrm{0}\le P^{\mathrm{PtG}}\le P^{\mathrm{PtG},\max },

(9)(X_{i\left(t-\mathrm{1}\right)}^{\mathrm{on}}-T_{it}^{\mathrm{on}})\left(I_{i\left(t-\mathrm{1}\right)}^{\mathrm{b}}-I_{it}^{\mathrm{b}}\right)\ge \mathrm{0},

(10)(X_{i\left(t-\mathrm{1}\right)}^{\mathrm{off}}-T_{it}^{\mathrm{off}})\left(I_{it}^{\mathrm{b}}-I_{i\left(t-\mathrm{1}\right)}^{\mathrm{b}}\right)\ge \mathrm{0},

P_{it}^{\mathrm{b}}-P_{i\left(t-\mathrm{1}\right)}^{\mathrm{b}}\le U{R}_i{I}_{i\left(t-\mathrm{1}\right)}^{\mathrm{b}}+P_i^{\min }\left(I_{it}^{\mathrm{b}}-I_{i\left(t-\mathrm{1}\right)}^{\mathrm{b}}\right)+P_i^{\max }\left(\mathrm{1}-I_{it}^{\mathrm{b}}\right),

P_{it}^{\mathrm{b}}-P_{i\left(t-\mathrm{1}\right)}^{\mathrm{b}}\le U{R}_i{I}_{i\left(t-\mathrm{1}\right)}^{\mathrm{b}}+P_i^{\min }\left(I_{it}^{\mathrm{b}}-I_{i\left(t-\mathrm{1}\right)}^{\mathrm{b}}\right)

(11)+P_i^{\max }\left(\mathrm{1}-I_{it}^{\mathrm{b}}\right),

P_{i\left(t-\mathrm{1}\right)}^{\mathrm{b}}-P_{it}^{\mathrm{b}}\le U{R}_i{I}_{it}^{\mathrm{b}}+P_i^{\min }\left(I_{i\left(t-\mathrm{1}\right)}^{\mathrm{b}}-I_{it}^{\mathrm{b}}\right)

(12)+P_i^{\max }\left(\mathrm{1}-I_{i\left(t-\mathrm{1}\right)}^{\mathrm{b}}\right),

(13)S{U}_{it}^{\mathrm{b}}\ge s{u}_i\left(I_{it}^{\mathrm{b}}-I_{i\left(t-\mathrm{1}\right)}^{\mathrm{b}}\right),\mathrm{\hspace{0.17em}}S{U}_{it}^{\mathrm{b}}\ge \mathrm{0},

(14)S{D}_{it}^{\mathrm{b}}\ge s{d}_i\left(I_{i\left(t-\mathrm{1}\right)}^{\mathrm{b}}-I_{it}^{\mathrm{b}}\right),\mathrm{\hspace{0.17em}}S{D}_{it}^{\mathrm{b}}\ge \mathrm{0},

(15)\sum _i{P}_i^{\max }{I}_{it}^{\mathrm{b}}\ge \sum _{\mathrm{d}}{P}_{\mathrm{d}t}^{\mathrm{b}}+S{R}_t,

(16)I_{it}^{\mathrm{b}}+I_t^{\mathrm{b},\mathrm{PtG}}\le \mathrm{1},\mathrm{if}\mathrm{\hspace{0.17em}}\mathrm{\hspace{0.17em}}{K}_{\mathrm{p}}\left(e,i\right)=K_{\mathrm{a}},\mathrm{}\mathrm{\hspace{0.17em}}\left(e,a\right)=\mathrm{1.}

2.2 Gas subsystem model

The simple natural gas system structure is depicted in Fig. 2, which normally contains a gas well, gas pipelines, a compressor, and gas loads [14]. The natural gas system model is established analogous to Dist-Flow in power systems, which includes the GasFlow model and the compressor model. The variables in this model include the node pressure, gas flow through each pipeline, compressor, and gas load demand.

Nodes and branches are two key points of the gas system. The node can be divided into the known-injection node and the known-pressure node, based on nodal pressure and gas injection [15]. Generally, the load nodes, source loads, and junctions with no gas injunctions belong to known-injection node, whose gas injection {\omega }_iis already known and pressure p_ineeds to be determined. Besides, the pressure of known-pressure node serves as a reference for the other nodal pressure, which is typically one of the source nodes.

For gas branches, the pipelines and compressors should be taken into account. The only variable of the pipeline flow model is the flow S_{ij}^{\mathrm{G}} which depends on the adjacent nodal pressure p_iand p_j. The pressure ratio between the outlet and the inlet gas pressure and the flow rate through the compressor are given in Section 2.2.2.

2.2.1 Steady-state gas flow model

In this paper, the Weymouth’s formula is applied to express the isothermal gas flow in a steady-state flow rate, where the gas flow S_{ij}^{\mathrm{G}} and direction signal sg{n}_{ij} are dependent on the node pressure. The general steady-state gas flow equation in transmission pipeline is given as

(17a)S_{ij}^{\mathrm{G}}=sg{n}_{ij}{K}_{\mathrm{G}}{\displaystyle \frac{T_{\mathrm{0}}}{p_{\mathrm{0}}}}\sqrt{sg{n}_{ij}{\displaystyle \frac{\left(p_i^{\mathrm{2}}-p_j^{\mathrm{2}}\right)D_{ij}^{\mathrm{5}}}{F_{ij}G{L}_{ij}{T}_{\mathrm{a}}Z}}},

(17b)K_{\mathrm{G}}=\sqrt{{\displaystyle \frac{{\mathrm{\pi }}^{\mathrm{2}}{R}^{\mathrm{air}}}{\mathrm{64}}}}=\mathrm{3.2387},

(17{{\rm\tf="fp1"\char"0063}})

where ρ_i and ρ_j are the nodal pressure at nodes i and j respectively; p_{\mathrm{0}} is the normal condition gas pressure; L_ij and D_ij are the length and internal diameter of pipeline i to j respectively; G is the specific gravity ratio as R^air/R^gas (usually G = 1/0.6), and K_G is the air constant in relation to P^{\mathrm{air}}; T_a and T₀ are the average temperature of gas and normal condition temperature of environment respectively; and Z and F_ij are the gas compressibility factor and dimensionless friction factor, respectively. The gas load can be expressed by the gas volumetric flow rate (m³/s) or the energy flow rate (TJ/h). 1 TJ/h equal to 0.278 GJ/s, which is approximately equal to 278 MW.

2.2.2 Compressor model

In the process of gas transmission, the gas flows lose the initial pressure due to the friction resistance in the gas pipeline and consumption by users. Therefore, compressor stations are installed to compensate the energy in the network to solve the pressure loss problem and transport the gas. A gas turbine is used to drive the compressor, which is operated by the gas extracted from the pipeline. Horsepower consumption is the energy consumed by a compressor, which depends on the pressure ratio between the outlet and the inlet gas pressure and the amount of gas flows through the compressor. The empirical equations are expressed as [16]

(18)S_{ij}^{\mathrm{G},\mathrm{cp}}=sgn\left(p_i,p_j\right){\displaystyle \frac{H{P}_{ij}}{{\omega }_{\mathrm{1},ij}-{\omega }_{\mathrm{2},ij}{\left[{\displaystyle \frac{\max \left(p_i,p_j\right)}{\min \left(p_i,p_j\right)}}\right]}^{{\omega }_{\mathrm{3},ij}}}},

where S_{ij}^{\mathrm{G},\mathrm{cp}}and HP_ij, are the gas flow through the compressor and the horsepower of the gas compressor, while ω_1,ij, ω_2,ij, and ω_3,ij are the coefficients of the gas compressor between node i and j.

2.2.3 Gas subsystem constraints

The gas system model has several parts. The constraints of each part include the node, the pipeline and the compressor, all of which need be considered [17]. In the gas system nodal and pipeline constraints part, the nodal gas flow balance constraints (Eq. (19)), the nodal pressure constraints (Eq. (20)), and the gas flow constraints (Eq. (21)) need to be satisfied. For the compressor part, the compressor pressure ratio r_i^{\mathrm{cp}} is bounded as Eq. (22). The amount of gas consumed by the compressor can be simply expressed by a fixed percentage of the gas flow as Eq. (23) and the gas flow through compressor constraints (Eq. (24)) need to be satisfied. The gas-fired unit equation is formulated as Eq. (25), and the ranges of gas supply are restricted as expressed in Eq. (26).

S_i^{\mathrm{G},\mathrm{S}}-S_i^{\mathrm{G},\mathrm{D}}-S_i^{\mathrm{G},\mathrm{F}}-{\sum }_{i\in N_{\mathrm{cp}}}{\tau }_i^{\mathrm{cp}}

(19)={\sum }_{i\in N_{\mathrm{G}}}{S}_{ij}^{\mathrm{G}}+{\sum }_{k\in N_{\mathrm{cp}}}{S}_{ij}^{\mathrm{G},\mathrm{cp}},

(20)p_i^{\min }\le p_i\le p_i^{\max },

(21)-S_{ij}^{\mathrm{G},\min }\le S_{ij}^{\mathrm{G}}\le S_{ij}^{\mathrm{G},\max },

(22)r_{i,\min }^{\mathrm{cp}}\le r_i^{\mathrm{cp}}={\displaystyle \frac{p_{i,\mathrm{out}}^{\mathrm{cp}}}{p_{i,\mathrm{in}}^{\mathrm{cp}}}}\le r_{i,\max }^{\mathrm{cp}},

(23){\tau }_i^{\mathrm{cp}}=C_i^{\mathrm{cp}}{S}_i^{\mathrm{G},\mathrm{cp}},

(24)\mathrm{0}\le S_i^{\mathrm{G},\mathrm{cp}}\le S_{i,\max }^{\mathrm{G},\mathrm{cp}},

(25)S_i^{\mathrm{G},\mathrm{F}}=C_{i-\mathrm{1}}^{\mathrm{G},\mathrm{F}}{P}_i^{\mathrm{G},\mathrm{F}}/Q,

(26)S_{i,\min }^{\mathrm{G},\mathrm{S}}\le S_i^{\mathrm{G},\mathrm{S}}\le S_{i,\max }^{\mathrm{G},\mathrm{S}},

where S_i^{\mathrm{G},\mathrm{S}}, S_i^{\mathrm{G},\mathrm{D}}, and S_i^{\mathrm{G},\mathrm{F}} are the gas source supply, the gas demand, and the gas-fired unit demand at gas branch i, respectively.

2.3 Thermal subsystem model

The thermal system normally consists of heat sources, heat medium, and heat loads [18]. The heat sources in this paper include the CHP unit and the electric boiler.

2.3.1 CHP unit and electric boiler

A typical CHP unit generates electricity by burning natural gas and captures heat from the turbine through a heat recovery device, which greatly enhances the coupling effect among power, thermal, and gas.

The feasible region of the CHP unit can be represented by the polygonal area in Fig. 3 [19]. The CHP unit can generate electricity and heat at the same time and the operating performance can be represented by the relationship between its electric power p_t^{\mathrm{CHP}}and thermal power {\varphi }_t^{\mathrm{CHP}}, which can be formulated by the extreme points in the feasible region of the CHP system at any time.

p_t^{\mathrm{CHP}}={{\displaystyle \sum }}_{k=\mathrm{1}}^K{\alpha }_t^k{P}^k,

(27){\varphi }_t^{\mathrm{CHP}}={{\displaystyle \sum }}_{k=\mathrm{1}}^K{\alpha }_t^k{\Phi }^k,

(28){\varphi }_t^{\mathrm{CHP}}=\eta p_t^{\mathrm{CHP}},

where P^k and {\Phi }^k are the electric power and thermal power at the extreme point k, \eta is an empirical value, (usually \eta = 2.58); K and α are the number of extreme operating points of the CHP unit at time t, which satisfy the following formula. The natural gas usage cost of the CHP unit is expressed as

(29){{\displaystyle \sum }}_{k=\mathrm{1}}^K{\alpha }_t^k=\mathrm{1},\mathrm{\hspace{0.17em}}\mathrm{\hspace{0.17em}}\mathrm{\hspace{0.17em}}\mathrm{\hspace{0.17em}}\mathrm{0}\le {\alpha }_t^k\le \mathrm{1},

(30)C_{\mathrm{g}}={\displaystyle \frac{C_{\mathrm{n}}g}{Q}}{\sum }^{\mathrm{}}{\displaystyle \frac{p_t^{\mathrm{CHP}}}{{\eta }^{\prime }}}\mathrm{\Delta }T,

where C_{\mathrm{n}}g andQ are the unit of natural gas price (\$/{\mathrm{m}}^{\mathrm{3}}) and the value of gross heating (usuallyQ = 9.7 kWh/m³).

The electric boiler is a normal heat source, which converts electricity into heat energy. The produced heat energy Q_i^{\mathrm{EB}}and the consumed electricity P_i^{\mathrm{EB}}satisfy

(31)Q_i^{\mathrm{EB}}={\eta }^{\prime }{P}_i^{\mathrm{EB}}.

The maximum output of the electric boiler is 3 MW and the value of {\eta }^{\prime } is 0.85 in this paper.

2.3.2 Thermal flow model

Water is normally used as the heat medium in civil heating because of its large specific heat capacity. Water supply and return pipelines are contained in a thermal system network. The thermal network system model is demonstrated in Fig. 4, where T_{\mathrm{s}}^{\mathrm{in}} and T_{\mathrm{s}}^{\mathrm{out}} are the inlet and the outlet water temperature of the water supply pipeline respectively while T_{\mathrm{r}}^{\mathrm{in}} and T_{\mathrm{r}}^{\mathrm{out}} are the inlet and the outlet water temperature of the water return pipeline respectively. There is a basic rule which should be met, that the inflow of water is equal to the outflow of water at each node when the hot water is flowing in the pipeline, which can be expressed by Eq. (32), where A_{\mathrm{n}l}is the node-branch incidence matrix of the thermal system structure (inflow is positive and outflow is negative), mrepresents the water flow through the pipeline, and m_nrepresents the water outflow of each node in the system.

(32)A_{\mathrm{n}l}m=m_n.

In thermal systems, normally the heat load power \Psi can be expressed as Eq. (33), where C_{\mathrm{p}}is the specific heat capacity of the heat medium (normally is water); T_i^{\mathrm{s}} and T_i^{\mathrm{r}} represent the mixed temperature at node i of the water supply and return system respectively. In the water supply and return part, the hot water from each pipeline is mixed at node i. The temperature of the hot water flowing out of node i is equaled to the mixed temperature of this node. Therefore, the temperature value can be expressed by Eqs. (34) and (35), whereU_i^{+}represents the set of pipelines with node i as the front section and U_i^{-} represents the set of pipelines with node i as the end section.

(33)\mathit{\Psi }=C_{\mathrm{p}}{m}_i\left(T_i^{\mathrm{s}}-T_i^{\mathrm{r}}\right),

(34){\sum }_{l\in U_i^{-}}\left(T_{\mathrm{s}}^{\mathrm{out}}m{s}_l\right)=T_i^{\mathrm{s}}{\sum }_{l\in U_i^{-}}m{s}_l,

(35){\sum }_{l\in U_i^{+}}\left(T_{\mathrm{r}}^{\mathrm{out}}m{r}_l\right)=T_i^{\mathrm{r}}{\sum }_{l\in U_i^{+}}m{r}_l.

2.3.3 Outlet water temperature considering temperature loss

In actual situations, due to the limited water flow rate and different lengths of heating pipelines, there will be a delay in the transmission process. However, the dynamic temperature model which considers thermal delay has a lot of nonlinear parts. Therefore, this paper ignores the thermal delay and simplifies the model. The temperature loss of the heat medium in the pipeline is affected by the environment and the transmission distance. Therefore, the outlet temperature of the pipeline considering the temperature loss can be expressed by Eq. (36).

(36)T_l^{\mathrm{out}}=\left(T_l^{\mathrm{in}}-T_{\mathrm{a}}\right){\mathrm{e}}^{\frac{-{\lambda }_{l}L}{C_{p}m{s}_l}}+T_{\mathrm{a}}^{\text{'}},

where T_l^{\mathrm{out}} and T_l^{\mathrm{in}} represent the outlet and the inlet water temperature of pipe l respectively, T_{\mathrm{a}}^{\mathrm{\text{'}}} is ambient temperature, and {\lambda }_l and L are heat transfer coefficient and length of pipeline l.

3 Integrated energy system network and objective function

3.1 Power to gas system

The power-to-gas (PtG) system concept was first proposed in Japan in the 1980s [20]. In recent years, with the increasing share of wind and solar power grow, the PtG technology has attracted a lot of attention from researchers. The PtG technology can achieve large-scale storage of energy through the efficient conversion of electricity and gas to maintain the stability of the system [21]. The structure of the PtG system is exhibited in Fig. 5.

The main process is to electrolyse water with abundant electricity to generate hydrogen, and the generated hydrogen is stored or used to produce methane with the methanation process [22]. Since this paper is focused on the optimal operation of integrated energy systems, here is only a brief introduction to the PtG system.

3.2 Power thermal gas integrated energy system

Power, thermal, and gas subsystems are coupled with each other through the energy conversion device to form an integrated energy system [23]. The energy conversion devices used in this paper include the gas-fired turbine (power and thermal system coupling points), electric boilers (power and thermal system coupling points), and CHP systems (three subsystems coupling point). Different forms of energy are related and transformed into each other in the coupling point of this system.

For the coupled power and gas subsystems, although the power subsystem and gas subsystem has been coupled with each other, the two subsystems still could be analyzed separately [24]. The place where the power generator buses and gas nodes coincide is used to integrate the power and gas network. The generator located in the coupled place is driven by the gas-fired turbine. The compressors are driven by the gas turbine where the gas is tapped from the inlet node. The optimal operation situation can be calculated by optimizing the entire system, while each energy system can meet the loads. The entire system structure in this paper can be found in Fig. 6 in the case study section, which is divided into three subsystems in terms of space when calculating [25].

3.3 Objective function

The objective function of optimal power flow aims to create the biggest social welfare, which means to minimize the entire system cost and maximize the system benefit. Thus the objective function of the optimal power flow in this power and gas coupling system can be expressed as

(37)\min \mathrm{\hspace{0.17em}}{C}_{\mathrm{E}}+C_{\mathrm{G}}+C_{\mathrm{CHP}}-B_{\mathrm{E}}-B_{\mathrm{G}},

where C_{\mathrm{E}}denotes the traditional electrical generation cost, whose equation is a quadratic function. The equation is linearized by the piecewise linearization method to construct a linear model and the process is similar to the linearization process of the natural gas pipeline flow equation. C_{\mathrm{G}} denotes the cost of gas well and C_{\mathrm{CHP}} is the operating cost of the CHP unit. B_{\mathrm{E}}denotes the benefit of electricity consumers and B_{\mathrm{G}}denotes the benefit of gas consumers except the gas-fired generator. The equations of these five terms are expressed as

(38a)C_{\mathrm{E}}={\sum }_{i\notin \mathit{\Psi }}\left[a_{\mathrm{1}}^{\text{'}}+b_{\mathrm{1}}^{\text{'}}{P}_i^{\mathrm{Ele}}+c_{\mathrm{1}}^{\text{'}}{\left(P_i^{\mathrm{Ele}}\right)}^{\mathrm{2}}\right],

(38b)C_{\mathrm{G}}={{\displaystyle \sum }}_{i=\mathrm{1}}^{N_{\mathrm{s}}}{c}_i{\omega }_{S_i},

(38c)C_{\mathrm{CHP}}={{\displaystyle \sum }}_{i=\mathrm{1}}^{N_{\mathrm{CHP}}}{ϵ}_{\mathrm{CHP}}\left({\gamma }_{\mathrm{P}}{P}_{it}^{\mathrm{CHP}}+{\gamma }_{\mathrm{H}}{\mathit{\Psi }}_{it}^{\mathrm{CHP}}\right)+C_{\mathrm{g}},

(38d)B_{\mathrm{E}}={{\displaystyle \sum }}_{i=\mathrm{1}}^{N_{\mathrm{L}}}\left(b_{\mathrm{1}}{P}_i^{\mathrm{Gas}}+c_{\mathrm{1}}{\left(P_i^{\mathrm{Gas}}\right)}^{\mathrm{2}}\right),

(38e)B_{\mathrm{G}}={{\displaystyle \sum }}_{i\notin \mathit{\Psi }}^{N_{\mathrm{G}}}\left(b_{\mathrm{2}}{\omega }_{L_i}+c_{\mathrm{2}}{\omega }_{L_i}^{\mathrm{2}}\right),

where c_irepresents the well-head gas price \$/\left({\mathrm{10}}^{\mathrm{9}}\mathrm{J}\right); {\omega }_{S_i} is the gas supply at node i; {\omega }_{L_i} is the gas demand at node i; \mathit{\Psi } represents the coupled node with gas-fired generators. The value of{ϵ}_{\mathrm{CHP}}, {\gamma }_{\mathrm{P}}, and {\gamma }_{\mathrm{H}} are 24.2, 0.31, and 2.40 in this paper. N_{\mathrm{s}} is the total gas nodes;N_{\mathrm{L}} is the total number of electric consumers; N_{\mathrm{H}}is the total number of heat consumers; a_{\mathrm{1}}^{\text{'}},b_{\mathrm{1}}^{\text{'}},c_{\mathrm{1}}^{\text{'}}represent the coefficients of the non-gas generator cost; and a_{\mathrm{1}},b_{\mathrm{1}},c_{\mathrm{1}}and a_{\mathrm{2}},b_{\mathrm{2}},c_{\mathrm{2}}represent the coefficients of the benefit of the electrical consumer, and the benefit of the gas consumer, respectively.

The total gas volume transformed is calculated by the Weymouth steady-state gas flow equation, and the compressor operation cost is embedded in the gas supply cost at the source node.

4 Methodology

In the optimization calculation process of the entire system, the piecewise linearization method is applied in this paper to facilitate calculation because the gas flow given by Weymouth’s formula is a quartic form [26]. Besides, after the power subsystem and the gas subsystem are calculated separately, the optimization results obtained are calculated using the ADMM method at the coupling point. Finally, the optimization results of the two subsystems at the coupling point tend to be the same under a certain tolerance.

4.1 Gas flow piecewise linearization

Since the gas flow in the pipeline given by Weymouth’s formula is a nonlinear function, the piecewise linearization method is applied in this paper to approximate the function [27]. For the convenience of piecewise linearization, first, the variable S_{ij}^{{\mathrm{G}}^{\mathrm{2}}}is used to represent the square of gas flow value. Then, Eq. (17a) is converted to Eq. (39). When the gas flows from node i to node j in the pipeline, the volume range of the gas flow is \left[\mathrm{0},S_{ij}^{\mathrm{G},\max }\right], where S_{ij}^{\mathrm{G},\max } is the maximum gas flow volume in the pipeline, and {\delta }_{ij}^l is a binary variable whose value is 0 or 1. The abscissa interval of the gas flow value \left[\mathrm{0},S_{ij}^{\mathrm{G},\max }\right]is divided into m parts, the gas flow calculation formula is converted into the linearized form of Eqs. (40)–(45), and the node pressure constraint is also converted into Eq. (46).

(39)C_{ij}^{\mathrm{2}}\left(p_i^{\mathrm{2}}-p_j^{\mathrm{2}}\right)={{\displaystyle \sum }}_{l=\mathrm{1}}^m\left[k_{ij}^l\left(S_{ij}^l-u_{ij}^l\right)+{\delta }_{ij}^l{v}_{ij}^l\right],

(40)\mathrm{0}\le S_{ij}^l\le {\delta }_{ij}^l{u}_{ij}^l\mathrm{\hspace{0.17em}}\mathrm{\hspace{0.17em}}\mathrm{\hspace{0.17em}}\mathrm{\hspace{0.17em}}l=\mathrm{1},

(41){\delta }_{ij}^l{u}_{ij}^{l-\mathrm{1}}\le S_{ij}^l\le {\delta }_{ij}^l{u}_{ij}^{l}l>\mathrm{1},

(42)v_{ij}^l=\mathrm{0}\mathrm{\hspace{0.17em}}\mathrm{\hspace{0.17em}}\mathrm{\hspace{0.17em}}\mathrm{\hspace{0.17em}}\mathrm{\hspace{0.17em}}l=\mathrm{1},

(43)v_{ij}^l={\left(u_{ij}^{i-\mathrm{1}}\right)}^{\mathrm{2}}l>\mathrm{1},

(44)u_{ij}^l={\displaystyle \frac{S_{ij}^{\mathrm{G},\max }l}{m}},

(45)p_i^{{\min }^{\mathrm{2}}}\le p_i^{\mathrm{2}}\le p_i^{{\max }^{\mathrm{2}}}.

4.2 ADMM method applied at the coupling point

The alternating direction of multiplier method (ADMM) is an algorithm that divides the convex optimization problem into smaller parts, and each problem that is broken down into smaller parts will be easier to calculate [28]. In this paper, the power subsystem and the gas subsystem are optimized based on their real-time conditions and produce the optimal value for the coupling point of the two systems respectively [29]. However, the two calculation results for the same coupling point are not equal because the two subsystems are optimized separately. At this time, the ADMM algorithm is used to increase the number of optimization calculation iterations of the two subsystems and narrow down the calculation result gap of the two subsystems. The iteration will stop when the two values tend to be equal within the tolerance, and the entire system will achieve an overall optimal operation.

In this paper, the gas turbine (GE) is selected as the coupling point energy conversion device to couple the power and gas system. When the two subsystems are optimized separately, the estimated value at the coupling point of the power system for the gas system is Bus j,Gand the estimated value of the gas system for the power system at the coupling point is Bus i,P. Meanwhile, the real values optimized by the two subsystems according to their conditions are Busj,P and Busi,G. The structure of the system coupling point is displayed in Fig. 6.

To solve this optimal power flow problem, the ADMM method is used to narrow down the deviation between the estimated value and the actual value by increasing iteration times to obtain the optimal value of the coupling point of the entire system.

4.3 Standard alternating direction multiplier method

The standard alternating direction multiplier method (ADMM) algorithm is used to minimize the augmented Lagrangian function, which is formulated by coupling the optimization problem and constraints through dual variables and adds the penalty term to make the function more convex [30]. Normally, the augmented Lagrangian function and ADMM update iteration form can be written as

L_{\rho }\left(x,z,\lambda \right)=f\left(x\right)+g\left(z\right)+{\lambda }^T\left(Ax+Bz-c\right)

(46)+{\displaystyle \frac{\rho }{\mathrm{2}}}{\left|\left|Ax+Bz-c\right|\right|}_{\mathrm{2}}^{\mathrm{2}},

(47a)x^{k+\mathrm{1}}={\mathrm{argmin}}_x{L}_{\rho }\left(x,z^k,{\lambda }^k\right),

(47b)z^{k+\mathrm{1}}={\mathrm{argmin}}_z{L}_{\rho }\left(x^{k+\mathrm{1}},z,z^k\right),

(47c){\lambda }^{k+\mathrm{1}}={\lambda }^k+\rho \left(A{x}^{k+\mathrm{1}}+B{z}^{k+\mathrm{1}}-c\right),

where \rho >\mathrm{0} is the penalty parameter to ensure the convergence of alternating the direction multiplier method.

The augmented Lagrangian function is formulated as Eqs. (48) and (49) by decomposing the model into the power subsystem cost and the gas subsystem cost, in which P_{it}^g and P_t^{g,ptg} represent the output power of the gas-fired unit i and the power consumption of the PtG unit at time t in the gas subsystem; and {\overline{P}}_{it}^eand {\overline{P}}_{it}^{e,ptg} represent the common variables in coupling power and gas systems. Thus, the decomposed power system function is formulated as Eq. (48) and the gas subsystem function is formulated as Eq. (49) respectively.

\min \mathrm{\hspace{0.17em}}{\sum }_t\{{\sum }_{i\notin \mathit{\Psi }}{C}_{\mathrm{E},t}+\left[{\lambda }_{it}^e\left(P_{it}^e-{\overline{P}}_{it}^e\right)+{\displaystyle \frac{{\rho }_i}{\mathrm{2}}}{\left(P_{it}^e-{\overline{P}}_{it}^e\right)}^{\mathrm{2}}\right]

(48)+C_{\mathrm{CHP},t}^{\mathrm{E}},

\min \mathrm{\hspace{0.17em}}{\sum }_t\{{\sum }_{i\in \mathit{\Psi }}{C}_{\mathrm{G},t}+\left[{\lambda }_{it}^g\left(P_{it}^g-{\overline{P}}_{it}^e\right)+{\displaystyle \frac{{\rho }_i}{\mathrm{2}}}{\left(P_{it}^g-{\overline{P}}_{it}^e\right)}^{\mathrm{2}}\right]

(49)+C_{\mathrm{CHP},t}^{\mathrm{G}}.

The developed 24-h day ahead optimal operation of the integrated energy system containing power thermal and gas subsystems is a quadratic optimization problem which is solved by the piecewise linearization, the ADMM method, and the GUROBI solver in the YALMIP modeling environment in this paper.

The detailed standard ADMM algorithm process is shown in Fig. 7, where k and ε are the iteration times and the convergence thresholds respectively. The iteration process will stop when the deviation between the estimated value and the real value is less than convergence thresholds, and the entire system reaches the optimal situation at that time.

5 Case study and discussion

This case study focuses on the comparison of the 24-h system operation cost and the output of each part of the integrated energy system in different scenarios. The power-thermal-gas integrated energy system includes the IEEE-33 bus power system, the six nodes thermal system, and the 20 nodes natural gas system. The structure of the integrated energy system after coupling is presented in Fig. 8. In this case study, the IEEE-33 bus radial network is used in response to the dist-flow method introduced above. For other power systems with different networks, conventional methods such as Newton’s method can be used to solve the power flow.

In this case study, the generators at nodes 23 and 30 in the IEEE-33 bus power subsystem are both set as the gas turbine, and the natural gas consumed is provided by node 9 and node 19 of the 20 node natural gas subsystem. The generator at the 15th bus is generated by the CHP unit, which is also connected at node1 of the 6-node thermal subsystem as a heat source. The 20-node natural gas subsystem contains 6 natural gas sources and 9 gas loads. The heat source in the thermal subsystem includes a gas-fired CHP unit and an electric boiler (EB). The energy consumed by the CHP and the EB is provided by the natural gas subsystem node 3 and the power subsystem bus 17 respectively.

The data of the IEEE-33 node power subsystem are taken from the standard data in the MATPOWER. The natural gas subsystem is a Belgian 20-node natural gas system. The basic parameters of the gas source are listed in Table A1, and the output cost of the natural gas source is set to 0.75$/m³ (1–3 gas sources) and 0.65$/m³ (4–6 gas sources) respectively. The maximum output of the gas turbine unit is 25 MW in this case study.

The thermal subsystem includes the CHP unit, the electric boiler, and the 3 thermal loads. The total load is 5 MW, the pyroelectricity ratio coefficient of the electric boiler η is taken as 0.8, and the basic parameters of the CHP unit are {ϵ}_{\mathrm{CHP}}=\mathrm{24.2},{\gamma }_{\mathrm{P}}=\mathrm{0.31} and {\gamma }_{\mathrm{H}}=\mathrm{2.40}. The power, thermal, and gas load curve is manifested in Fig. 9. Due to the large gap between the thermal load and the other two loads, the thermal load curve in Fig. 9 is 100 times as large as the real thermal load value.

Based on this integrated energy system model, three scenarios are considered in which the thermal subsystem load is met by the natural gas subsystem, the thermal subsystem load is met by the power subsystem, and the thermal subsystem load is shared by the power and gas subsystem. The output of each part of the integrated energy system and the overall operation cost are optimized and analyzed, and the output of power system generators, the gas turbines, the CHP unit, and the electric boilers in different scenarios is shown in Fig. 10.

The output of power system generators in three scenarios is shown in Fig. 10. Since natural gas is used as the fuel for gas turbines and the gas-fired CHP unit in the power subsystem and thermal subsystem, the output of gas sources is different in different scenarios. The gas source output in different scenarios obtained after optimization of this case study is shown in Fig. 11.

For the gas turbine, the CHP unit, and the electric boiler, the system thermal load is met by the gas-fired CHP unit in Scenario 1, the system thermal load is met by the electric boiler in Scenario 2. In Scenario 3, due to the operation of the gas turbine, the thermal load is met by the electric boiler and the CHP system together. Therefore, when the integrated energy system reaches its optimal operation, the gas turbine is only used in Scenario 3, the CHP unit is not used in Scenario 2, and the electric boiler is not used in Scenario 1. In Scenarios 1 and 3, the electricity output will be increased based on meeting the thermal demand. Meanwhile, the operation of the CHP unit is restricted by its feasible region and the output of the CHP will be adjusted to the most efficient operating point in the feasible region. The electric power output and thermal power output of the CHP unit, the thermal power of the electric boiler as well as electric power output of the gas turbine are shown in Figs. 12, 13, and 14.

Table 1 summarizes and compares the overall operation cost of the integrated energy system and the output of each part in the three scenarios. The results in Table 1 indicate that compared to using the power subsystem and the gas sub-system to meet the thermal load of the system individually, the integrated energy system that couples the power, thermal, and gas subsystems has the lowest overall cost.

6 Conclusions

A steady-state model of the natural gas subsystem has been established by analogy with the power subsystem model and applying the dist-flow method to quickly solve the power flow value of the power subsystem in this paper. In addition, the method of piecewise linearization is applied to solve the Weymouth quadratic equation. Moreover, the ADMM algorithm is applied to narrow down the deviation between the different optimal results to the tolerance and calculate the optimal value at the coupling point of the power and gas subsystem through multiple iterations to make sure that the entire system can reach its optimal operating state. The power, thermal, and gas integrated energy system can coordinate various energy distributions while meeting various load requirements and enhance the flexibility and stability of the system as well as reducing system operation costs to some extent.

Future research will focus on the establishment of transient models for natural gas flow and the privacy protection of each subsystem in the integrated energy system.

References

[1] An S, Li Q, Gedra T W. Natural gas and electricity optimal power flow. In: 2003 IEEE PES Transmission and Distribution Conference and Exposition, Dallas, TX, USA, 2003

[2] Baran M, Wu F F. Optimal sizing of capacitors placed on a radial distribution system. IEEE Transactions on Power Delivery, 1989, 4(1): 735–743

[3] Chen S, Wei Z, Sun G, Identifying optimal energy flow solvability in electricity-gas integrated energy systems. IEEE Transactions on Sustainable Energy, 2017, 8(2): 846–854

[4] Clegg S, Mancarella P. Integrated modeling and assessment of the operational impact of power-to-gas (P2G) on electrical and gas transmission networks. IEEE Transactions on Sustainable Energy, 2015, 6(4): 1234–1244

[5] Erseghe T. Distributed optimal power flow using ADMM. IEEE Transactions on Power Systems, 2014, 29(5): 2370–2380

[6] Guo J, Hug G, Tonguz O K. A case for nonconvex distributed optimization in large-scale power systems. IEEE Transactions on Power Systems, 2017, 32(5): 3842–3851

[7] He C, Wu L, Liu T, Robust co-optimization scheduling of electricity and natural gas systems via ADMM. IEEE Transactions on Sustainable Energy, 2017, 8(2): 658–670

[8] Kantor I, Robineau J L, Bütün H, A mixed-integer linear programming formulation for optimizing multi-scale material and energy integration. Frontiers in Energy Research, 2020, 8: 49

[9] Khani H, Farag H E Z. Optimal day-ahead scheduling of power-to-gas energy storage and gas load management in wholesale electricity and gas markets. IEEE Transactions on Sustainable Energy, 2018, 9(2): 940–951

[10] Ko W, Kim J. Generation expansion planning model for integrated energy system considering feasible operation region and generation efficiency of combined heat and power. Energies, 2019, 12(2): 226

[11] Lejarza F, Baldea M. Economic model predictive control for robust optimal operation of sparse storage networks. Automatica, 2021, 125: 109346

[12] Li C, Dong Z, Chen G, Data-driven planning of electric vehicle charging infrastructure: a case study of Sydney, Australia. IEEE Transactions on Smart Grid, 2021, 12(4): 3289–3304

[13] Li C, Liu C, Deng K, Data-driven charging strategy of PEVs under transformer aging risk. IEEE Transactions on Control Systems Technology, 2018, 26(4): 1386–1399

[14] Li C, Xu Y, Yu X, Risk-averse energy trading in multienergy microgrids: a two-stage stochastic game approach. IEEE Transactions on Industrial Informatics, 2017, 13(5): 2620–2630

[15] Li C, Yu X, Huang T, Distributed optimal consensus over resource allocation network and its application to dynamical economic dispatch. IEEE Transactions on Neural Networks and Learning Systems, 2018, 29(6): 2407–2418

[16] Liu C, Shahidehpour M, Wang J. Coordinated scheduling of electricity and natural gas infrastructures with a transient model for natural gas flow. Chaos (Woodbury, N.Y.), 2011, 21(2): 025102

[17] Liu N, Tang Q, Zhang J, A hybrid forecasting model with parameter optimization for short-term load forecasting of micro-grids. Applied Energy, 2014, 129: 336–345

[18] Liu T, Zhang W, Meng K, A power-to-gas integrated microgrid optimal operation strategy based on rolling horizon. In: 2019 IEEE Milan PowerTech, Milan, Italy, 2019

[19] Luo S, Yang L, Zhang X, A fully linear-constrained optimal electricity-gas flow in an integrated energy system. In: 2018 2nd IEEE Conference on Energy Internet and Energy System Integration (EI2), Beijing, China, 2018

[20] Mühlpfordt T, Dai X, Engelmann A, Distributed power flow and distributed optimization—formulation, solution, and open source implementation. Sustainable Energy. Grids and Networks, 2021, 26: 100471

[21] Qiu J, Dong Z, Zhao J, Multi-stage flexible expansion co-planning under uncertainties in a combined electricity and gas market. IEEE Transactions on Power Systems, 2015, 30(4): 2119–2129

[22] Qiu J, Zhao J, Yang H, Optimal scheduling for prosumers in coupled transactive power and gas systems. IEEE Transactions on Power Systems, 2018, 33(2): 1970–1980

[23] Rubio R, Ojeda-Esteybar D, Ano O, Integrated natural gas and electricity market: a survey of the state of the art in operation planning and market issues. In: 2008 IEEE/PES Transmission and Distribution Conference and Exposition: Latin America, Bogota, Colombia, 2008

[24] Shabanpour-Haghighi A, Seifi A R. Energy flow optimization in multicarrier systems. IEEE Transactions on Industrial Informatics, 2015, 11(5): 1067–1077

[25] Sun L, Jin Y, Shen J, Sustainable residential micro-cogeneration system based on a fuel cell using dynamic programming-based economic day-ahead scheduling. ACS Sustainable Chemistry & Engineering, 2021, 9(8): 3258–3266

[26] Sun L, Li G, You F. Combined internal resistance and state-of-charge estimation of lithium-ion battery based on extended state observer. Renewable & Sustainable Energy Reviews, 2020, 131: 109994

[27] Unsihuay C, Lima J W M, de Souza A C Z. Modeling the integrated natural gas and electricity optimal power flow. In: 2007 IEEE Power Engineering Society General Meeting, Tampa, FL, USA, 2007

[28] Wang C, Dong S, Xu S, Impact of power-to-gas cost characteristics on power-gas-heating integrated system scheduling. IEEE Access: Practical Innovations, Open Solutions, 2019, 7: 17654–17662

[29] Wang J, Xu J, Liao S, Coordinated optimization of integrated electricity-gas energy system considering uncertainty of renewable energy output. Automation of Electric Power Systems, 2019, 43(15): 2–9 (in Chinese)

[30] Zhou Q, Shahidehpour M, Li Z, Two-layer control scheme for maintaining the frequency and the optimal economic operation of hybrid AC/DC microgrids. IEEE Transactions on Power Systems, 2019, 34(1): 64–75

RIGHTS & PERMISSIONS

Higher Education Press