Delft University of Technology
A graph-based framework for steady-state load flow analysis of multi-carrier
Markensteijn, A.S.; Romate, Johan; Vuik, C.
Publication date 2019
Document Version Final published version Citation (APA)
Markensteijn, A. S., Romate, J., & Vuik, C. (2019). A graph-based framework for steady-state load flow analysis of multi-carrier. (Reports of the Delft Institute of Applied Mathematics; Vol. 19-01). Delft University of Technology.
Important note
To cite this publication, please use the final published version (if applicable). Please check the document version above.
Copyright
Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy
Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.
DELFT UNIVERSITY OF TECHNOLOGY
REPORT 19-01
A graph-based framework for steady-state load flow analysis
of multi-carrier energy networks
A.S. Markensteijn, J.E. Romate, C. Vuik
ISSN 1389-6520
Copyright c
2019 by Delft Institute of Applied Mathematics, Delft, The Netherlands.
No part of the Journal may be reproduced, stored in a retrieval system, or transmitted, in
any form or by any means, electronic, mechanical, photocopying, recording, or otherwise,
without the prior written permission from Delft Institute of Applied Mathematics, Delft
A graph-based framework for steady-state load flow analysis
of multi-carrier energy networks
A.S. Markensteijn, J.E. Romate, C. Vuik
March 14, 2019
Abstract
Energy systems are becoming more complex due to increased coupling between different networks, resulting in multi-carrier energy networks. Conventional models for the separate networks are not able to capture the full extend of the coupling. Recently, different models for multi-carrier networks have been proposed, either using the energy hub concept or using a case specific approach. Although the energy hub concept can be applied to a general integrated net-work, it is unclear how the energy hub should be represented in the graph of the multi-carrier network. This paper presents a graph-based framework for steady-state load flow models of multi-carrier energy systems. Furthermore, the effect of coupling on the integrated system of equations is investigated. The proposed framework is tested on two small multi-carrier net-works, for comparison with models in literature. Results show that our framework is applicable to a general system, and that it generalizes both the energy hub concept and the case specific approaches.
1
Introduction
Multi-carrier energy systems (MES) have become more important over the years as the need for efficient, reliable and low carbon energy systems increases. In MES, different energy carriers, such as electricity and heat, interact with each other leading to one combined system. They have higher performance than classical single-carrier energy systems due to increased flexibility, reliability, use of renewables and distributed generation, and reduced carbon emission. Because these multi-carrier energy systems integrate two or more separate energy systems they are sometimes called integrated energy systems. An overview of MES is given by Mancarella [2014].
An important tool for designing and operating energy systems is steady-state energy flow anal-ysis (also called power flow or load flow analanal-ysis) of the energy transportation and distribution networks. This analysis determines the network state parameters and the flow of energy in the network. Load flow models for single-carrier network have been well studied. Recently, different models have been proposed to model load flow for MES. Two types of models can be distinguished. The first uses the energy hub concept, the other a more ad hoc approach.
The energy hub concept was first introduced by Geidl and Andersson [2007] to model the rela-tion between different energy carriers, and to optimize MES. The input and output energy of the hub are related through a coupling matrix. Unidirectional flow from input to output is assumed,
such that the coupling matrix is constant or a function of the input power only. Within the energy hub, the transmission of energy carriers is not taken into account. Wasilewski [2015] and Long et al. [2017] extend the energy hub concept to allow bidirectional flow, based on graph and network theory. They provide different detailed graph representations of energy flow within an energy hub. However, connecting the energy hub to the rest of the energy network is not described, such that the network state parameters in the single-carrier parts of the network are not modeled. Ayele et al. [2018] extend the energy hub concept by explicitly modeling the connection between the en-ergy hub and the rest of the enen-ergy network. All local enen-ergy generation is included in the hub. Their model determines network state parameters in the single-carrier parts. However, the out-put power of the hub is assumed to be a linear function of the inout-put power. Moreover, the energy hub is not seen as a node or an edge, such that the graph representation of the energy hub is unclear. The second type of model combines the existing model equations of the single-carrier networks into one system of equations. This way of modeling provides both the energy flow through a net-work and the netnet-work state parameters, and allows the use of detailed models for the conversion and transmission of energy. An et al. [2003] introduced such a model for a combined gas and elec-tricity network, and Liu et al. [2016] introduced one for a combined heat and elecelec-tricity network. Martinez-Mares and Fuerte-Esquivel [2012] introduced a model for a combined electricity and gas network, using a distributed slack node approach in the electricity network, and taking into account temperature effects in the gas network. Pan et al. [2016] introduced a model for a combined elec-tricity and heat network, taking the different time scales in the heat network into account. More recently, Shabanpour-Haghighi and Seifi [2015] and Abeysekera and Wu [2015] both gave models for a combined gas, electricity, and heat network. Furthermore, Liu and Mancarella [2016] extended the previous work of Liu et al. [2016] to also include gas. However, these models are case specific and are therefore difficult to apply to a general integrated energy network.
Jalving et al. [2017] proposed a general graph-based computational framework for optimizing com-bined networks, which they apply to a comcom-bined electricity and gas network. The load flow equa-tions are part of the constraints, such that steady-state load flow analysis is considered a special case of their proposed optimization model. They assume that a feasible solution, satisfying these constraints, exists. However, some coupling models and single-carrier load flow equations can lead to a combined system which has no (feasible) solution.
To the best of the authors’ knowledge the current load flow models provide a case specific ap-proach for a load flow problem of MES, or a general model for the conversion of one carrier into another. Most do not state how the graphs of single-carrier networks can be combined into one multi-carrier network. Moreover, they do not consider the effect of this combination on the load flow model for a MES. Usually, the coupling models introduce more unknowns than equations to the system. Therefore, additional equations or boundary conditions must be added for the system to be (uniquely) solvable. Depending on how the single-carrier networks are connected to form an integrated network, and depending on the choice of the additional equations, the resulting sys-tem of load flow equations might have none, one, or infinite solutions. A syssys-tematic analysis of the single-carrier load flow models and graphs, and of possible graph representations of the multi-carrier system, is important to identify and understand these ill-posed systems. Therefore, we propose a general load flow model framework for multi-carrier energy networks, based on graph and network theory.
We propose a graph representation for MES, and provide guidelines for combining the single-carrier load flow models with conversion models, into one multi-single-carrier load flow model. The graph representation of single-carrier and multi-carrier networks is presented in Section 2. The integrated load flow model is discussed in Section 3. A case study consisting of a gas network, power grid (electrical network), and a heating network is provided in Section 4. This case study is based on the ones by Shabanpour-Haghighi and Seifi [2015] and Ayele et al. [2018] for comparison. Furthermore, it provides insight into the difficulties of combining single-carrier load flow models and conversion models. Some concluding remarks are given in Section 5.
2
Graph-representation
Energy networks all have their own terminology. For instance, an electrical network is usually called a power grid. This specific terminology also extends to the elements of a network. For instance, a basic power grid consists of power lines connected by buses, whereas a gas pipeline network consists of pipelines connected by junctions. Mathematically, all networks are abstracted to a graph. Some terms and definitions of graphs and energy networks are given below.
2.1
Terms and definitions
An (undirected) graph G is a pair (V, E ) where V is a set of nodes or vertices vi and E is a set of
links or edges ek. An edge is itself a set of two nodes such that ek= {vi, vj}. We denote the size of
a set S by |S|. A directed graph G is a pair (V, E ) where V is again a set of nodes. E is a set of arcs ek, where an arc is an ordered pair of nodes ek= (vi, vj). Thus, an arc is a link with a predefined
direction. The actual direction of flow through the link could be opposite to this direction. Arcs are also referred to as edges or links.
Inflow and outflow of a network are usually represented by special nodes called sources and sinks respectively. However, this representation is difficult to use for energy networks. First, it is some-times impossible to know up front if a node represents an inflow or outflow of energy. Secondly, in modeling it is convenient to see the inflow and outflow of a node as a flow through a link. Therefore, we introduce terminal nodes and terminal links.
A terminal node is a node that can act as both a source or a sink, depending on the direction of the terminal link. A terminal link is a special type of link that has no physical model, and is only connected to one node. It is sometimes called a half-edge and can be denoted by tl= {vi}. A
terminal link is simply a representation of flow entering or leaving the network. Hence, by defini-tion, a terminal link can only be connected to a terminal node. Moreover, connecting a terminal link to any node turns that node into a terminal node. Conversely, a node without a terminal link connected to it is not a terminal node. One node can have more than one terminal link connected to it. We denote the set of terminal links in the graph by T , and the set of terminal links connected to a node v ∈ V by T (v).
A network can then be represented as a collection of nodes, (directed) links, and terminal links, such that N = {V, E , T } = {G, T }. If G is a directed graph the network is said to be directed, and it is said to be undirected if G is an undirected graph.
One property of graphs that is used for modeling energy networks is the incidence matrix. For an undirected graph G = (V, E ) the elements of the |V| × |E | incidence matrix A are defined as
Aik=
(
1, if ek = {vi, vj}
0, otherwise (1)
for every link ek ∈ E and every node vi ∈ V. Similarly, for a directed graph the elements of the
incidence matrix are defined as
Aik= 1, if ek= (vj, vi) −1, if ek = (vi, vj) 0, otherwise (2)
for every link ek ∈ E and every node vi∈ V.
2.2
Single-carrier energy systems
Usually, a gas network and a heat network are represented by a directed graph, whereas a balanced AC power grid is represented by an undirected graph. The representation of the heat network as a directed graph is less straightforward and is explained below in more detail.
The physical pipeline system of a heat network consists of a supply part (supplying warm wa-ter to demands) and a return part (returning cold wawa-ter to sources). These two parts are connected to each other through the (heat) loads and sources. Heat is then injected or consumed through heat exchangers Frederiksen and Werner [2014]. Figure 1a gives a model representation of a source and a load connected with pipes. This means that the hydraulic part of the heat network is a closed system. For most classical heat networks it can be assumed that the water flow through the return lines is opposite in direction, but equal in size, to the water flow in the supply line [Kuosa et al., 2013]. This means that the return line does not have to be modeled explicitly, so that a source connected with a pipe to a load can be abstracted to a directed graph, with terminal links, as shown in Figure 1b. Note that in this representation the (hydraulic part) of the heat network is no longer a closed system.
For every network, variables are associated with the links and nodes in the graph. For basic steady-state load flow, these are pressure p and flow q for a gas network, voltage V , current I, and complex power S for a power grid, and head h, water flow m, heat power ϕ, supply temperature Ts,
return temperature Tr, and outflow temperature Tofor a heat network. Table 1 gives an overview
of these variables and the graph element they are associated with. The variables associated with a terminal link can be seen as variables of the terminal node. To distinguish between (terminal) link and terminal node variables, the ones on the terminal node are called injected. If a node has more than one terminal link connected to it, then the injected flow or power is the sum of the flows and powers of the terminal links.
For the basic load flow equations (see Section 3) it is convenient to associate some of the link and node variables with the beginning or end of a link, directly next to a node. Figures 1b, 2a, and 2b show the graph representation for each of the three single-carrier networks. They show a source
ms ij Ts j ϕj ϕi mj To j Tijs,start Tijr,end To i Tr i mi mr ji load source j Ts i Tr i Ts j Tr j Tijr,start Tijs,end j i i
(a) Model representation.
m
i, T
ioϕ
im
j, T
joϕ
jm
ijT
ijr,endT
ijs,startT
s,end ijT
s j, T
jrT
s i, T
iri
j
T
ijr,start (b) Graph representation.Figure 1: A heat network with a source, represented by node i, connected to a load, represented by node j, using one heat pipe. (a) shows a more realistic model representation with both the supply and return line. (b) shows the graph representation with only the supply line, represented by an arc, and the load and source, represented by terminal links. T is the temperature, m the mass flow, and ϕ the heat power.
Table 1: Variables for a gas, heat, and electrical network, per graph element they are associated with.
Network Node Link Terminal node
Gas pressure p flow q injected flow q
Heat head h flow m injected flow m
supply temperature Ts outflow temperature To
return temperature Tr heat power ϕ
Electricity voltage V current I injected current I
i j pj qij
pi
qi qj
(a) Gas network.
Iij Sij Sji Iji j i Sj Vj Vi Si Ii Ij (b) Power grid.
Figure 2: Graph representation of a source, represented by node i, connected to a load, represented by node j, by one link. The variables used in load flow analysis are shown at the location in the graph where they are modeled. (a) shows a gas network with p pressure and q volume flow. (b) shows a power grid with I current, V voltage, and S complex power.
connected with one link to a sink, and the variables at the location in the graph where they are modeled. A line in a power grid (see Figure 2b) does not have one single link current, but one at each side of a link. For most transmission line models it holds that Iij 6= −Iji, whereas qij = −qji
and mij = −mji for pipelines.
2.3
Coupling of energy systems
Wasilewski [2015] proposes a detailed graph representation of an energy hub, allowing for bidirec-tional flow. It uses directed edges and terminal nodes to form energy converter graphs and energy storage graphs. Long et al. [2017] also extend the energy hub concept to allow bidirectional flow by introducing a MES node, consisting of an energy conversion model, a storage model, and a prosumer model. The conversion model is represented by a graph consisting of terminal nodes, sum nodes, and transmitter nodes. Nodes are associated with an energy carrier, such that conversion of one energy carrier to another is represented by a directed edge. However, both only consider energy flow and do not describe how the energy hub should be connected to the rest of the network. Ayele et al. [2018] connect the energy hub with the rest of the network using a node called the point of interconnection. The energy hub, consisting of demand, local generation, edges, and a coupling system, is not seen as a node or an edge. However, a graph is a collection of nodes and links, which holds for the abstraction of any energy system, both single-carrier and multi-carrier. Coupling single-carrier networks to form a multi-carrier network can therefore only be done through (a combination of) nodes or links.
There are three main options to couple two nodes: connect the two nodes by a link, merge the two nodes into one node, or introduce an additional node and connect the two nodes to it. Since a link is a component that has two flow connections, it is difficult to use a link to model a coupling that involves more than two energy carriers such as a combined heat and power plant. Furthermore, the physical interpretation of a coupling link is not straightforward. Coupling two nodes by merg-ing them into one is difficult because of the parameters associated with the nodes. For instance, suppose we want to connect two gas networks by merging one node of each network with each other into one new gas node. Both nodes have a pressure associated with them, which means some combined pressure must be defined for the new gas node, or gas nodes with multiple pressures must be allowed. When we want to couple two networks of different energy carriers this problem becomes even more clear. Additionally, merging two nodes of a different energy carrier would require some adaptation of Kirchhoff’s law, for instance to distinguish between the links connected to the node (see Section 3 for more details on Kirchhoff’s law). Therefore, we choose to couple networks by
power
grid
gas network
i
heat network
h
T
s,endi
ci
gi
eQ
cT
r,startT
r,endP
cq
cT
s,startm
cT
o,cϕ
cFigure 3: Coupling node ic (gray), connected by dummy links to node ig of a gas network (solid green), node ie of a power grid (dashed red), and node ih of a heat network (dotted-dashed blue). The coupling parameters are shown next to the (terminal) links they are associated with. Here, q is gas volume flow, P is active power, Q is reactive power, T is temperature, m is water mass flow, and ϕ heat power.
introducing an additional node, called a coupling node.
For our model framework, no parameters are associated with a coupling node. Therefore, the coupling node does not belong to any of the single-carrier networks. We distinguish a homogeneous coupling node, which couples networks of the same energy carrier, and a heterogeneous coupling node, which couples networks with different energy carrier. Nodes and links of a single-carrier net-work are also said to be homogeneous. A netnet-work with one or more heterogeneous nodes is called heterogeneous, and conversely, a network without any heterogeneous node is called homogeneous. We focus on coupling networks of different energy carriers. A heterogeneous coupling node can have (terminal) links of different energy carriers connected to it. In theory, this could be any link, but since we do not define parameters on the coupling node, some links can not be connected to it. For instance, a gas link representing a pipeline has a (steady-state) flow equation associated with it, which gives a relation between the gas flow through the pipe and the pressures of the start and end node of the link (see Section 3 for more details). Since the coupling node has no pressure associated with it, such a gas link cannot be connected to it. We choose to connect the coupling node to homogeneous (single-carrier) nodes by dummy links. These links do not represent any physical component, they merely show a connection between nodes. However, they are ho-mogeneous links and as such have the same parameters associated with them as any hoho-mogeneous link. The graph representation of a coupling node, connected to a single node of a gas network, a power grid, and a heat network is shown in Figure 3. The direction shown on the gas and heat dummy links and on the heat terminal links are the predefined directions of the links, which could be different from the actual direction of flow. Hence, the coupling node allows for bidirectional flow. The coupling node can be used to represent something as relatively simple as a junction node
in gas network, or something as complex as a complete (heterogeneous) network. Representing another network, for instance the energy hub graph introduced by Wasilewski [2015], will lead to an hierarchical network. Such an hierarchical approach is used by Jalving et al. [2017], in their proposed computation framework. Furthermore, the proposed coupling node only needs to have inflow and outflow that can be represented by dummy links to connect it to the rest of the network. Due to the choice for dummy links, the integrated network could be separated easily into single-carrier parts by cutting the dummy links into two terminal links. This could be used to solve the integrated network in a decoupled way. This framework for coupling single-carrier networks allows great flexibility when connecting multiple networks, making it applicable to general multi-carrier energy networks.
3
Load flow models
The steady-state energy flow problem (also called power flow or load flow problem) is formulated by collecting the model equations of different physical components, which are in turn related to elements in the graph. There are two equations that hold for any (single-carrier) energy network: Kirchhoff’s first and second laws. For every node, Kirchhoff’s first law states that the total in- and outflow must sum up to zero. In a power grid this law is usually called Kirchhoff’s current law, whereas in a gas and heat network it is conservation of mass or conservation of flow. For every loop in a network, Kirchhoff’s second law states that the sum of potential differences over every link must be zero. In a power grid this law is usually called Kirchhoff’s voltage law, whereas in a gas and heat network it is called the loop pressure equation. Additionally, most links are modeled by an equation that gives a relation between the link variable and the nodal variables. Different models are used, depending on what physical component is associated to the link. For instance, Ohm’s law is used for a transmission line in a power grid. Some networks have additional node or link equations, such as the complex power equation in a power grid, which gives a relation between total injected current, total injected complex power, and nodal voltage.
3.1
Gas network
Assuming a steady-state approximation, a basic gas network can be completely described by a combination of the conservation of mass and a link equations that gives the relation between nodal pressures and link flow [Osiadacz, 1987]. For a general gas network, the conservation of flow is given by a linear system of equations:
Agq − qinj= 0 (3)
with Ag the |Vg| × |Eg| incidence matrix of the gas network (2), q the vector of length |Eg| of link
flows, and qinj the vector of length |Vg| of injected flows. A link equation is then used to give the relation between nodal pressures and link flow:
fkg(qk, pi, pj) = 0 (4)
for every link egk ∈ Eg. These link equations are generally non-linear in the pressures. Combining
|Vg| of nodal pressures, gives the following non-linear system of equations: Fg(xg) =A gq − qinj fg(q, p) = 0 with xg= qinj q p (5)
This system consists of |Vg| + |Eg| equations and 2|Vg| + |Eg| unknowns. To reduce the number
of unknowns, some nodal variables are assumed known. We will refer to this as the boundary conditions of the network. For basic load flow models, either the pressure or the total injected flow is given at every node. Furthermore, at least one of the equations of the conservation of mass (3) is a linear combination of the other equations. This means that if the injected flow is given at every node, the system of equations will generally not have a solution. If the pressure is not given for at least one node, there are infinite solutions, if a solution exists at all, because the link equations (4) are functions of the pressure drop. Hence, at one node the pressure is given while the injected flow is not known. This pressure is called the reference pressure and this node is called the reference node. A node where the injected flow is given while the pressure is not known is called a load node. Nodes where both the pressure and the injected flow are given, or nodes where neither are given, can also occur. We call them a reference load node and a slack node respectively. The system of equations can be reformulated in different ways, based on for instance algebraic substitutions. The most common ones are the nodal formulation, the loop formulation, and the nodal-loop formulation. See for instance Osiadacz [1987] for details, advantages, and disadvantages of each formulation.
3.2
Power grid
Assuming an AC steady-state approximation, a power grid is completely described by a combination of Kirchhoff’s current law, a link equation that gives a relation between nodal voltage and link current (e.g. Ohm’s law), and the complex power equation [Schavemaker and Van der Sluis, 2008]. For a general power grid, Kirchhoff’s current law is given by a linear system of equations:
AeI − Iinj= 0 (6)
with Ae the incidence matrix of the power grid (1), I the vector of complex currents, and Iinj the vector of injected complex currents. A link equation is then used to give the relation between nodal voltages Vi and Vj, and link current Ik:
fke(Ik, Vi, Vj) = 0 (7)
for every link ee k ∈ E
e. Finally, the complex power equation is given by a non-linear equation in I
and V :
Sinj= V Iinj∗
(8) where [·]∗denotes the complex conjugate. Combining Kirchhoff’s current law (6), the link equations (7), and the complex power equations (8) gives the following non-linear system of equations:
Fe(xe) = AeI − Iinj fe(I, V ) Sinj− V Iinj∗ = 0 with x e= Iinj I V Sinj (9)
Because all parameters associated with an AC power grid are complex, this system consists of 4|Ve| + 2|Ee| equations and 6|Ve| + 2|Ee| variables. The complex power S is usually split in its real
and imaginary part. The real part is called the active power P , and imaginary part the reactive power Q. The voltage V is split in its angle δ and amplitude |V |. Again, boundary conditions are introduced to reduce the number of unknowns. For classic load flow, two variables are given at each node: Pinj and Qinjfor a load node, Pinj and |V | for a generator node, or |V | and δ for a reference node. The reference node is introduced for similar reasons as for the gas network; for at least one node the injected complex power must be unknown, and for at least one node the voltage must be given. Mathematically, any other combination of known or unknown nodal variables could be used, for instance a Qδ-node with known injected reactive power Q and nodal voltage angle δ. However, from a technical perspective, only the first three nodes are the feasible node types [Schavemaker and Van der Sluis, 2008].
The system of equations with boundary conditions can be reformulated in different ways, but the most common is the complex power formulation (using polar coordinates) (e.g. Schavemaker and Van der Sluis [2008], Idema et al. [2010]). For other possible formulations, see for instance [Stott, 1974] or [Sereeter et al., 2017].
3.3
Heat network
A heat network consist of a hydraulic part and a thermal part. Assuming a steady-state approxi-mation, the hydraulic part is completely described by a combination of conservation of mass and a link equation. The thermal part is completely described by a nodal mixing-rule, a link equation to relate mass flow to temperatures, and a heat power equation [Bordin et al., 2016, Liu et al., 2016]. For a general heat network, the conservation of mass is given by a linear system of equations:
Ahm − minj= 0 (10)
with Ah the incidence matrix of the heat network (2), m the vector of mass flows, and minj the vector of injected mass flows. As for gas, a link equation is used to give the relation between nodal heads hi and hj, and link mass flow mk:
fkh(mk, hi, hj) = 0 (11)
for every link ehk ∈ Eh. Again, these link equations are generally non-linear. The conservation of
mass (10) combined with the link equations (11) gives the hydraulic model for the heat network. The thermal model governs the temperatures in the supply lines, the return lines, and directly after a source or load, and the heat power injected by a source or consumed by a sink. We will assume that nodes are a source, a sink, or a junction. A source is a node with water inflow and heat power injection, and conversely a sink (or load) has water outflow and heat power consumption. A junction is a node in the network where the flow is redistributed, there is no connection between return and supply line, such that there is no water in- or outflow. For a component connected to a source or sink (that is, for a terminal link) the heat power equation holds:
fi,lϕ ϕi,l, mi,l, Tis, T r i, T
o
i,l = 0 (12)
for every terminal link tl∈ Th(vi) connected to node vi∈ Vh. Here, ϕi,lis the injected or consumed
vi and Tir its return temperature, and Ti,lo is the temperature directly after the source or sink at
terminal link tl. The total injected mass flow of node vi is then given by minji =
P
tl∈T (vi)
ml,i. At
every node vi ∈ Vh, the temperature in the supply and return line are determined by the mixing
rule, which states that the nodal temperature is the weighted average of the inflow temperatures: fiTs =XmsoutTis−XmsinTins= 0
fiTr=XmroutTir−XmrinTinr= 0
(13)
Here,P ms
out is used to denote the sum of all the outgoing flow of node i in the supply line, both
on the edges and on the terminal links. Similarly,P ms
in denotes the sum of all ingoing flows of
node i in the supply line,P mr
out the sum of all outgoing flows of node i in the return line, and
P mr
in the sum of all ingoing flows of node i in the return line. It holds thatP m r
in=P m s out and
P mr
out =P msin. For every supply and every return pipeline, a link equation is used to give the
relation between the temperatures at the beginning and end of the pipe and the pipe mass flow: fkψ Tkstart, T
end
k , mk = 0 (14)
The ‘start’ and ‘end’ of the pipeline are defined for the actual direction of flow, not with respect to the predefined direction of the directed link. Hence, Tks,start= Tis if the flow through pipe k is from node i to node j, and Tks,start= Ts
i if the direction of flow is in the opposite direction.
The heat power equation (12) combined with the mixing rule (13) and the link equation (14) gives the thermal model for the heat network. The hydraulic model and thermal model can be combined to one hydraulic-thermal model, given by the following, generally non-linear, system of equations: Fh xh = Ahm − minj fh(h, m) fϕ(ϕl, ml, Ts, Tr, To) fTs(m, m l, Ts, To) fTr(m, m l, Tr, To) fψ Ts,start, Ts,end, m fψ Tr,start, Tr,end, −m = 0 with xh= ml m h Ts Ts,end Tr Tr,end To ϕ (15)
This system consists of 3|Vh| + 3|Eh| + |Th| equations and 3|Vh| + 3|Eh| + 3|Th| variables. Boundary
conditions are introduced to reduce the number of unknowns. Generally, at each source or sink node (at each terminal link) two variables are given: ϕl and Tlo. These nodes are called a source
node and a load node respectively. Because a circulation pump is usually located at a source, at one of the source nodes Ts and h are given. This head is then the reference head, and since the
injected heat power is not known, this node is called a source reference slack node. We assume such a node has only one component (i.e. one terminal link) connected to it. Otherwise, for instance ϕl
for much the same reasons as for the gas network and the power grid; for at least one node the in-jected heat power and mass flow must be unknown, and for at least one node the head must be given. Again, the system of equations with boundary conditions can be reformulated in different ways. See for instance Liu et al. [2016] who consider the difference between the hydraulic-thermal model and separate hydraulic and thermal models, or Arsene et al. [1989] who compare the nodal and loop formulation of the hydraulic model.
3.4
Multi-carrier energy systems
3.4.1 Model
We couple homogeneous node(s) to a heterogeneous coupling node using dummy links. The dummy links do not represent any physical component, they only show a connection between nodes. How-ever, they could be seen as lossless link, which means we could define the dummy links for the three single-carrier networks as follows. A gas dummy link does not have a link equation, it only has a gas flow q. A power dummy link has a complex current I for which it holds that Iji= −Iij.
It also has active and reactive powers at the start and end of the line, for which it holds that Pji = −Pij := P and Qji = −Qij := Q. Note that P and Q are independent of I and V (also
because V is undefined for a heterogeneous node). A heat dummy link has a water flow m and supply and return temperature at the start and end of the link. Seeing it as a lossless line, it follows that Tstart
k = T end
k for both supply and return. The coupling variables, denoted by [·]
c, and their
associated location in the graph are shown in Figure 3.
A heterogeneous coupling node has one or more node equations associated with it, that relate the different energy carriers to each other. We assume these coupling equations to be of the form
fc(qc, Pc, ϕc) = 0 (16)
with qc the vector of gas flows of all gas dummy links connected to the coupling node, Pc the
vector of active powers of all power dummy links connected to the coupling node, and ϕcthe vector
of all heat powers of all heat terminal links connected to the coupling node. The heat terminal link represents some physical component that injects or consumes heat (with respect to the heat network), which means that there is a heat power equation (12) associated with it. Combining the coupling equations (16) with the heat power equation (12) gives the system of equation for all coupling nodes: Fc(xc) = f c(qc, Pc, ϕc) fϕ(ϕc, mc, Ts, Tr, To,c) ! = 0 with xc = qc Pc Qc mc ϕc To,c (17)
Since the dummy links used for coupling are homogeneous links, the system of equations for the single-carrier parts of the total integrated network are slightly altered (compared with before cou-pling). That is, the incidence matrices Ag, Ah, and Agalso include the dummy links. And, for heat,
the dummy links are also included in the mixing rules. However, the parameters of the dummy links are included in xc, and not in xg, xe, or xh.
Combining the systems of equations for the gas (5), power (9), heat (15), and coupling (17), gives a non-linear system of equations for the total integrated multi-carrier energy network:
F (x) = Fg(xg, xc) Fe(xe, xc) Fh xh, xc Fc xc, xh = 0 with x = xg xh xe xc (18)
Due to the heterogeneous coupling node and the dummy links, this system of equations gets a specific structure. First, the equations belonging to the single-carrier parts do not directly depend on variables of the other single-carrier networks. All dependencies are incorporated through the coupling. Furthermore, Fg and Fe are linear in the coupling variables xc (while Fh is non-linear
in xc because of the mixing rules (13)). Finally, the coupling system Fc depends on the heat
pa-rameters xh, because the heat power equation of the terminal link depends on Ts (for a sink) or
Tr(for a source) of the heat node connected to the coupling node (node ihin Figure 3).
3.4.2 Node types
A coupling will usually introduce more variables than equations, such that (additional) boundary conditions must be introduced. One option is to impose these on the coupling part of the system, see for instance Shabanpour-Haghighi and Seifi [2015] who assume all heat powers to be known. However, assuming the coupling energy flows known will essentially decouple the system. For cou-pling components that produce heat, it is possible to assume To,c known, as is done for any heat
source or sink, without decoupling the systems. Another option is to impose the additional bound-ary conditions in the single-carrier parts of the network. This will introduce new node types. An overview of the new and standard node types is shown in Table 2.
To analyze the effect that coupling has on the node types, we assume that the systems of equations of the single-carrier networks are formulated such that (i.e., the boundary condition are chosen such that) they are solvable. Because a dummy link only represent an (energy) flow going into or out of a node, it could be seen as a terminal link from the perspective of the single-carrier networks. However, we will adopt the perspective of the total integrated multi-carrier network, in which the dummy links are just links. Node types are based on the boundary conditions, that is, based on which nodal variables are assumed known and which are assumed unknown. Therefore, having a dummy link connected to a homogeneous node does not change the node type. Adding additional boundary conditions in the single-carrier networks will.
A gas terminal node has only two nodal variables, pressure p and injected flow q. Assuming p, q, or both known at a node where it was originally unknown changes the node type. For instance, assuming the pressure known at a node that was originally a load node turns it into a reference load node. However, no new node types are found, since all possible node types are already used
in a single-carrier gas network.
For a power grid, the technically feasible node types are a load node, a generator node, or a reference node. However, at power nodes that are connected to a coupling node, the injected active or reactive power can be assumed known, since the unknown coupling power can serve as slack for that node. So a reference node can become a P V δ, QV δ, or P QV δ node, and a generator node can become a P QV node. Similar nodes were introduced by Martinez-Mares and Fuerte-Esquivel [2012], who consider a generator P V node with variable active power, and a generator P Q node with variable active power. Note that, unlike in the gas network, the additional boundary condi-tions cannot be imposed ‘far’ from the coupling in the power grid.
A heat terminal node has six nodal variables, meaning that many combinations of known and unknown variables exist, some of which are not technically feasible. We will only look at the effect of the coupling on the four most commonly used node types. The resulting node types will therefore not be the complete set of possibilities. Due to technical reasons it is uncommon to assume the injected mass flow to be known. Furthermore, it is unlikely that both Ts and Tr are known for
one node. Furthermore, Tsfor a source node, or Trfor a load node, can only be assumed known if
To,cis unknown. Then, To,c will serve as ‘slack’ for the mixing rule. The original node types are
then affected as follows. At a source node, assuming the head known leads to a source reference node, and assuming the supply temperature known leads to a source temperature node. Assuming the head known for a load node leads to a load reference node and for a junction to a reference node. Usually To and ϕ are assumed known for every heat component (i.e. every terminal link), ex-cept for the source reference slack node. Assuming To,cknown for the terminal heat link connected to the coupling node does not decouple the system of equations, whereas assuming ϕcknown would. Therefore, we will consider this a possible boundary condition on a coupling node.
Not all combinations of node types (that is, not all combinations of boundary conditions) lead to a solvable system. Even for a small single-carrier network, some combinations of node types will lead to a system with no solutions or a system with multiple solution. Consider a gas network consisting of a source connected with a pipe to a sink, such as shown in Figure 2a. If both node i and node j are load nodes, the system will have no solution if the boundary conditions are imposed such that qi 6= qj, and infinite solutions if qi = qj, because only the pressure drop can be
deter-mined. Assuming the node types in single-carrier networks were such that the load flow problem was (uniquely) solvable, imposing the additional boundary conditions after coupling could lead to an unsolvable system. For instance, one single-carrier network could become overdetermined while another becomes underdetermined. However, a necessary condition for a solvable system, with respect to the number of each node type, can be derived by requiring that the number of equations equals the number of variables.
3.4.3 Newton-Raphson method
Newton-Raphson iteration can be used to solve a non-linear system of equations. The iteration scheme for multiple dimensions is given by:
x(k+1)= xk− J xk−1
Table 2: Node types for an integrated energy network, specified per energy carrier. New node types are shown in bold. The last column shows the notation for the set of node types, which are subsets of the set of nodes in the graph.
Network Node type Specified parameters Unknown parameters Number of nodes
Gas reference node p q |VRg|
load node q p |VLg|
slack node p, q |VSg|
reference load node p, q |VRLg |
Electricity slack bus |V |, δ P, Q |Ve
S|
generator bus (PV-node) P, |V | Q, δ |Ve
G|
load bus (PQ-node) P, Q |V |, δ |Ve
L| PVδ-node P, |V |, δ Q |Ve P V δ| QVδ-node Q, |V |, δ P |Ve QV δ| PQVδ-node P, Q, |V |, δ |Ve P QV δ| PQV-node P, Q, |V | δ |Ve P QV|
Heat source reference slack node Ts, h Tr, To, ϕ, ˙m |Vh
RS| source node To, ϕ Tr, Ts, h, ˙m |Vh S| load node To, ϕ Tr, Ts, h, ˙m |Vh L| junction node m = 0 Tr, Ts, h |Vh J|
source reference node To, ϕ, h Tr, Ts, m |Vh
SR|
source temperature node To, Ts, ϕ Tr, h, m |Vh
ST|
load reference node To, ϕ, h Tr, Ts, m |Vh
LR|
reference node m = 0, h Tr, Ts |Vh
R|
with k the iteration and J (x) the Jacobian matrix. Due to the choice for the heterogeneous coupling node connected with homogeneous dummy links, the Jacobian matrix of the integrated system of equations (18) is of the form
J (x) = Jgg Jge Jgh Jgc Jeg Jee Jeh Jec Jhg Jhe Jhh Jhc Jcg Jce Jch Jcc = Jgg 0 0 Jgc 0 Jee 0 Jec 0 0 Jhh Jhc 0 0 Jch Jcc (20)
where the submatrices are defined as Jαβ:=
∂Fα
∂xβ with α, β ∈ {g, e, h, c} (21)
Since the coupling parameters, except possibly To,c, are assumed unknown, the coupling generally
introduces more variables than equations. The submatrix Jcc is then not square, having more
columns than rows. Subsequently, as the additional boundary conditions are imposed in the single-carrier parts, the submatrices Jgg, Jee, and Jhh will also not be square, having more rows than
columns. This means we will solve the system as a whole, since solving the system blockwise is not possible.
However, if the (output) energy flows of the coupling components are known, the system becomes decoupled. The required additional boundary conditions can be imposed on the coupling variables by assuming (some of) them known. Imposing these boundary condition such that a (unique) so-lution exists means that Jccis square. Furthermore, if all the heat powers are assumed known, the
heat equation can be replaced with an equation ϕc− (ϕc)known
= 0, such that Jch = 0. Then the
coupling part can be solved for xg, independent of the single-carrier parts. The coupling variables can then be substituted into the boundary conditions of the single-carrier networks, which can then be solved independent of each other.
4
Case studies
To validate and illustrate the proposed model framework, we consider a small integrated energy sys-tem, based on a case study by Shabanpour-Haghighi and Seifi [2015]. Ayele et al. [2018] considered an adapted version of this case study, using an extended energy hub approach. For comparison, we consider two different ways of coupling the single-carrier networks, one similar to Shabanpour-Haghighi and Seifi [2015] and the other similar to Ayele et al. [2018]. Figure 4 shows the graph representation for both cases. The single-carrier networks consist of three nodes each. Gas node 0g is assumed to be connected to the national gas grid, such that this node will provide the energy for the rest of the integrated energy system. Coupling occurs at node 0 and node 2 of each network, whereas node 1 of each network is a sink. There is a compressor directly after node 1g on the link
to node 2gin the gas network.
The choice of node types in the integrated network determines if the system is solvable. For both case studies we derive a condition such that the system has equal number of equations and variables. We consider a set of node types, per case, for which the system is (uniquely) solvable. For both cases we use the same models for the single-carrier parts. In the gas network, links (0g, 1g) and (0g, 2g) represent a basic pipeline with
qk = Cksign p2i − p2j q p2 i − p2j (22)
where qkis the gas flow in standard m3/s and piis the nodal pressure in Pa. Ckis the pipe constant
of pipe k given by [Osiadacz, 1987]:
Ck = π Tn pn r Rs,air 16 s D5 k fkLkZT S (23) with Tn standard temperature in K, pn standard pressure, Rs,air the specific gas constant of air,
Dk the inner diameter of the pipe in m, Lk the length of the pipe in m, Z the gas compressibility
factor, S the relative density of gas to air, and fk the Darcy friction factor of the pipeline. The
friction factor is determined from the implicit Colebrook equation for the turbulent regime: 1 fk = −2 log10 k 3.7Dk + 2.51 Re√fk (24)
gas network electricity network heat network 2g 0 1e 0 1 2 2 1 1 0 2 3 4 4 1g 2e 2h 1h 0h 0e 0g 1c 2c 0c 3 3 4 5 (a) Case 1 gas network heat network electricity network 2g 0 1e 0 1 2 2 1 1 0 2 3 3 3 4 4 4 1g 2e 2h 1h 0h 0e 1c 0c 0g (b) Case 2
Figure 4: Network topology of the two case studies; (a) case 1 for comparison with Shabanpour-Haghighi and Seifi [2015]; (b) case 2 for comparison with Ayele et al. [2018]. In case 1 (a) node 0g is connected to a gas-fired generator, represented by node 0c, and to a gas-boiler, represented by node 1c. In case 2 (b), these two components are modeled by one energy hub, represented by node 0c. Similarly, in case 1 (a) node 2c represents a CHP, whereas in case 2 (b) node 119 c represents an
with k the pipe roughness in m and Re the Reynolds number given by:
Re = 4q
πDν (25)
with ν the kinematic viscosity in m2/s. Substituting T
n = 273.15K, pn = 1.01325 · 105Pa, and rescaling gives Ck= 1.2913 · 10−2 s D5 k fkLkZT S (26) with Dk in mm and Lk in km. Taking p in bar in equation (22) then gives qk is in standard m3/h.
Link (1g, 2g) represents a basic pipeline with a compressor connected at the beginning of the pipe,
such that qk = Cksign p2i − p 2 j r (Hpi) 2 − p2 j (27)
with Ck given by equation (26), and H = ppout
in the compression ratio.
We use the nodal formulation, and only use the equations for nodes where the injected gas flow is known, that is only load nodes and reference load nodes. The system of equations for the gas part is then given by Fg(xg) = Ag0q − qinj0= 0 with xg = p (28) with Ag0 the (|Vg L| + |V g RL|) × |E
g| reduced incidence matrix, qinj0 the reduced vector of injected
flows, and qk is given by the link equation (22) for non-dummy links. This reduced system consist
of |VLg| + |VRLg | equations and |VLg| + |VSg| variables.
In the power grid, all (non-dummy) links are represented by a short line model [Schavemaker and Van der Sluis, 2008], such that:
Iij= yij(Vi− Vj) (29)
where yij = gij+ ιbij is the admittance of the line, with yij = yji. Substituting the link equation
(29) in the complex power equation leads to
Pij =
gij|Vi|2− |Vi||Vj| (gijcos δij+ bijsin δij) for non-dummy links
Pk for dummy links
0 otherwise Qij =
−bij|Vi|2− |Vi||Vj| (gijsin δij− bijcos δij) for non-dummy links
Qk for dummy links
0 otherwise
(30)
We use the complex power formulation, and only take the equations for the nodes where the injected active or reactive power is known. The system of equations for the electrical part is then given by:
Fe(xe) = P j, j6=i Pij− Pinj0 P j, j6=i Qij− Qinj0 = 0 with x e= δ |V | (31)
Here, P
j, j6=i
Pij and P j, j6=i
Qij are the reduced vectors with calculated injected active and reactive
power respectively, Pinj0 and Qinj0 are the reduced vectors of known injected active and reactive
energy flows. This reduced system consists of |Ve
G| + 2|VLe| + |VP V δe | + |VQV δe | + |VP QV δe | + |VP QVe |
equations and |Ve
G| + 2|VLe| + |VP QVe | variables.
In the heat network, all (non-dummy) links represent a pipeline with a basic head loss equation [Ayele et al., 2018]:
hi− hj− Kk|mk|mk = 0 (32)
with Kk the pipe constant given by:
Kk =
8Lkfk
π2gρ2D5 (33)
with Dkthe inner diameter of the pipe in m, Lkthe length of the pipe in m, ρ the density of the water
in kg/m3, and f
k the Darcy friction factor of the pipeline. The friction factor is determined from
the implicit Colebrook equation for the turbulent regime (24). For the thermal model we consider the temperatures with respect to the ambient temperature. That is, we consider T0 := T − Ta,
with Ta the ambient temperature in◦C. For notational simplicity, we drop the prime and simply
denote this temperature as T . The pipeline has an exponential temperature drop for both supply and return:
Tkend− ψ (mk) Tkstart= 0 with ψ (mk) := exp
−λkLk
Cp|mk|
(34) where λkis the heat transfer coefficient of the pipe in W/(m K), and Cpis the specific heat of water
in J/(kg K). For the terminal links the heat power equation is given by
ϕi,l= Cpmi,l Tis− To i,l if node i is a sink Cpmi,l Ti,lo − Tr i if node i is a source (35)
The mass flow of a terminal link l connected to node i can then be found by
mi,l= ϕi,l Cp To i,l− Tir if vi∈ V h S ∪ VSRh ∪ VSTh ϕi,l Cp Ts i − T o i,l if vi∈ V h L∪ VLRh 0 if vi∈ VJh∪ V h R X Eh ms if vi∈ VRSh (36)
The equation for the source reference slack node is based on the assumption that there is only one component connected to it, and that ϕi,lis unknown. Substituting the temperature drop (34) and
the components mass flow (36) in the mixing rule (13) leads to FTs = 0 for the supply line, with FiTs = X Eh msout ! Tis−X j |mij|ψijTjs− X l −mi,lTi,lo if vi∈ VSh∪ V h SR∪ V h ST X Eh msout+ X l mi,l ! Tis− X j |mij|ψijTjs if vi∈ VLh∪ V h LR X Eh msout ! Tis−X j |mij|ψijTjs if vi∈ VJh∪ VRh X Eh msout ! Tis−X j |mij|ψijTjs− (−miTio) if vi∈ VRSh (37)
Similarly, the mixing rule (13) for the return temperature is then given by FTr= 0, with
FiTr = X Eh msin+X l −mi,l ! Tir−X j |mij|ψijTjr if vi∈ VSh∪ VSRh ∪ VSTh X Eh msin ! Tir−X j |mij|ψijTjr− X l mi,lTi,lo if vi∈ VLh∪ V h LR X Eh msin ! Tir− X j |mij|ψijTjr if vi∈ VJh∪ V h R X Eh msin− mi ! Tir−X j |mij|ψijTjr if vi∈ VRSh (38)
While for the gas network the nodal formulation is commonly used, this is not done for the heat network. The nodal formulation would make the Jacobian matrix complicated due to the depen-dence of edge mass flow on heads and, in turn, the dependepen-dence of the mixing rules on these edge flows. One of the commonly used formulations is the loop formulation which is for instance used by Liu et al. [2016] and Shabanpour-Haghighi and Seifi [2015]. However, this requires loops to be found in the network. Furthermore, it is common to substitute conservation of mass (10) in the heat equation (12), resulting in a ‘conservation of heat’. However, this formulation is inconvenient when defining dummy links, and the interpretation for a node with multiple components is unclear. Therefore, we formulate the system based on Ayele et al. [2018] for the hydraulic part and on Abeysekera [2016] for the thermal part. Essentially, we do not reduce the size of the hydraulic part and only apply algebraic substitutions in the thermal part. The conservation of flow and the supply mixing rule are not taking into Fhfor source reference slack nodes. Similarly, the head and supply
temperature are not part of xh for source reference slack nodes. The head is also not part of xhfor
source reference, load reference or reference nodes. And the supply temperature is not part of xh
for source temperature nodes. The system of equations for the heat part is then given by
Fh xh = Ah0m − minj0 AhT h − Diag (K) Diag (|m|) m FTs FTr = 0 with xh= m h0 Ts0 Tr (39)
Table 3: Values of parameters for both case studies, per network. Here, | · |bdenotes the base value,
and p.u. stands for per unit.
Network Network parameters Carrier parameters Link Link parameters gas |Eg|b[MW] Tn[◦C] pn105Pa GHV10−2MW/m3 Z [−] S [−] ν10−6m/s L [km] D [mm] [mm] H [−] 100 273.15 1.01325 1.19034 0.8 0.6106 0.288 0-1 30 150 0.05 0-2 30 150 0.05 1-2 30 150 0.05 1.3
power |S|b[MW] x [p.u.] r [p.u.]
100 0-1 0.5 0.05 0-2 0.5 0.05 1-2 0.5 0.05 heat |ϕ|b[MW] Tref high[◦C] T a [◦C] Cp10−3MJ K/kg ρkg/m3 ν10−6m/s L [km] D [mm] [mm] λ [W K/m] 100 130 10 4.182 960 0.294 0-1 30 150 1.23 0.2 0-2 30 150 1.23 0.2 1-2 30 150 1.23 0.2 with minji =P l
mi,l, Ah0 the reduced incidence matrix, minj0 the reduced vector of injected flows,
and Ts0 and Tr0 the reduced vector of supply and return temperatures respectively. Diag (v) is used to denote a diagonal matrix with vector v on the diagonal. This reduced system con-sists of 3|Vh S| + 3|V h L| + 3|V h J| + |V h RS| + 3|V h SR| + 3|V h ST| + 3|V h LR| + 3|V h R| + |E h N D| equations and 3|Vh S| + 3|V h L| + 3|V h J| + |V h RS| + 2|V h SR| + 2|V h ST| + 2|V h LR| + 2|V h R| + |E h N D|variables. Here, E h N Ddenotes
the set of non-dummy heat links in the network.
To solve the total system of equations (of the form (18)), we adopt a per unit formulation for the single-carrier networks and the coupling part. This normalization is commonly used in power grids. Appendix A gives details about the formulation for all three networks. Table 3 gives the values for the network parameters used in both cases.
4.1
Case 1
Case 1 is used for comparison with the case study of Shabanpour-Haghighi and Seifi [2015]. Node 0c represents a gas-fired generator, node 1c a gas boiler, and node 2c a combined heat and power
plant (CHP). For the coupling components we use the same models as Shabanpour-Haghighi and Seifi [2015]: f0c(q3, P3) = aP32+ bP3+ c + d sin e Pmin− P3 − GHVq3= 0 f1c(q4, ϕ1c) = ϕ1c− ηGBGHVq4= 0 f2c(q5, P4, ϕ2c) = P4+ ϕ2c− ηCHPGHVq5= 0 (40)
with a, b, c, d, and e parameters of the gas-fired generator and Pminthe minimum produced power, ηGB the efficiency of the gas boiler, and ηCHP of the CHP. We take a = 0.002 931, b = 1.1724,
c = 43.965, d = 4.3965, e = 0.5, ηGB= 0.88, and ηCHP = 0.88.
equations for the outflow temperature gives the system of equations for the coupling part: Fc(xc) = P3− ηGGGHVq3 ϕ1c− ηGBGHVq4= 0 P4+ ϕ2c− ηCHPGHVq5 ϕ1c− Cpm3 T1oc− T0rh ϕ2c− Cpm4 T2oc− T2rh To 1c− (T1oc) known T2oc− (T2oc) known = 0 with xc= q3 q4 q5 P3 P4 Q3 Q4 m3 m4 ϕ1c ϕ2c To 1c T2oc (41)
A system like this, where every coupling node has only one coupling equation fc and a maximum
of one dummy link per energy carrier, consists of |Vc
ge| + 2|Vgh,Sc | + 3|V c gh,T| + 2|V c eh,S| + 3|V c eh,T| + 2|Vc geh,S|+3|V c
geh,T| equations and 3|V c
ge|+4|Vghc |+5|V c eh|+6|V
c
geh| variables. Here, V c αβand V
c αβγ
de-note the set of coupling nodes which are linked to networks with energy carriers α, β, γ ∈ {g, e, h}. The subscript S denotes a standard coupling node, with unknown outflow temperature To,c, while
the subscript T indicates a temperature coupling node with known outflow temperature. For the latter, an equation To,c− (To,c)known
= 0 is added to the system of equations.
Combining the systems of equations (28), (31), (39), and (41), results in a total system of equations of the form (18). This system has the same number of equations as variables, that is |F (x) | = |x| if |Vg RL| + |V e P V δ| + |V e QV δ| + 2|V e P QV δ| + |V e P QV| + |V h SR| + |V h ST| + |V h LR| + |V h R| =|V g S| + 2|V c ge|2|V c gh,S| + |V c gh,T| + 3|V c eh,S| + 2|Vc
eh,T| + 4|Vgeh,Sc | + 3|Vgeh,Tc |
(42) Assuming the node types in the single-carrier were such that they were solvable, these couplings require 8 additional boundary conditions to be added. We will consider the set of nodes types shown in Table 4a. The power grid and the heat network do not have any external inflow; all the energy is provided by the gas network. This means that that node 0h is a junction. Node 0e has
a heat pump connected to it in the case study by Shabanpour-Haghighi and Seifi [2015], such that the total injected active power for node 0e is equal to the power supplied to that heat pump. Both
these nodes will be taken as reference nodes, which means that 0eis a P QV δ-node and node 0his a
reference node. We assume the outflow temperature of the gas-boiler and the CHP to be known, as is done for any heat source or sink. Then, three more additional boundary conditions are needed; one is imposed on node 1g, one on node 2e, and the last on node 2h.
Table 4: Node types for both cases. The x and F column show the amount of variables and functions the node contributes to the system. For the total system of equations, all non-dummy heat links contribute one entry to x and one to F .
(a) Case 1. Node Type x F 0g ref. 0 0 1g load 1 1 2g ref. load 0 1 0e PQVδ 0 2 1e load 2 2 2e PQV 1 2 0h ref. 2 3 1h load 3 3 2h load ref. 2 3 0c ge 3 1 1c gh temp. 4 3 2c geh temp. 6 3 Total 27 27 (b) Case 2. Node Type x F 0g ref. 0 0 1g load 1 1 2g load 1 1 0e PQVδ 0 2 1e load 2 2 2e PQV 1 2 0h ref. 2 3 1h load 3 3 2h load 3 3 0c geh temp. 6 4 1c geh temp. 6 4 Total 28 28
gas network heat network electricity network GB GG
q
3ϕ
3P
3 1− ν0 ν0(a) Energy hub 0
CHP
P
4
ϕ
4q
4(b) Energy hub 1
Figure 5: Schematic representation of the energy hubs. (a) shows energy hub 0 consisting of a gas-fired generator (GG) and a gas boiler (GB). ν0denotes the fraction of q3provided to the generator.
(b) shows energy hub 1 consisting of a CHP. Here, q is gas flow, P is active power, and ϕ is heat power. Energy hubs 0 and 1 are represented by nodes 0c and qc in Figure 4b respectively.
4.2
Case 2
Case 2 is used for comparison with the case study of Ayele et al. [2018]. Both nodes 0c and 1c represent an energy hub, called energy hub 0 and energy hub 1 respectively. Energy hub 0 consists of a gas-fired generator and a gas-boiler, energy hub 1 consists only of a CHP. Figure 5 shows a schematic representation of both energy hubs.
For both energy hubs we use a coupling matrix C as suggested by Geidl and Andersson [2007]. The entry cαβ of C gives the conversion factor from energy carrier α to carrier β. Based on the
schematic representation of the energy hubs we find: GHVq3 P3 ϕ0c out = 0 0 0 ηGGν0 0 0 ηGB(1 − ν0) 0 0 GHVq3 P3 ϕ0c in := C0 GHVq3 P3 ϕ3 in GHVq4 P4 ϕ1c out = 0 0 0 ηCHPν1 0 0 ηCHP(1 − ν1) 0 0 GHVq4 P4 ϕ1c in := C1 GHVq4 P4 ϕ1c in (43)
The coupling matrix itself allows for bidirectional flow. If for instance ceg, chg6= 0, these hubs could
also convert power and heat back to gas. For this case, unidirectional flow is assumed, since the physical components can only convert gas to power and heat. The coupling equations are then given by: P3 ϕ0c = ηGGν0 ηGB(1 − ν0) GHVq3 P4 ϕ1c = ηCHPν1 ηCHP(1 − ν1) GHVq4 (44)
We take ηGG= 0.45, ηGB= 0.88, ηCHP = 0.88, ν0≈ 0.775 05, and ν1≈ 0.266 34
The system of equations for the coupling part is then given by
Fc(xc) = P3− c0geGHVq3 ϕ0c− c0ghGHVq3 P4− c1geGHVq4 ϕ1c− c1ghGHVq4 ϕ0c− Cpm3 T0oc− T0rh ϕ1c− Cpm4 T1oc− T1rh To 0c− (T0oc) known To 1c− (T1oc) known = 0 with xc = q3 q4 P3 P4 Q3 Q4 m3 m4 ϕ0c ϕ0c T0oc To 1c (45)
A system like this, where every coupling node has an energy hub matrix equation as coupling equa-tion, and a maximum of one dummy link per energy carrier, consists of |Eout
ge ||Vgec |+ |Eout gh,S| + 1 |Vc gh,S|+ |Eout gh,T| + 2 |Vc gh,T|+ |Eout eh,S| + 1 |Vc eh,S|+ |Eout eh,T| + 2 |Vc eh,T|+ |Eout geh,S| + 1 |Vc geh,S|+ |Eout geh,T| + 2 |Vc geh,T| equations and 3|Vc ge| + 4|Vghc | + 5|V c eh| + 6|V c geh| variables.
Combining the systems of equations (28), (31), (39), and (45), results in a total system of equations of the form (18). Assuming the node types in the single-carrier networks were such that, before coupling, those networks were solvable, these couplings require 6 additional boundary conditions. We will consider the set of node types shown in Table 4b. Similar to case 1, 0eis a P QV δ-node, 1e
a P QV -node, and node 0h is a reference node. We assume the outflow temperature of both energy
hubs to be known, as is done for any heat source or sink.
Compared with case 1, 6 additional boundary conditions are needed instead of 8. This is because the energy hub concept specifies the ratio between the input powers for hub 0, and the output powers for hub 1. Instead of imposing 2 additional boundary conditions in case 1, we could also add 2 equations in accordance with the energy hubs. The first would give the ratio between q3and
q4, the second between P4and ϕ4. Conversely, if ν0and ν1are unknown, or functions of the in- or
output power, the ratio between the powers in the energy hub becomes unknown. Then, a total of 8 additional boundary conditions would be required in case 2.
4.3
Results
We use a basic Newton-Raphson method to solve the system of equations for both cases. With a tolerance of 10−6p.u. the method converged in 4 iterations for case 1 and 5 iterations for case 2.
Table 5: Results for the gas network, (a) for case 1 with node types set 1, and (b) for Shabanpour-Haghighi and Seifi [2015].
Node p [bar] qinj 103standard m3/h
Link q103standard m3/h
(a) (b) (a) (b) (a) (b)
0 50 50.0000 -46.715 -47.1083 0-1 18.2327 29.8308
1 29.1021 40.8160 10.8649 10 0-2 16.4078 4.8098
2 34.0770 49.7827 20 20 1-2 7.3678 18.9658
Table 6: Results for the electrical network, (a) for case 1 with node types set 1, and (b) for Shabanpour-Haghighi and Seifi [2015].
Node |V | [p.u.] δ [◦] Sinj[MW]
(a) (b) (a) (b) (a) (b)
0 1.06 1.0600 0 0 0.1451+0ι 0+0ι
1 0.9801 0.9800 -6.9888 -7.0219 30+15ι 30.0000+15.0000ι
2 1 1.000 -6.048 -6.1156 30.1360+15ι 30.1360 +15.0000ι
Link Sij[MW] Sji[MW] Sijloss[MW]
(a) (b) (a) (b) (a) (b)
0-1 26.8615+15.8007ι 26.9806+15.8106ι -26.4293-11.4789ι -26.5454-11.4588ι 0.4322+4.3218ι 0.4352+4.3518ι 0-2 23.4916+11.5508ι 23.7405+11.5524ι -23.1867-8.5013ι -23.4303-8.4505ι 0.3049+3.0495ι 0.3102+3.1019ι
1-2 -3.5707-3.5211ι -3.4546-3.5412ι 3.5838+3.6520ι 3.4673+3.6686ι 0.0131+0.1309ι 0.0127+0.1274ι
Total 0.7502+7.5022ι 0.7581+7.5811ι
Table 7: Results for hydraulic part of the heat network, (a) for case 1 with node types set 1, and (b) for Shabanpour-Haghighi and Seifi [2015].
Node h [m] minj [kg/s] Link m [kg/s]
(a) (b) (a) (b) (a) (b)
0 5517 - 0 - 0-1 64.6877 64.6943
1 224.9319 - 121.2256 - 0-2 31.4083 31.4533
