- Email: [email protected]
Matrix modeling of energy hub with variable energy efficiencies
Electrical Power and Energy Systems 119 (2020) 105876
Contents lists available at ScienceDirect
Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes
Matrix modeling of energy hub with variable energy efficiencies a
a,⁎
b
c
c
a,⁎
Wujing Huang , Ning Zhang , Yi Wang , Tomislav Capuder , Igor Kuzle , Chongqing Kang a b c
T
Department of Electrical Engineering, Tsinghua University, China Department of Information Technology & Electrical Engineering, ETH Zurich, Switzerland Faculty of Electrical Engineering and Computing, University of Zagreb, Croatia
ARTICLE INFO
ABSTRACT
Keywords: Multi-energy systems Energy hub Matrix modeling Piecewise linearization Variable efficiency Part-load performance Operation optimization
The modeling of multi-energy systems (MES) is the basic task of analyzing energy systems integration. The variable energy efficiencies of the energy conversion and storage components in MES introduce nonlinearity to the model and thus complicate the analysis and optimization of MES. In this paper, a standardized matrix modeling approach is proposed to automatically model MES with variable energy efficiencies based on the energy hub (EH) modeling framework. Piecewise linearization is used to approximate the variable energy efficiencies; as a result, a component with variable efficiency is equivalent to several parallel components with constant efficiencies. Splitters and concentrators are proposed as standardized components to facilitate the split and merge of energy flows imposed by piecewise linearization. The nonlinear energy conversion and storage relationship in EH can thus be further modeled under a linear modeling framework using matrices. Such matrix modeling approach makes the modeling of an arbitrary EH with nonlinear energy components highly automated by computers. The proposed modeling approach can further facilitate the operation and planning optimization of EH with variable efficiencies. Case studies are executed in MATLAB and presented to show how the nonlinear approximation accuracy and calculation efficiency can be balanced using the proposed model in the optimal operation of EH. For the optimal operation of a five-component EH with energy storage, the proposed approach reduces the approximation error from 13.7% to 0.1% with only 0.2 s computation time increase compared to the exiting constant efficiency model.
1. Introduction The interactions between electricity, gas and heat/cooling systems enable higher energy efficiency and utilization of intermittent renewable energy sources [1–3]. The modeling of multi-energy systems (MES) has gained increasing importance in recent years as the foundation of the operation and planning of MES, which can unlock the flexibility of shifting across multiple energy sectors and result in reduced costs and lower emissions as compared to separate operation and planning of energy systems. The coupling relationship among different energy sectors in MES and the interaction characteristics between MES and external networks/loads are the key issues that the modeling approach needs to address. The energy hub (EH) concept was introduced to model MES as a unit in which multiple energy carriers can be converted, conditioned and stored. EHs import energy at their input ports, which are connected to external energy infrastructures, e.g., power grid and gas network. They provide required energy services, such as electricity, heat and cooling, at their output ports [4]. The compactness and effectiveness of EH-based model facilitate the calculation of energy
⁎
flows in MES of different levels, such as trigeneration systems [5], residential consumers [6], commercial buildings [7], industrial parks [8], national energy systems [9] and interconnected countries [10]. In addition, EH can be integrated with smart technologies especially smart meters, becoming smart energy hub (SEH), to handle a large volume of information in MES [11]. While a substantial number of papers have discussed the optimal planning [12,13], operation [14,15] and carbon emission [16,17] of MES based on the EH model, very few papers have studied the standardized and automatic modeling of MES. Based on the concept of EH, a nonlinear framework for modeling of energy systems with conversion and storage of multiple energy carriers is proposed in [18]. In this modeling framework, the converter and storage coupling matrices are manually formulated according to the given EH layout. Reference [5] proposes an approach to automatically generate coupling matrices for small-scale trigeneration systems using dispatch factors [13]. The formulation involves products of dispatch factors and energy flows, leading to nonlinear rather than linear optimization problems. As a result, the approach is more applicable to small-scale and simple-
Corresponding authors at: Room 3-206, West Main Building, Tsinghua University, Haidian, Beijing, 100084, China. E-mail addresses: [email protected] (N. Zhang), [email protected] (C. Kang).
https://doi.org/10.1016/j.ijepes.2020.105876 Received 16 July 2019; Received in revised form 20 December 2019; Accepted 19 January 2020 0142-0615/ © 2020 Elsevier Ltd. All rights reserved.
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
Nomenclature AB CCHP CERG CHP EH HP
HS MES MILP MIMO NP PWL SIMO SISO
Gas-fired auxiliary boiler Combined cooling, heat and power Compression electric refrigerator group Combined heat and power Energy hub Electric heat pump
structure MES. Reference [19] introduces a linear and automatic modeling tool, Hubert, which leverages a concise ASCII description format for EH with an energy “input, input storage, converter, output storage, output” structure. Reference [20] improves Hubert by adding an “export” branch between “converter” block and “output storage” block to facilitate selling PV-generated energy to the external grid. Reference [21] proposes the idea of a linear coupling relationship of an EH by introducing the concept of augmented variables. The standardized procedure for modeling the dispatch-factor-free coupling matrix is then designed to facilitate the computerized calculation for arbitrary EH configurations. Reference [22] proposes an automatic and standardized matrix modeling method for EH with arbitrary configuration. The method uses graph theory to cast the topology and the characteristics of the energy converters into a matrix form and uses linear equations for the dispatch of energy flows without using dispatch factors. However, all the energy efficiencies are regarded as constants in the above research. In reality, energy efficiencies of energy conversion and storage components vary with their operating conditions. The products of variable efficiencies and energy flows bring nonlinearity to the components and systems. Energy conversion efficiencies of converters can reduce significantly when they do not work in rated condition; therefore, the constant efficiency approximation model of MES with high nonlinearity is very likely to be inaccurate and impractical. There are already some methods taking into consideration the variable efficiencies of energy conversion and storage components. In [23], nonlinear energy converters are directly modeled with highly nonlinear part-load efficiency curves. The resulting model is a nonlinear programming (NP) problem which gives no guarantee that the global optimum can be achieved and few numerically robust solvers for NP problems exist. To avoid NP problems, several linearization approaches are used to model nonlinear converters. The first one is step-wise approximation based on binary variables which guarantees that only one operating point is selected in each time step (e.g., [24,25]). This formulation carries the disadvantage that it only allows single tabled values to be selected, with no interpolation between them. The second one is piecewise linear (PWL) approximation using special-ordered-set (SOS) variables (e.g., [26]). A special ordered set is a set of consecutive variables in which not more than one (named SOS1) or two (named SOS2) adjacent members may be non-zero in a feasible solution. By using SOS2 variables as weight variables, any operating point can be written as a linear combination of two pre-defined operating points. Although these methods are able to address variable efficiencies, the formulation can hardly be automated by computers. Reference [27] improves Hubert by enabling the modeling of nonlinear converters using PWL approximation. The operating points are defined as the sum of a set of consecutive variables related to different PWL segments, which has been widely used to approximate quadratic fuel cost and transmission loss in power system operation and planning. However, this model can still only deal with the proposed structure of MES where the location of energy storage components is fixed and only one energy conversion block is considered in an EH. In addition, these methods cannot handle the nonlinear energy converters with multiple outputs and complex operation characteristics, e.g., combined heat and power (CHP) in extraction condensing operation mode in which the ratio between electricity and heat production is adjustable within a certain
Heat storage Multi-energy systems Mixed integer linear programming Multi-input multi-output Nonlinear programming Piecewise linear Single-input multi-output Single-input single-output
range. In this paper, we propose a standardized matrix modeling approach for EH with variable energy efficiencies. Using piecewise linearization, the nonlinear energy conversion or charging/discharging process in an energy component is divided into several linear processes with constant efficiencies. Splitters and concentrators are proposed as standardized components to facilitate the split and merge of energy flows imposed by piecewise linearization; as a result, a nonlinear energy component is converted into a “splitter, linear energy component, concentrator” three-layer structure. The nonlinear energy conversion and storage relationship in EH can thus be further modeled under a linear and matrixbased modeling framework. Such approach facilitates an automatic modeling of an arbitrary EH with nonlinear energy components. Compared to the constant efficiency model and the nonlinear model, the proposed standardized matrix model is proven to significantly improve the accuracy and calculation efficiency, respectively, of the operation optimization of EH. The contributions of this paper are threefold: (1) Providing a standardized matrix modeling approach for EH with nonlinear energy conversion and storage components, which makes the modeling of an arbitrary EH highly automated by computers. Such technique fills the gap that current nonlinear/linearized modeling approaches can hardly automatically formulate an arbitrary EH with nonlinear energy components. (2) Proposing a generalized linearized model that is capable of accurately modeling any energy conversion or storage component with constant or variable efficiencies under one linear modeling framework. Especially, the generalized linearized model is able to handle the nonlinear energy converters with multiple outputs and complex operation characteristics, e.g., CHP in extraction condensing operation mode. (3) Providing analysis of how the nonlinear approximation accuracy and calculation efficiency can be balanced using the proposed model in the optimal operation of EH. The rest of this paper is organized as follows: Section 2 provides the basic backgrounds of the standardized matrix model of EH. Section 3 proposes the standardized matrix modeling approach to tackle the challenges of variable energy efficiencies of energy components. Section 4 presents the standardized data structure of the proposed modeling approach. Section 5 conducts a comprehensive case study and applies the proposed modeling approach to the operational optimization problem. Section 6 provides a discussion on the proposed modeling approach. Section 7 draws the conclusion. 2. Background: Standardized matrix model of EH with constant efficiencies First, a standardized matrix model of EH with constant energy efficiencies [22] is briefly introduced to facilitate understanding of the proposed matrix model of EH with variable efficiencies and the relationship between these two models. From the viewpoint of graph theory, each energy conversion or storage component in an EH can be treated as a node, and each energy 2
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
flow to or from a component can be treated as a branch, as shown in Fig. 1. The inputs and outputs of an EH are treated as special nodes. They are also related to one or more energy flows but not convert or store energy. Each node has several input and output ports through which it exchanges energy with other nodes.
conversion equation of the EH can be obtained:
ZV = 0 D. Comprehensive energy flow equations for EH
The m -dimensional energy input vector Vin and the n -dimensional energy output vector Vout are defined if an EH has m inputs and n outputs. The input incidence matrix X is an m × B matrix that relates the energy inputs to the branch energy flows. Similarly, the output incidence matrix Y is an n × B matrix that relates the energy outputs to the branch energy flows:
A. Formulating port-branch incidence matrices A port-branch incidence matrix A g of node g defines the connections between the ports of node g and all the branches in the EH. If node g has K g ports and the EH has B branches, the port-branch incidence matrix A g has dimensions K g × B and is defined as follows:
Ag (k , b) =
1, branch b is connected to input port k of node g - 1, branch b is connected to output port k of node g 0, otherwise
(1)
where k is the serial number of input and output ports of a node.
(7) (8)
Given the port-branch incidence matrix A g and the converter characteristic matrix Hg , the nodal energy conversion matrix Zg for node g can be calculated: (3)
Zg = Hg Ag
The energy conversion matrix Zg for node g defines the relationship between the energy flows in the branches connected to node g . This matrix has as many rows as there are energy conversion processes in node g and as many columns as there are branches in the EH. Suppose that an EH has N nodes. The system energy conversion matrix Z combines the nodal energy conversion matrix of all nodes in the EH: T
ZTN ]
Energy conversion efficiency is taken as an example here. (1) Type 1: SISO energy converters
The energy flows in all branches form a B × 1 vector V . According to the energy conversion characteristic of each node, the energy
v1
Port
(9)
A. Piecewise linearization of energy efficiency
(4)
Node
T
T T Vout 0T ] Y T Z T ]T V = [ Vin
When the energy efficiencies of energy conversion and storage components vary with their input or output powers, the constant efficiency model of EH is no longer accurate. Fig. 2 shows the output power curve of a nonlinear converter and its linear approximation which has a significant error at a certain load level. The approximations of variable efficiencies of production units, particularly CHP units, have a high impact on short-term schedules [28,29]. The products of variable efficiencies and energy flows also bring about nonlinearity to the optimization problem. In this paper, EH with variable efficiencies is modeled using a piecewise linearization method. Through piecewise linearization, a nonlinear energy conversion or storage component is divided into several equivalent components with constant efficiencies. The modeling of two types of energy converters, i.e., single-input single-output (SISO) and single-input multi-output (SIMO) energy converters, and the modeling of energy storage components are presented. The modeling of energy storage components is similar to that of SISO energy converters where the input port models charging and the output port discharging.
C. Formulating energy conversion matrices
Gas
1, branch j is connected to output port i 0, otherwise
3. Proposed standardized matrix modeling approach
is the energy conversion efficiency of energy conversion pro-
vin ,1
Y (i , j ) =
The comprehensive energy flow equations of the EH are:
(2)
Input
(6)
[ XT
s , input port k is related to energy conversion process s Hg (s, k ) = 1, output port k is related to energy conversion process s 0, otherwise
Z = [ Z1T Z2T
1, branch j is connected to input port i 0, otherwise
Vin = XV, Vout = YV
The converter characteristic matrix Hg of node g defines the energy conversion characteristics of the node. If node g has Sg energy conversion processes and K g ports, the converter characteristic matrix Hg has dimensions Sg × K g . If each energy conversion process in node g is related to only one input and one output, then Hg is defined as follows
s
X (i , j ) =
The input incidence and output incidence equations are:
B. Formulating converter characteristic matrices
where cess s .
(5)
Electricity
v2
Branch
v3
#1
v4
Heat
#2
v5
Cooling
Fig. 1. Definitions of basic elements of an EH in terms of graph theory. 3
Output
vout,1
Output
vout,2
Output
vout,3
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
Fig. 2. Nonlinear component variable efficiency linear approximation and piecewise linear approximation.
Suppose that node g has a single input vgin and a single output vgout and that its energy conversion efficiency g (vgin) is a function of the input power:
vgout =
in in g (vg )·vg
k
where v gin, k is the range of the k -th segment and
k= 1
vgin, k
vgin
=
k 1 i=1
v gin, i 0,
and vgin, k (k = 1, 2,
,
k 1 i=1
vgin
k i=1
v gin, i
v gin, k = v gin . Secondary
vgout ,k =
k 1 i=1
g,k
v gin, i
< vgin < v gin, i
k i=1
(14)
s ) are introduced and
vgout ,k
(15)
In the k -th (k = 1, 2, , s ) segment, the variable efficiency g (vgin) is replaced by a fixed efficiency g , k determined by the secant line of the power output curve in this segment:
variables vgin, k (k = 1, 2, s ) are introduced to represent the amounts of input power that fall in each segment:
v gin, k , vgin
,s
k=1
(11) s
v gin, i , k = 1, 2,
i=1
s
vgout =
[0, v gin,1], [v gin,1, v gin,1 + v gin,2], [v gin,1 + v gin,2, v gin,1 + v gin,2 + v gin,3], +v gin, s]
i=1
k
Similarly, secondary variables vgout , k (k = 1, 2, subject to
As shown in Fig. 2, the range of vgin (denoted by v gin ) can be divided into s segments:
+ v gin, s 1, v gin,1+
v gin, i ·
g
i=1
(10)
, [v gin,1+
k
v gout ,i =
in g , k ·vg , k
, k = 1, 2,
,s
in = v gout , k vg,k
(17)
Binary variables ug ,1, ug ,2 , continuity of vgin and vgout :
v gin, i
ug,1·v gin,1 (12)
s ) are subject to
(16)
vgin,1
, ug , s
1
are introduced to guarantee the
v gin,1
ug ,2· v gin,2
vgin,2
ug ,1· v gin,2
u g,3·v gin,3
vgin,3
ug ,2· v gin,3
s
vgin =
vgin, k k=1
ug , s 1·v gin, s
(13)
0
The range of vgout can be divided into the same amount of segments as that of vgin :
vgin, s
1
vgin, s
1
ug , s 2·v gin, s
1
ug , s 1·v gin, s
If vgin is in the k -th (1 < k < s ) segment, then: 4
(18)
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
ug ,1, ug ,2, , ug , k 1 = 1 ug , k , u g, k + 1, , u g, s 1 = 0
In this case, the input and output ports of the component should be split into multiple ports. Splitters and concentrators are thus introduced to facilitate the split and merge of branches as shown in Fig. 4. In this paper, the original branches in the EH are named primary branches, and the branches split from primary branches are named secondary branches. If a nonlinear SIMO converter has one input and b outputs, then one splitter, b concentrators and one linear mapping component should be introduced, as shown in Fig. 5. The energy flows in all secondary branches form a vector V (energy flows in secondary branches connected to the input ports of each node come first in V ), and the primary branch vector V should be expanded as:
(19)
Eq. (18) can be written in matrix form to make it easily automated by computers:
V in g Ug , a
V in g
V in g Ug , b
(20)
where in in V in g = diag (v g ,1, v g ,2,
, v gin, s ) T
in in V in g = [ vg,1 vg ,2
vgin, s ]
Ug , a = [ Ug 0 ]T , Ug , b = [1 Ug ]T ug , s 1 ] Ug = [ug,1 ug,2
(21)
C. Formulating port-branch incidence matrices
(2) Type 2: SIMO energy converters
First, consider the simplest case where the EH has only one nonlinear SISO energy converter, which is represented by node g , and its input and output ports are connected to only one branch respectively. Using the piecewise linearization method, node g is converted into an s -input s -output converter, and 2s secondary branches are introduced. A portbranch incidence matrix A g of node g defines the connections between the 2s ports and all the branches (including primary and secondary branches). If the EH has B primary branches, then the port-branch incidence matrix A g has dimensions 2s × (B + 2s ) and is defined as follows:
If the energy output of each output port is proportional to the energy input, e.g., CHP in backpressure operation mode, the energy conversion processes from the input port to each output port can be piecewise linearized separately in a similar way to SISO energy converters. A more complicated situation is when the proportion of energy directed to output ports can be flexibly adjusted, e.g., CHP in extraction condensing operation mode. The input-output characteristics become a multivariate quadratic function:
F = a ·P 2 + b ·Q 2 + c·PQ + d· P + e ·Q + f
(22)
Ag =
where F , P , Q denote the input gas, output electricity and output heat, respectively, and a, b, c, d , e , f are all constants. The characteristics of most SIMO energy converters can be modeled into a multivariate quadratic function, as shown in (22). We conduct a linear mapping to the input-output characteristics:
0s × B Is × s 0s × s 0s × B 0s × s Is × s
(1) If node g represents a nonlinear energy converter, A g is defined as follows:
2
(23)
2
1, branch b is secondary branch and is connected to input
where:
F P Q
1 0 0 = 0 1 c 2a 0 0 1
(26)
where the first B columns correspond to primary branches and the last 2s columns correspond to secondary branches. Then, consider the situation in which the EH has several nonlinear energy converters and several linear energy converters:
F = F1 + F2 F1 = a ·P + d ·P + f1 2 F = b ·Q + e~· Q + f 2
(25)
V =[ V T V T ]T
port k of node g
F P , Q
Ag (k , b) =
c2
b =b 4a ~ e = e cd 2a f = f1 + f2
- 1, branch b is secondary branch and is connected to output port k of node g 0, otherwise
(24)
After this mapping, the energy converter can be equivalently seen as two parallel independent ones. Curves F1 and F2 can be separately piecewise linearized in a similar way to SISO energy converters. A broader sense of “output energy” is introduced, i.e., P = P + (c 2a) Q , which is a linear combination of the original output energy types. Energy conversion with more output ports can be handled similarly.
(27) (2) If node g represents a linear energy converter, A g is defined as follows:
Ag = [ A g 0 ]
(28)
where A g is the port-branch incidence matrix of node g before piecewise linearization and the number of columns of the 0 matrix is equal to the length of V .
B. Introducing Splitter/Concentrator Through piecewise linearization, a nonlinear energy conversion or charging/discharging process in a SISO component is divided into several parallel processes with constant efficiencies. From a component level viewpoint, a nonlinear SISO component with input vgin and output vgout can be equivalent to a linear multi-input multi-output (MIMO) vgin, k (k = 1, 2, , s ) component with inputs and outputs out vg, k (k = 1, 2, , s ) , as shown in Fig. 3.
D. Formulating Converter Characteristic Matrices (1) SISO energy converters Using the method of piecewise linearization, a nonlinear SISO converter is converted into an s -input s -output converter. The converter characteristic matrix Hg has dimensions s × 2s and is defined as Fig. 3. Piecewise linearization of a nonlinear SISO component.
5
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
v
v
in g
in g
Primary Branch
Secondary Branch
Primary Branch
Primary Branch
Nonlinear Component η g (vgin )
Secondary Branch
Linear Component
Splitter
vgout
Concentrator
Primary Branch
vgout
η g ,k (k = 1, 2, , s ) Fig. 4. Splitter and concentrator.
Primary Branch
vgin
Primary Branch
Nonlinear Component
v
Primary Branch
vgout _ 2
Secondary Branch
Secondary Branch
in g
vgout _1
Primary Branch
Concentrator Linear Mapping
Linear Component
Splitter
vgout _1 vgout _ 2
Concentrator
Fig. 5. Piecewise linearization of a nonlinear SIMO energy converter.
follows:
Hg = [
Is× s ]
g, s × s
The converter characteristic matrices of linear energy converters remain the same.
(29)
where g, s × s = diag ( g ,1, g ,2, , g , s ) and conversion efficiency in the k -th segment.
g,k
E. Formulating energy balance matrices
is the constant energy
Given the port-branch incidence matrix A g and the component characteristic matrix Hg , the nodal energy balance matrix for node g can be calculated:
(2) SIMO energy converters Using piecewise linearization, a nonlinear 1-input b -output converter is converted into an s -input b·s -output converter. If the energy output of each output port is proportional to the energy input, the converter characteristic matrix Hg has dimensions b·s × (b + 1)· s and is defined as follows: 1, s × s
Hg =
2, s × s b, s × s
Is × s 0s × s 0s × s I s× s 0s × s
Suppose that EH has N nodes. The system energy balance matrix Z combines the nodal energy balance matrix of all nodes in the EH:
Z = Z1T Z 2T
0s × s 0s × s 0s × s I s×s
T
ZN
T
(33)
The energy balance matrix for each node defines the relationship between the energy flows in the branches connected to the node:
(30)
ZV = 0
If the proportion of energy directed to each output port can be flexibly adjusted, the converter characteristic matrix Hg has dimensions s × (b + 1)·s and is defined as follows:
(34)
F. Formulating splitter/concentrator characteristic matrices Splitter/concentrator characteristic matrix Wg of node g defines the relationship between the primary branches and secondary branches related to node g . First, consider the simplest case where the EH has only one nonlinear SISO energy conversion or storage component, which is
1, input port k is related to energy conversion process i Hg (i , k ) =
(32)
Zg = H g A g
1/ i , output port k is related to energy conversion process i 0, otherwise (31) 6
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
represented by node g , and its input and output ports are connected to only one branch respectively. A g is a matrix whose elements are the absolute value of those in A g (the port-branch incidence matrix of node g before piecewise linearization). If the number of segments in piecewise linearization is s , then:
11× s 01× s 01× s 11× s
A g = Ag 0
V =
(38)
W=
(40)
vgdisch
D g
T
(44)
C
D
(45)
1]
(46)
g
The input incidence matrix X is an m × B (B is the length of V ) matrix that relates the energy inputs to the branch energy flows. Similarly, the output incidence matrix Y is an n × B matrix that relates the energy outputs to the branch energy flows. We define: (47)
X=[ X 0 ], Y=[ Y 0 ]
where the number of columns of 0 matrix equals to the length of V . Then,
(41)
(48)
Vin = XV , Vout = YV
where and denote the charging and discharging efficiencies of component g , respectively. If the charging and discharging efficiencies of energy storage component g vary with its charging and discharging powers respectively, its efficiency curves can be also piecewise linearized in a similar way to the energy converters. In this case, secondary variables vgchar and vgdisch ,k ,k (k = 1, 2, s ) are introduced to represent the amounts of charging and discharging powers that fall in each segment and then Eg is: C g
Eg
H. Comprehensive energy flow equations for EH
The modeling of energy storage components is similar to that of SISO energy converters where the input port models charging and the output port discharging. The only difference is that energy can be stored in energy storage components and state of charge (SOC) is used to describe the ratio of stored energy to storage capacity. Let vgchar and vgdisch denote the charging and discharging powers of energy storage component g respectively, the rate of change of the SOC of this component is: C g
(43)
1
The formulation of the splitter/concentrator characteristic matrices of energy storage components is the same as that of energy converters.
G. Formulating matrices for energy storage components
Eg = vgchar ·
and augmented branch vector V :
0
T
Zg = H g A
(39)
WV = 0
g
where C = [ gC,1, gC,2, , gC, s ] and D = [1 gD,1, 1 gD,2, , 1 gD, s]. The nodal energy balance matrix for energy storage component g can be calculated as follows:
T T WG
V
Hg = [
The splitter/concentrator characteristic matrix for general situations can be formulated in a similar way. Suppose that EH has G nonlinear components in the EH. The system splitter/concentrator characteristic matrix W combines the nodal splitter/concentrator characteristic matrix of all nonlinear components: T W2
D g
Based on the above, the comprehensive energy flow equations for EH are: T
T T [ X T Y T Z T W T ]T V = [ Vin Vout 0T 0T ]
v2 vin ,1
v1
v3
CHP
ηe ηh
(49)
Except for the difference in formulation of the input/output incidence matrix and energy balance matrix, equations (49) have extra variables for energy flows in secondary branches and equations for
vout ,1
Gas
(42)
where A g can be established in a similar way to a SISO energy converter. Considering the energy balance of energy storage component g , i.e., (42), its component characteristic matrix is:
Then,
T W1
k=1
D g, k
D g, k
branch incidence matrix A
(37)
Wg [ V T V T ]T = Wg V = 0
vgdisch ,k
where and are the constant charging and discharging efficiencies of energy storage component g in the k -th segment. To model the SOC for energy storage component g under the proposed standardized matrix modeling framework, we introduce a virtual output port and a virtual branch Eg and construct augmented port-
(36)
A g Wg ]
C g, k
C g, k
where A g V denotes the energy flows in the primary branches related to node g . 11 × s (01 × s ) denotes a 1 × s vector whose elements are all ones (zeros). We define splitter/concentrator characteristic matrix Wg of node g as:
Wg = [
s
vgchar ,k · k=1
(35)
A g V = Wg V
Wg =
s
Eg =
v4
vout ,2
ηc WARG
Fig. 6. A CCHP-based EH. 7
v5
vout ,3
Electricity
Heat
Cooling
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
splitters/concentrators compared with (9). A simple example based on a combined cooling, heating and power (CCHP)-based EH shown in Fig. 6 is given to illustrate the proposed standardized matrix modeling approach of EH with variable efficiencies. The EH consists of a CHP unit and a water absorption refrigerator group (WARG). The input and output vectors of this CCHP-based EH are:
Vin = [vin,1]
(50)
Vout = [ vout ,1 vout ,2 vout ,3 ]T
(51)
and WARG can be calculated:
W2 = [ (52)
v 6 ]T
v5 v 1
are
A1 = [ A1 0 ] =
A2 =
0 1 0
03 × 5 I3 × 3 03 × 3 03 × 5 03 × 3 I3× 3
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
=
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
1 0 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0
e
0 0 0 0 0
H2 = [
c,3 × 3
h
I
]=
3×3
c,1
0 c,2
0 0
0
c,3
1 0 0 0 1 0 0 0 1
Gas
c,3
0 1 0
0 0 1
(59)
0 0 0 0 0 0
1 0
0 1 1 1 0 0 0 1 0 0 0 1 1 1
(60)
(61)
0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
(62)
0 0 1 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 1 0
0 0 0 1 0 0 0 0 0 0 1
0 0 0 0 0 0 c,1
0 0 0 0 0 0 0
0 0 1 0
c,2
0 0 0 0 0 0 0 0
0 1 0
1 0
c,3
0 0 0 0 0 0 1 0 0 0 1
0 0 0 0 0 0 0 1 0 0 1
0 0 0 0 0 0 0 0 1 0 1
v1 v2 v3 v4 v5 v1 = v2 v3 v4 v5 v6
vin,1 vout ,1 vout ,2 vout ,3 0 0 0 0 0 0 0
4. Standardized data structure This section presents a standardized data structure of the proposed modeling approach, which supports an automatic procedure for establishing the matrix model of arbitrary EH with variable energy efficiencies and related computations. The proposed data structure consists of two basic tables, i.e., node table and branch table, which act as the standardized input for MES information.
(57)
Given the port-branch incidence matrices and the component characteristic matrices, the nodal energy balance matrix for the CHP
vin ,1 v1
0
1 0 0
Eq. (63) provides a matrix model for EH with variable efficiencies. The proposed matrix modeling approach can be highly automated by computers and facilitates a highly efficient optimization in operation and planning.
(56)
0 0
0 0
(63)
0 0 0 0 0 1
1 0 0 1
A2 W2 ] =
0 1 0 0 1 0 0 0 0 0 0
h
(55)
e
c ,2
The input incidence matrix and output incidence matrix of the EH
1 0 0 0
The component characteristic matrices for the CHP and WARG are:
H1 = H1 =
0
c,1
0 0
(58)
Based on the above, the comprehensive energy flow equations for the EH are:
(54)
0 0 0 0 0 0 1 0 0 1 0 0
0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Y = [ Y 0] =
(53)
0 0 1
0 1
X = [ X 0 ] = [ 1 0 0 0 0 0 0 0 0 0 0]
Assume that the CHP is Node #1 and the WARG is Node #2, the port-branch incidence matrices for the CHP and WARG are:
1 0 0
1 0
The splitter/concentrator characteristic matrix of the WARG is:
The CHP is assumed to operate in backpressure mode (the production of electricity and heat are both proportional to the gas input) and produce electricity and heat with fixed efficiencies e and h , respectively. The WARG is assumed to convert its heat input into cooling with a variable efficiency c (v4) . Using the piecewise linearization method described above, the domain of function c (v4) can be divided into 3 segments (taking 3 segments as an example here) with fixed efficiencies c,1 , c,2 , and c,3 . As a result, the SISO WARG is converted to a 3-input 3output converter, and one splitter, one concentrator and 6 secondary branches vk (k = 1, 2, , 6) should be introduced. The EH after piecewise linearization is shown in Fig. 7. The branch vector V should be expanded:
V = [ V V ]T = [ v1
h
Z2 = H2 A2 =
The set of branches (energy flows) can be written as:
V = [ v1 v2 v3 v4 v5 ]T
e
Z1 = H1A1 =
(1) Node table
vout ,1
v2 CHP
v3
ηe ηh
vout ,2
v4 v1′ ~v3′
v4′ ~v6′
WARG
Splitter
Concentrator
ηc ,1 , ηc,2 , ηc,3 Fig. 7. A CCHP-based EH after piecewise linearization. 8
v5
vout ,3
Electricity
Heat
Cooling
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
In the node table, each row describes a single energy conversion/ storage component or an input/output of the EH, including:
Table 3 Parameters of the energy components in Fig. 8.
• Node No.: Numbering all nodes in EH, including energy conversion/storage components and inputs/outputs of EH, in order and starting from 1. • Node type: Identifying the type and operating mode of each node. • • •
1 2
3 1
3 4 5 6
−1 0 0 0
e A c
0 0 0 0
h B c
0 0 0 0
0 C c
0 0 0 0
Node Capacity
Segment Number
CCHP CWARG
0 3
0 0 0 0
0 0 0 0
Branch Type
Source Node No.
Sink Node No.
1 2 3 4 5
4 1 3 3 2
3 1 1 1 2
1 4 5 2 6
400
Out = 0.00003041· In3 + 0.01901·In2 + 0.2593· In * Out = 3·In
400 El: 300 Therm:420 900 800 (3.2 MWh)
A 2 B c = c ·v4 + c · v4 + WARG, respectively.
El: Out = 0.0001150·In2 + 0.2305·In
Therm: Out = 0.0001611·In2 + 0.3228·In ** Out = 0.8·In Charging Efficiency = 0.00005 In + 0.93 Discharging Efficiency = 0.00005 Out + 0.93
C c ).
CCHP and CWARG are capacities of CHP and
(2) Branch table In the branch table, each row describes one energy flow (only for primary branches) in the EH, including:
• Branch No.: Numbering all energy flows in the EH in order and starting from 1. • Branch Type: Identifying the type of each energy flow. For example, • •
the electricity, cooling, heat and gas flows are numbered as 1, 2, 3 and 4, respectively. Source Node No.: Node No. of the source node of the branch. Sink Node No.: Node No. of the sink node of the branch.
Taking the EH shown in Fig. 6 as example, the branch table is shown as Table 2. The component characteristic matrix Hg for each energy conversion or storage component is established based on the node table. The portbranch incidence matrix A g for each component, input incidence matrix X and output incidence matrix Y are established based on the branch table. Splitter/concentrator characteristic matrix Wg is established based on both the node table and branch table.
Table 2 Branch table of the EH in Fig. 6. Branch No.
CERG
* Part-load performance of the existing chiller with specific settings of condensing pressure and outdoor temperature; ** Part-load performance of natural gas-powered microturbine with a heat exchanger.
Table 1 Node table of the EH in Fig. 6. Energy Efficiency Parameter
Performance [5,26,31,33]
AB HS
The inputs/outputs of EH do not have data of Energy Efficiency Parameter, Node Capacity and Segment Number. Taking the EH shown in Fig. 6 as example, the node table is shown as Table 1. The energy conversion characteristics of CHP and WARG are described as v2 = e ·v1, (v3 + v4) = h · v1, v5 = cA ·v43 + cB ·v42 + cC · v4 (i.e.,
Node Type
Capacity (kW)
HP CHP
For example, −1 for input, 0 for output, 1 for WARG, 2 for CHP unit in extraction condensing operation mode, 3 for CHP unit in back pressure operation mode and so on. A component library can be constructed for index. Energy Efficiency Parameter: Providing information on energy conversion characteristics of energy conversion components and energy charging/discharging characteristics of energy storage components. The number of parameters depends on the type of node. Node Capacity: Providing information on capacity of each energy conversion or storage component, which is used when modeling the operating constraints of the EH. Segment Number: Determining the number of segments for piecewise linearization. 0 is for energy conversion and storage components with constant efficiencies.
Node No.
Converter
5. Case study This section illustrates the application of the proposed standardized matrix modeling approach in optimal operation of EH. A. System description Fig. 8 illustrates an EH which consists of a CHP in backpressure
Fig. 8. An EH with three nonlinear energy components. 9
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
operation mode (heat-to-power ratio is a constant), a compression electric refrigerator group (CERG), an electric heat pump (HP), an auxiliary boiler (AB) and a heat storage component (HS). The energy inputs are electricity and gas, while the outputs consist of cooling, electricity and heat. The CHP, CERG and HS are considered to exhibit considerable efficiency variations with their input or output powers changing. The electric efficiency of the CHP in partial loads can degrade more than 20% compared to the full-load efficiency [30]. The specific type of CERG in Fig. 8 is an air-cooled centrifugal chiller whose highest coefficient of performance (COP) occurs at a part-load ratio of 0.71–0.84 rather than at full load [31]. The efficiencies of the HP and AB do not differ significantly under part-load conditions [26,32]; therefore, their efficiencies are assumed to be constant. Table 3 lists the parameters of the energy components of this system. The hourly electricity, heat and cooling demand patterns refer to a hospital site in a Mediterranean area [5] and are indicated in Fig. 9. The gas price is set at 20 euro/MWh and is considered constant during the daily analysis period [5]. The variation in hourly electricity prices for the selected day corresponding to a real case in Italy [5] is also shown in Fig. 9. Using the proposed modeling approach described above, the piecewise linearized EH with splitters/concentrators is shown in Fig. 10 (taking 2 segments as an example here). Given the topology of EH, parameters of energy components and the pre-defined number of segments for piecewise linearization (taking 2 segments as an example here), the computer can automatically formulate the comprehensive energy flow equations of EH as shown in (64).
given time horizon. By replacing the comprehensive energy flow equations (9) with (49) and introducing piecewise linearization constraints (e.g., (20)–(21)), a revised optimal operation model based on the model in [22] is used. The decision variables of this problem are the continuous variables V , the continuous variables Vin , the continuous variables E (state of charge of energy storage components) and the binary variables Ug in each time period, making the problem a mixed integer linear programming (MILP) problem. C. Optimization results and sensitivity analysis The objective of the EH operational optimization in this case study is to minimize the daily operating cost of the whole system. Three cases are considered: (1) Case 1: The exact performances of CHP, CERG and HS are considered, making the optimal operation problem a NP problem; (2) Case 2: The efficiencies of the CHP, CERG and HS are approximated as constants (e.g., efficiencies of converters are calculated by dividing maximum output by maximum input), making the optimal operation problem a linear programming problem; (3) Case 3: The proposed standardized matrix modeling approach is used. The CHP, CERG and HS are piecewise linearized, making the optimal operation problem an MILP problem. In addition, subcases of Case 3 with different numbers of segments in piecewise linearization are considered. When the number of segments is large en-
(64)
B. Multi-period optimal operation model for EH
ough (e.g., 300), the optimization results are considered accurate and taken as the reference compared with the other optimization results.
The optimal energy flow problem of a single EH determines how energy flows should be dispatched among different components to minimize the cost or maximize the profit of the whole system over a
Case studies are executed in MATLAB with YALMIP [34] on a PC 10
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
Fig. 9. Hourly load patterns and electricity prices of the selected day.
with an Intel Core i7-7500U 2.70 GHz CPU and 16 GB RAM. Cases 2 and 3 are solved using GUROBI. Case 1 is solved using fmincon with the sequential quadratic programming (SQP) method. The maximum number of iterations is set at 3000. The optimization results of Case 2 are used as initial values for Case 1. The optimal operating costs, computation times and relative errors of optimal values for these three cases are summarized in Table 4 and Fig. 11, in which s denotes the number of segments in piecewise linearization. For Case 2, the relative error of the optimal operating cost is over 13%, which indicates that the constant efficiency approximation does not adequately model the EH. With an increase in the number of nonlinear components, nonlinearity of components and complexity of EH, the constant efficiency model of EH results in larger errors due to approximations, making its optimization results impractical. With the proposed standardized matrix modeling approach, the relative error of Case 3 is below 1.0% and 0.1% when s is equal to 12 and 36, respectively. With no evident calculation time increase (less than 0.2 s) compared to Case 2, the proposed approach reduces the approximation error from 13.7% to 0.1% in Case 3. Using the results provided by Case 2 as initial values, Case 1 obtains much more accurate results than Case 2. Nevertheless, the NP problem cannot guarantee a global optimal solution, and its solution is time-consuming. The proposed standardized
v2
Electricity
Gas
Splitter
vin ,1
v1′, v2′
CERG
v7′ , v8′
Table 4 Optimal operating costs, computation times and relative errors for the three given cases. Case Case Case Case Case Case Case Case Case Case Case Case Case Case Case Case
(s (s (s (s (s (s (s (s (s (s (s (s (s
= = = = = = = = = = = = =
2) 4) 6) 8) 10 ) 12 ) 20 ) 36 ) 50 ) 76 ) 100 ) 200 ) 300 )
Computation Time (s)
Relative Error* (%)
313.628 264.440 280.842 289.797 297.182 301.089 303.020 304.222 305.683 306.230 306.325 306.372 306.393 306.416 306.420
292.73 0.43 0.47 0.61 0.60 0.50 0.50 0.49 0.53 0.58 0.67 0.88 1.29 2.73 3.39
2.35 13.70 8.35 5.43 3.01 1.74 1.11 0.72 0.24 0.06 0.03 0.02 0.01 0.001 –
* Relative errors for each case/subcase are calculated according to the optimization results of Case 3 (s = 300 ).
Concentrator
vout ,1
v6
v1 v3
vin ,2
1 2 3 3 3 3 3 3 3 3 3 3 3 3 3
Optimal Value (Euro)
v4
HP
′ v9′ , v10
Concentrator
CHP
Splitter
v3′ , v4′
′ , v12 ′ v11
Concentrator
v11 v7
v8
Splitter
v5′ , v6′
AB
v10
Fig. 10. The piecewise linearized EH with splitters/concentrators. 11
HS
Electricity
′ , v14 ′ v13 Concentrator
v9
v5
vout ,2
Cooling
v12
vout ,3
Heat
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
Fig. 11. Computation times and relative errors for different numbers of segments in Case 3.
matrix modeling approach makes the optimal operation problem an MILP problem that can be efficiently solved by many existing commercial software (e.g., GUROBI, CPLEX, MOSEK). It should also be noted that the inefficiency of the operation is better approximated when the number of segments in piecewise linearization increases, which leads to a higher operating cost. As shown in Fig. 11, relative error decreases rapidly when s increases from 2 to 20, and computation time increases relatively slowly when s increases from 30 to 70; consequently, the optimal number of s is between 30 and 70 considering the balance between nonlinear approximation accuracy and calculation efficiency. Figs. 12 and 13 show the operating states of each energy component in the EH over a 24-hour period obtained in Case 1, Case 2 and Case 3
(s = 100 ). The electricity and heat demands are provided by the CHP, HS and AB during the daytime because the price of electricity is high. Load demands and electricity price decrease during the night; therefore, it is inefficient for the CHP to work at a low load level, and it is more cost-effective to purchase electricity from the distribution system to meet the electricity demand and provide power for the HP to produce heat. HS contributes to reducing the heating cost by storing heat when electricity price is low during the night and discharging it in the daytime. The AB is used as a supplemental source for heat supply. The constant efficiency approximation in Case 2 leads to a wrong and uneconomical operation strategy as shown in Fig. 12: (1) The CHP still works at a relatively low load level during the night without regard to its part-load performance. (2) The electricity input of the CERG
CERG: Case 3 HP: Case 3 CHP: Case 3 AB: Case 3 HS: Case 3
CERG: Case 2 HP: Case 2 CHP: Case 2 AB: Case 2 HS: Case 2
Fig. 12. Operating states of the energy components in the EH of Case 2 and Case 3 (s = 100 ). 12
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
CERG: Case 3 HP: Case 3 CHP: Case 3 AB: Case 3 HS: Case 3
CERG: Case 1 HP: Case 1 CHP: Case 1 AB: Case 1 HS: Case 1
Fig. 13. Operating states of the energy components in the EH of Case 1 and Case 3 (s = 100 ).
decreases because the efficiency is supposed to remain high at any load level. (3) The HS discharges heat with higher power between 10:00 and 14:00 without considering the increasing power loss with increasing discharging power. The NP problem of Case 1 results in a sub-optimal operation strategy by obtaining the local optimal solution, which is highly impacted by the initial values of decision variables. Comparing Fig. 12 and Fig. 13, it can be found that the operation strategy obtained by Case 1 is very similar to that of Case 2, which is because the optimization results of Case 2 are used as initial values for Case 1. Though Case 1 considers the exact performances of CHP, CERG and HS, its optimization results are misled by Case 2 which ignores the part-load performances. In addition, between 8:00 and 22:00, the input powers of CHP and CERG in Case 1 are almost the same as those in Case 3 (s = 100 ), but different from those in Case 2. This indicates that the part-load performances of components are well modeled by the proposed approach used in Case 3 (s = 100 ).
condition for the operation model and it can be obtained in advance by the weather forecast. Different efficiency curves can be set in advance for different temperature differences during different time periods in the optimal operation problem. In other words, the efficiency variation caused by external factor, e.g., environmental factor, can be obtained as conditional efficiency in advance based on the forecast of the factor. In addition, for the optimal operation problem, the external factor does not depend on the operation strategy. Therefore, the operation model under a certain external factor, e.g., temperature difference, still has a linear formulation. By contrary, the efficiency variation and nonlinearity caused by partial load is the key difficulty in the optimal operation problem, which is handled using the proposed modeling approach. The piecewise linearization process can be further improved using self-adaptive technique, i.e., the widths and number of piecewise linearization segments can be unevenly chosen according to the nonlinearity of efficiency curves. For example, the widths of different segments for one component can be different. Narrower segments are used to approximate the part with higher nonlinearity on the curve, which further improves the approximation accuracy. In addition, different energy conversion and storage components have different nonlinearities and they can use different numbers of segments in piecewise linearization to reduce the number of segment variables without sacrificing modeling accuracy. The proposed approach can effectively model efficiency curves with all kinds of shapes, irrespective of whether they are convex or nonconvex. E.g., linear, quadratic, cubic, higher order polynomial and exponential efficiency curves can all be easily handled by the proposed piecewise linearization based modeling approach.
6. Discussion The proposed matrix modeling approach of EH can be effectively applied to a more complex case because the proposed approach can be highly automated by computers. In practical MES, more linear and nonlinear components, more complicated system structure and a longer time horizon would lead to an increase in the size of the matrices and the number of continuous and binary variables in the operation model. The matrices obtained, though very large in scale, can be easily handled by computers using effective techniques such as sparse technique. Binary variables in the proposed approach are tightly constrained by (20)–(21). Once the value of one binary variable for one segment is determined, the values of the other binary variables for the other segments can be immediately determined. So the number of binary variables does not significantly impact the calculation time which is proved in the case study (e.g., Fig. 11). The energy efficiencies of components may depend on several parameters besides input or output power, e.g., temperature difference; however, it does not introduce nonlinearity to the operation model. For example, the efficiency of heat pump varies with the temperature difference, e.g., the heating efficiency is high when temperature difference is small during the daytime while it is low when temperature difference is large at night. However, the outdoor temperature is a boundary
7. Conclusion The variable energy efficiencies of the energy conversion and storage components in MES complicate the modeling of EH in two ways: (1) The nonlinearity caused by variable efficiencies complicates the analysis and optimization of EH; (2) The nonlinear energy components can hardly be automatically modeled by computers. To jointly address these issues, we propose a standardized matrix modeling approach to automatically model EH with variable energy efficiencies. Using piecewise linearization, a nonlinear energy component is converted into a “splitter, linear energy component, concentrator” three-layer structure. 13
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
The nonlinear energy conversion and storage relationship in EH can thus be further modeled under a linear modeling framework using matrices. Such modeling approach facilitates an automatic modeling of an arbitrary EH with various kinds of nonlinear energy components. Using the proposed modeling approach, the optimal operation of a fivecomponent EH with energy storage shows a great improvement on the accuracy and calculation efficiency compared to the exiting constant efficiency model and nonlinear model, respectively. Future work includes extending the modeling framework to consider the dynamics of energy conversion in different time scales.
[8] Tian X, Zhao R. Energy network flow model and optimization based on energy hub for big harbor industrial park. J Coastal Res 2015;73:298–303. [9] Krause T, Kienzle F, Liu Y, Andersson G. Modeling interconnected national energy systems using an energy hub approach. In: Proc IEEE powertech, 2011, Trondheim, Norway; Jun. 2011. p. 1–7. [10] Guler B, Çelebi E, Nathwani J. A ‘Regional Energy Hub’for achieving a low-carbon energy transition. Energy Policy 2018;113:376–85. [11] Mohammadi M, Noorollahi Y, Mohammadi-ivatloo B, Hosseinzadeh M, Yousefi H, Khorasani ST. Optimal management of energy hubs and smart energy hubs–a review. Renew Sustain Energy Rev 2018;89:33–50. [12] Wang Y, Zhang N, Zhuo Z, Kang C, Kirschen DS. Mixed-integer linear programmingbased optimal configuration planning for energy hub: Starting from scratch. Appl Energy 2018;210:1141–50. [13] Sadeghi H, Rashidinejad M, Moeini-Aghtaie M, Abdollahi A. The energy hub: an extensive survey on the state-of-the-art. Appl Therm Eng 2019:114071. [14] Geidl M, Andersson G. Optimal power flow of multiple energy carriers. IEEE Trans Power Syst 2007;22(1):145–55. [15] Biglia A, Caredda FV, Fabrizio E, Filippi M, Mandas N. Technical-economic feasibility of CHP systems in large hospitals through the Energy Hub method: the case of Cagliari AOB. Energy Build 2017;147:101–12. [16] Cheng Y, Zhang N, Wang Y, Yang J, Kang C, Xia Q. Modeling carbon emission flow in multiple energy systems. IEEE Trans Smart Grid July 2019;10(4):3562–74. [17] Waibel C, Evins R, Carmeliet J. Co-simulation and optimization of building geometry and multi-energy systems: Interdependencies in energy supply, energy demand and solar potentials. Appl Energy 2019;242:1661–82. [18] Geidl M, Andersson G. Optimal coupling of energy infrastructures. Proc IEEE Lausanne Power Tech 2007:1398–403. [19] Almassalkhi M, Hiskens I. Optimization framework for the analysis of large-scale networks of energy hubs. In: Proc 17th power syst comput conf; 2011. p. 1–7. [20] Beccuti G, Demiray T, Batic M, Tomasevic N, Vranes S. Energy hub modelling and optimisation: an analytical case-study. In: PowerTech 2015 IEEE Eindhoven; 2015. p. 1–6. [21] Wang Y, Cheng J, Zhang N, Kang C. Automatic and linearized modeling of energy hub and its flexibility analysis. Appl Energy Feb. 2018;211:705–14. [22] Wang Y, Zhang N, Kang C, Kirschen DS, Yang J, Xia Q. Standardized matrix modeling of multiple energy systems. IEEE Trans Smart Grid 2019;10(1):257–70. [23] Fabrizio E, Filippi M, Virgone J. An hourly modelling framework for the assessment of energy sources exploitation and energy converters selection and sizing in buildings. Energy Build 2009;41(10):1037–50. [24] Evins R, Orehounig K, Dorer V, Carmeliet J. New formulations of the ‘energy hub’ model to address operational constraints. Energy 2014;73:387–98. [25] Sakti A, Gallagher KG, Sepulveda N, Uckun C, Vergara C, de Sisternes FJ, et al. Enhanced representations of lithium-ion batteries in power systems models and their effect on the valuation of energy arbitrage applications. J Power Sources 2017;342:279–91. [26] Milan C, Stadler M, Cardoso G, Mashayekh S. Modeling of non-linear CHP efficiency curves in distributed energy systems. Appl Energy 2015;148:334–47. [27] Almassalkhi MR, Towle A. Enabling city-scale multi-energy optimal dispatch with energy hubs. In: Presented at 2016 power syst comput conf (PSCC), Genoa; 2016. p. 1–7. [28] Holjevac N, Capuder T, Zhang N, Kuzle I, Kang C. Corrective receding horizon scheduling of flexible distributed multi-energy microgrids. Appl Energy 2017;207:176–94. [29] Shabanpour-Haghighi A, Seifi AR. Energy flow optimization in multi-carrier systems. IEEE Trans. Ind. Informat. 2015;11(5):1067–77. [30] Bianchi M, Pascale AD, Melino F, Peretto A. Performance prediction of micro-CHP systems using simple virtual operating cycles. Appl Therm Eng 2014;71(2):771–9. [31] Yu FW, Chan KT. Part load performance of air-cooled centrifugal chillers with variable speed condenser fan control. Build Environ 2007;42(11):3816–29. [32] Anonymous. Boiler efficiency guide: facts about firetube boilers and boiler efficiency. CleaverBroocks; Mar. 2010. [Online]. Available: http://cleaverbrooks.com/ reference-center/insights/Boiler%20Efficiency%20Guide.pdf. [33] Zhou Z, Liu P, Li Z, Ni W. An engineering approach to the optimal design of distributed energy systems in China. Appl Therm Eng 2013;53:387–96. [34] Lofberg J. YALMIP: a toolbox for modeling and optimization in MATLAB. In: presented at proc IEEE int'l symp comput aided control syst des, Taipei; 2004. p. 284–289.
CRediT authorship contribution statement Wujing Huang: Methodology, Software, Validation, Investigation, Data curation, Writing - original draft, Visualization. Ning Zhang: Conceptualization, Writing - review & editing, Supervision. Yi Wang: Writing - review & editing. Tomislav Capuder: Writing - review & editing. Igor Kuzle: Writing - review & editing. Chongqing Kang: Writing - review & editing. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgement This work was supported in part by China-Croatia InterGovernmental S&T Cooperation Project Flexible Integrated Demand Response in Multi-Energy System (FIRE), National Natural Science Foundation of China (No. 71961137004), Tsinghua University Initiative Scientific Research Program (20193080026), and EU-CHINA Research and Innovation Partnership Project Instigation of Research and Innovation Partnership on Renewable Energy, Energy Efficiency and Sustainable Energy Solutions for Cities (IRES-8), contract number ICI+/2014/347-910. References [1] Mancarella P. MES (multi-energy systems): an overview of concepts and evaluation models. Energy 2014;65:1–17. [2] Holjevac N, Capuder T, Kuzle I. Adaptive control for evaluation of flexibility benefits in microgrid systems. Energy 2015;92(4):487–504. [3] Ceseña EAM, Capuder T, Mancarella P. Flexible distributed multienergy generation system expansion planning under uncertainty. IEEE Trans Smart Grid 2016;7(1):348–57. [4] Geidl M, Koeppel G, Favre-Perrod P, Klckl B, Andersson G, Frhlich K. Energy hubs for the future. IEEE Power Energ Mag 2007;5(1):24–30. [5] Chicco G, Mancarella P. Matrix modelling of small-scale trigeneration systems and application to operational optimization. Energy Mar. 2009;34(3):261–73. [6] Rastegar M, Fotuhi-Firuzabad M, Lehtonen M. Home load management in a residential energy hub. Elect Power Syst Res 2015;119:322–8. [7] Bahrami S, Safe F. A financial approach to evaluate an optimized combined cooling heat and power system. Energy Power Eng J 2013;5(5):352–62.
14
Contents lists available at ScienceDirect
Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes
Matrix modeling of energy hub with variable energy efficiencies a
a,⁎
b
c
c
a,⁎
Wujing Huang , Ning Zhang , Yi Wang , Tomislav Capuder , Igor Kuzle , Chongqing Kang a b c
T
Department of Electrical Engineering, Tsinghua University, China Department of Information Technology & Electrical Engineering, ETH Zurich, Switzerland Faculty of Electrical Engineering and Computing, University of Zagreb, Croatia
ARTICLE INFO
ABSTRACT
Keywords: Multi-energy systems Energy hub Matrix modeling Piecewise linearization Variable efficiency Part-load performance Operation optimization
The modeling of multi-energy systems (MES) is the basic task of analyzing energy systems integration. The variable energy efficiencies of the energy conversion and storage components in MES introduce nonlinearity to the model and thus complicate the analysis and optimization of MES. In this paper, a standardized matrix modeling approach is proposed to automatically model MES with variable energy efficiencies based on the energy hub (EH) modeling framework. Piecewise linearization is used to approximate the variable energy efficiencies; as a result, a component with variable efficiency is equivalent to several parallel components with constant efficiencies. Splitters and concentrators are proposed as standardized components to facilitate the split and merge of energy flows imposed by piecewise linearization. The nonlinear energy conversion and storage relationship in EH can thus be further modeled under a linear modeling framework using matrices. Such matrix modeling approach makes the modeling of an arbitrary EH with nonlinear energy components highly automated by computers. The proposed modeling approach can further facilitate the operation and planning optimization of EH with variable efficiencies. Case studies are executed in MATLAB and presented to show how the nonlinear approximation accuracy and calculation efficiency can be balanced using the proposed model in the optimal operation of EH. For the optimal operation of a five-component EH with energy storage, the proposed approach reduces the approximation error from 13.7% to 0.1% with only 0.2 s computation time increase compared to the exiting constant efficiency model.
1. Introduction The interactions between electricity, gas and heat/cooling systems enable higher energy efficiency and utilization of intermittent renewable energy sources [1–3]. The modeling of multi-energy systems (MES) has gained increasing importance in recent years as the foundation of the operation and planning of MES, which can unlock the flexibility of shifting across multiple energy sectors and result in reduced costs and lower emissions as compared to separate operation and planning of energy systems. The coupling relationship among different energy sectors in MES and the interaction characteristics between MES and external networks/loads are the key issues that the modeling approach needs to address. The energy hub (EH) concept was introduced to model MES as a unit in which multiple energy carriers can be converted, conditioned and stored. EHs import energy at their input ports, which are connected to external energy infrastructures, e.g., power grid and gas network. They provide required energy services, such as electricity, heat and cooling, at their output ports [4]. The compactness and effectiveness of EH-based model facilitate the calculation of energy
⁎
flows in MES of different levels, such as trigeneration systems [5], residential consumers [6], commercial buildings [7], industrial parks [8], national energy systems [9] and interconnected countries [10]. In addition, EH can be integrated with smart technologies especially smart meters, becoming smart energy hub (SEH), to handle a large volume of information in MES [11]. While a substantial number of papers have discussed the optimal planning [12,13], operation [14,15] and carbon emission [16,17] of MES based on the EH model, very few papers have studied the standardized and automatic modeling of MES. Based on the concept of EH, a nonlinear framework for modeling of energy systems with conversion and storage of multiple energy carriers is proposed in [18]. In this modeling framework, the converter and storage coupling matrices are manually formulated according to the given EH layout. Reference [5] proposes an approach to automatically generate coupling matrices for small-scale trigeneration systems using dispatch factors [13]. The formulation involves products of dispatch factors and energy flows, leading to nonlinear rather than linear optimization problems. As a result, the approach is more applicable to small-scale and simple-
Corresponding authors at: Room 3-206, West Main Building, Tsinghua University, Haidian, Beijing, 100084, China. E-mail addresses: [email protected] (N. Zhang), [email protected] (C. Kang).
https://doi.org/10.1016/j.ijepes.2020.105876 Received 16 July 2019; Received in revised form 20 December 2019; Accepted 19 January 2020 0142-0615/ © 2020 Elsevier Ltd. All rights reserved.
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
Nomenclature AB CCHP CERG CHP EH HP
HS MES MILP MIMO NP PWL SIMO SISO
Gas-fired auxiliary boiler Combined cooling, heat and power Compression electric refrigerator group Combined heat and power Energy hub Electric heat pump
structure MES. Reference [19] introduces a linear and automatic modeling tool, Hubert, which leverages a concise ASCII description format for EH with an energy “input, input storage, converter, output storage, output” structure. Reference [20] improves Hubert by adding an “export” branch between “converter” block and “output storage” block to facilitate selling PV-generated energy to the external grid. Reference [21] proposes the idea of a linear coupling relationship of an EH by introducing the concept of augmented variables. The standardized procedure for modeling the dispatch-factor-free coupling matrix is then designed to facilitate the computerized calculation for arbitrary EH configurations. Reference [22] proposes an automatic and standardized matrix modeling method for EH with arbitrary configuration. The method uses graph theory to cast the topology and the characteristics of the energy converters into a matrix form and uses linear equations for the dispatch of energy flows without using dispatch factors. However, all the energy efficiencies are regarded as constants in the above research. In reality, energy efficiencies of energy conversion and storage components vary with their operating conditions. The products of variable efficiencies and energy flows bring nonlinearity to the components and systems. Energy conversion efficiencies of converters can reduce significantly when they do not work in rated condition; therefore, the constant efficiency approximation model of MES with high nonlinearity is very likely to be inaccurate and impractical. There are already some methods taking into consideration the variable efficiencies of energy conversion and storage components. In [23], nonlinear energy converters are directly modeled with highly nonlinear part-load efficiency curves. The resulting model is a nonlinear programming (NP) problem which gives no guarantee that the global optimum can be achieved and few numerically robust solvers for NP problems exist. To avoid NP problems, several linearization approaches are used to model nonlinear converters. The first one is step-wise approximation based on binary variables which guarantees that only one operating point is selected in each time step (e.g., [24,25]). This formulation carries the disadvantage that it only allows single tabled values to be selected, with no interpolation between them. The second one is piecewise linear (PWL) approximation using special-ordered-set (SOS) variables (e.g., [26]). A special ordered set is a set of consecutive variables in which not more than one (named SOS1) or two (named SOS2) adjacent members may be non-zero in a feasible solution. By using SOS2 variables as weight variables, any operating point can be written as a linear combination of two pre-defined operating points. Although these methods are able to address variable efficiencies, the formulation can hardly be automated by computers. Reference [27] improves Hubert by enabling the modeling of nonlinear converters using PWL approximation. The operating points are defined as the sum of a set of consecutive variables related to different PWL segments, which has been widely used to approximate quadratic fuel cost and transmission loss in power system operation and planning. However, this model can still only deal with the proposed structure of MES where the location of energy storage components is fixed and only one energy conversion block is considered in an EH. In addition, these methods cannot handle the nonlinear energy converters with multiple outputs and complex operation characteristics, e.g., combined heat and power (CHP) in extraction condensing operation mode in which the ratio between electricity and heat production is adjustable within a certain
Heat storage Multi-energy systems Mixed integer linear programming Multi-input multi-output Nonlinear programming Piecewise linear Single-input multi-output Single-input single-output
range. In this paper, we propose a standardized matrix modeling approach for EH with variable energy efficiencies. Using piecewise linearization, the nonlinear energy conversion or charging/discharging process in an energy component is divided into several linear processes with constant efficiencies. Splitters and concentrators are proposed as standardized components to facilitate the split and merge of energy flows imposed by piecewise linearization; as a result, a nonlinear energy component is converted into a “splitter, linear energy component, concentrator” three-layer structure. The nonlinear energy conversion and storage relationship in EH can thus be further modeled under a linear and matrixbased modeling framework. Such approach facilitates an automatic modeling of an arbitrary EH with nonlinear energy components. Compared to the constant efficiency model and the nonlinear model, the proposed standardized matrix model is proven to significantly improve the accuracy and calculation efficiency, respectively, of the operation optimization of EH. The contributions of this paper are threefold: (1) Providing a standardized matrix modeling approach for EH with nonlinear energy conversion and storage components, which makes the modeling of an arbitrary EH highly automated by computers. Such technique fills the gap that current nonlinear/linearized modeling approaches can hardly automatically formulate an arbitrary EH with nonlinear energy components. (2) Proposing a generalized linearized model that is capable of accurately modeling any energy conversion or storage component with constant or variable efficiencies under one linear modeling framework. Especially, the generalized linearized model is able to handle the nonlinear energy converters with multiple outputs and complex operation characteristics, e.g., CHP in extraction condensing operation mode. (3) Providing analysis of how the nonlinear approximation accuracy and calculation efficiency can be balanced using the proposed model in the optimal operation of EH. The rest of this paper is organized as follows: Section 2 provides the basic backgrounds of the standardized matrix model of EH. Section 3 proposes the standardized matrix modeling approach to tackle the challenges of variable energy efficiencies of energy components. Section 4 presents the standardized data structure of the proposed modeling approach. Section 5 conducts a comprehensive case study and applies the proposed modeling approach to the operational optimization problem. Section 6 provides a discussion on the proposed modeling approach. Section 7 draws the conclusion. 2. Background: Standardized matrix model of EH with constant efficiencies First, a standardized matrix model of EH with constant energy efficiencies [22] is briefly introduced to facilitate understanding of the proposed matrix model of EH with variable efficiencies and the relationship between these two models. From the viewpoint of graph theory, each energy conversion or storage component in an EH can be treated as a node, and each energy 2
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
flow to or from a component can be treated as a branch, as shown in Fig. 1. The inputs and outputs of an EH are treated as special nodes. They are also related to one or more energy flows but not convert or store energy. Each node has several input and output ports through which it exchanges energy with other nodes.
conversion equation of the EH can be obtained:
ZV = 0 D. Comprehensive energy flow equations for EH
The m -dimensional energy input vector Vin and the n -dimensional energy output vector Vout are defined if an EH has m inputs and n outputs. The input incidence matrix X is an m × B matrix that relates the energy inputs to the branch energy flows. Similarly, the output incidence matrix Y is an n × B matrix that relates the energy outputs to the branch energy flows:
A. Formulating port-branch incidence matrices A port-branch incidence matrix A g of node g defines the connections between the ports of node g and all the branches in the EH. If node g has K g ports and the EH has B branches, the port-branch incidence matrix A g has dimensions K g × B and is defined as follows:
Ag (k , b) =
1, branch b is connected to input port k of node g - 1, branch b is connected to output port k of node g 0, otherwise
(1)
where k is the serial number of input and output ports of a node.
(7) (8)
Given the port-branch incidence matrix A g and the converter characteristic matrix Hg , the nodal energy conversion matrix Zg for node g can be calculated: (3)
Zg = Hg Ag
The energy conversion matrix Zg for node g defines the relationship between the energy flows in the branches connected to node g . This matrix has as many rows as there are energy conversion processes in node g and as many columns as there are branches in the EH. Suppose that an EH has N nodes. The system energy conversion matrix Z combines the nodal energy conversion matrix of all nodes in the EH: T
ZTN ]
Energy conversion efficiency is taken as an example here. (1) Type 1: SISO energy converters
The energy flows in all branches form a B × 1 vector V . According to the energy conversion characteristic of each node, the energy
v1
Port
(9)
A. Piecewise linearization of energy efficiency
(4)
Node
T
T T Vout 0T ] Y T Z T ]T V = [ Vin
When the energy efficiencies of energy conversion and storage components vary with their input or output powers, the constant efficiency model of EH is no longer accurate. Fig. 2 shows the output power curve of a nonlinear converter and its linear approximation which has a significant error at a certain load level. The approximations of variable efficiencies of production units, particularly CHP units, have a high impact on short-term schedules [28,29]. The products of variable efficiencies and energy flows also bring about nonlinearity to the optimization problem. In this paper, EH with variable efficiencies is modeled using a piecewise linearization method. Through piecewise linearization, a nonlinear energy conversion or storage component is divided into several equivalent components with constant efficiencies. The modeling of two types of energy converters, i.e., single-input single-output (SISO) and single-input multi-output (SIMO) energy converters, and the modeling of energy storage components are presented. The modeling of energy storage components is similar to that of SISO energy converters where the input port models charging and the output port discharging.
C. Formulating energy conversion matrices
Gas
1, branch j is connected to output port i 0, otherwise
3. Proposed standardized matrix modeling approach
is the energy conversion efficiency of energy conversion pro-
vin ,1
Y (i , j ) =
The comprehensive energy flow equations of the EH are:
(2)
Input
(6)
[ XT
s , input port k is related to energy conversion process s Hg (s, k ) = 1, output port k is related to energy conversion process s 0, otherwise
Z = [ Z1T Z2T
1, branch j is connected to input port i 0, otherwise
Vin = XV, Vout = YV
The converter characteristic matrix Hg of node g defines the energy conversion characteristics of the node. If node g has Sg energy conversion processes and K g ports, the converter characteristic matrix Hg has dimensions Sg × K g . If each energy conversion process in node g is related to only one input and one output, then Hg is defined as follows
s
X (i , j ) =
The input incidence and output incidence equations are:
B. Formulating converter characteristic matrices
where cess s .
(5)
Electricity
v2
Branch
v3
#1
v4
Heat
#2
v5
Cooling
Fig. 1. Definitions of basic elements of an EH in terms of graph theory. 3
Output
vout,1
Output
vout,2
Output
vout,3
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
Fig. 2. Nonlinear component variable efficiency linear approximation and piecewise linear approximation.
Suppose that node g has a single input vgin and a single output vgout and that its energy conversion efficiency g (vgin) is a function of the input power:
vgout =
in in g (vg )·vg
k
where v gin, k is the range of the k -th segment and
k= 1
vgin, k
vgin
=
k 1 i=1
v gin, i 0,
and vgin, k (k = 1, 2,
,
k 1 i=1
vgin
k i=1
v gin, i
v gin, k = v gin . Secondary
vgout ,k =
k 1 i=1
g,k
v gin, i
< vgin < v gin, i
k i=1
(14)
s ) are introduced and
vgout ,k
(15)
In the k -th (k = 1, 2, , s ) segment, the variable efficiency g (vgin) is replaced by a fixed efficiency g , k determined by the secant line of the power output curve in this segment:
variables vgin, k (k = 1, 2, s ) are introduced to represent the amounts of input power that fall in each segment:
v gin, k , vgin
,s
k=1
(11) s
v gin, i , k = 1, 2,
i=1
s
vgout =
[0, v gin,1], [v gin,1, v gin,1 + v gin,2], [v gin,1 + v gin,2, v gin,1 + v gin,2 + v gin,3], +v gin, s]
i=1
k
Similarly, secondary variables vgout , k (k = 1, 2, subject to
As shown in Fig. 2, the range of vgin (denoted by v gin ) can be divided into s segments:
+ v gin, s 1, v gin,1+
v gin, i ·
g
i=1
(10)
, [v gin,1+
k
v gout ,i =
in g , k ·vg , k
, k = 1, 2,
,s
in = v gout , k vg,k
(17)
Binary variables ug ,1, ug ,2 , continuity of vgin and vgout :
v gin, i
ug,1·v gin,1 (12)
s ) are subject to
(16)
vgin,1
, ug , s
1
are introduced to guarantee the
v gin,1
ug ,2· v gin,2
vgin,2
ug ,1· v gin,2
u g,3·v gin,3
vgin,3
ug ,2· v gin,3
s
vgin =
vgin, k k=1
ug , s 1·v gin, s
(13)
0
The range of vgout can be divided into the same amount of segments as that of vgin :
vgin, s
1
vgin, s
1
ug , s 2·v gin, s
1
ug , s 1·v gin, s
If vgin is in the k -th (1 < k < s ) segment, then: 4
(18)
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
ug ,1, ug ,2, , ug , k 1 = 1 ug , k , u g, k + 1, , u g, s 1 = 0
In this case, the input and output ports of the component should be split into multiple ports. Splitters and concentrators are thus introduced to facilitate the split and merge of branches as shown in Fig. 4. In this paper, the original branches in the EH are named primary branches, and the branches split from primary branches are named secondary branches. If a nonlinear SIMO converter has one input and b outputs, then one splitter, b concentrators and one linear mapping component should be introduced, as shown in Fig. 5. The energy flows in all secondary branches form a vector V (energy flows in secondary branches connected to the input ports of each node come first in V ), and the primary branch vector V should be expanded as:
(19)
Eq. (18) can be written in matrix form to make it easily automated by computers:
V in g Ug , a
V in g
V in g Ug , b
(20)
where in in V in g = diag (v g ,1, v g ,2,
, v gin, s ) T
in in V in g = [ vg,1 vg ,2
vgin, s ]
Ug , a = [ Ug 0 ]T , Ug , b = [1 Ug ]T ug , s 1 ] Ug = [ug,1 ug,2
(21)
C. Formulating port-branch incidence matrices
(2) Type 2: SIMO energy converters
First, consider the simplest case where the EH has only one nonlinear SISO energy converter, which is represented by node g , and its input and output ports are connected to only one branch respectively. Using the piecewise linearization method, node g is converted into an s -input s -output converter, and 2s secondary branches are introduced. A portbranch incidence matrix A g of node g defines the connections between the 2s ports and all the branches (including primary and secondary branches). If the EH has B primary branches, then the port-branch incidence matrix A g has dimensions 2s × (B + 2s ) and is defined as follows:
If the energy output of each output port is proportional to the energy input, e.g., CHP in backpressure operation mode, the energy conversion processes from the input port to each output port can be piecewise linearized separately in a similar way to SISO energy converters. A more complicated situation is when the proportion of energy directed to output ports can be flexibly adjusted, e.g., CHP in extraction condensing operation mode. The input-output characteristics become a multivariate quadratic function:
F = a ·P 2 + b ·Q 2 + c·PQ + d· P + e ·Q + f
(22)
Ag =
where F , P , Q denote the input gas, output electricity and output heat, respectively, and a, b, c, d , e , f are all constants. The characteristics of most SIMO energy converters can be modeled into a multivariate quadratic function, as shown in (22). We conduct a linear mapping to the input-output characteristics:
0s × B Is × s 0s × s 0s × B 0s × s Is × s
(1) If node g represents a nonlinear energy converter, A g is defined as follows:
2
(23)
2
1, branch b is secondary branch and is connected to input
where:
F P Q
1 0 0 = 0 1 c 2a 0 0 1
(26)
where the first B columns correspond to primary branches and the last 2s columns correspond to secondary branches. Then, consider the situation in which the EH has several nonlinear energy converters and several linear energy converters:
F = F1 + F2 F1 = a ·P + d ·P + f1 2 F = b ·Q + e~· Q + f 2
(25)
V =[ V T V T ]T
port k of node g
F P , Q
Ag (k , b) =
c2
b =b 4a ~ e = e cd 2a f = f1 + f2
- 1, branch b is secondary branch and is connected to output port k of node g 0, otherwise
(24)
After this mapping, the energy converter can be equivalently seen as two parallel independent ones. Curves F1 and F2 can be separately piecewise linearized in a similar way to SISO energy converters. A broader sense of “output energy” is introduced, i.e., P = P + (c 2a) Q , which is a linear combination of the original output energy types. Energy conversion with more output ports can be handled similarly.
(27) (2) If node g represents a linear energy converter, A g is defined as follows:
Ag = [ A g 0 ]
(28)
where A g is the port-branch incidence matrix of node g before piecewise linearization and the number of columns of the 0 matrix is equal to the length of V .
B. Introducing Splitter/Concentrator Through piecewise linearization, a nonlinear energy conversion or charging/discharging process in a SISO component is divided into several parallel processes with constant efficiencies. From a component level viewpoint, a nonlinear SISO component with input vgin and output vgout can be equivalent to a linear multi-input multi-output (MIMO) vgin, k (k = 1, 2, , s ) component with inputs and outputs out vg, k (k = 1, 2, , s ) , as shown in Fig. 3.
D. Formulating Converter Characteristic Matrices (1) SISO energy converters Using the method of piecewise linearization, a nonlinear SISO converter is converted into an s -input s -output converter. The converter characteristic matrix Hg has dimensions s × 2s and is defined as Fig. 3. Piecewise linearization of a nonlinear SISO component.
5
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
v
v
in g
in g
Primary Branch
Secondary Branch
Primary Branch
Primary Branch
Nonlinear Component η g (vgin )
Secondary Branch
Linear Component
Splitter
vgout
Concentrator
Primary Branch
vgout
η g ,k (k = 1, 2, , s ) Fig. 4. Splitter and concentrator.
Primary Branch
vgin
Primary Branch
Nonlinear Component
v
Primary Branch
vgout _ 2
Secondary Branch
Secondary Branch
in g
vgout _1
Primary Branch
Concentrator Linear Mapping
Linear Component
Splitter
vgout _1 vgout _ 2
Concentrator
Fig. 5. Piecewise linearization of a nonlinear SIMO energy converter.
follows:
Hg = [
Is× s ]
g, s × s
The converter characteristic matrices of linear energy converters remain the same.
(29)
where g, s × s = diag ( g ,1, g ,2, , g , s ) and conversion efficiency in the k -th segment.
g,k
E. Formulating energy balance matrices
is the constant energy
Given the port-branch incidence matrix A g and the component characteristic matrix Hg , the nodal energy balance matrix for node g can be calculated:
(2) SIMO energy converters Using piecewise linearization, a nonlinear 1-input b -output converter is converted into an s -input b·s -output converter. If the energy output of each output port is proportional to the energy input, the converter characteristic matrix Hg has dimensions b·s × (b + 1)· s and is defined as follows: 1, s × s
Hg =
2, s × s b, s × s
Is × s 0s × s 0s × s I s× s 0s × s
Suppose that EH has N nodes. The system energy balance matrix Z combines the nodal energy balance matrix of all nodes in the EH:
Z = Z1T Z 2T
0s × s 0s × s 0s × s I s×s
T
ZN
T
(33)
The energy balance matrix for each node defines the relationship between the energy flows in the branches connected to the node:
(30)
ZV = 0
If the proportion of energy directed to each output port can be flexibly adjusted, the converter characteristic matrix Hg has dimensions s × (b + 1)·s and is defined as follows:
(34)
F. Formulating splitter/concentrator characteristic matrices Splitter/concentrator characteristic matrix Wg of node g defines the relationship between the primary branches and secondary branches related to node g . First, consider the simplest case where the EH has only one nonlinear SISO energy conversion or storage component, which is
1, input port k is related to energy conversion process i Hg (i , k ) =
(32)
Zg = H g A g
1/ i , output port k is related to energy conversion process i 0, otherwise (31) 6
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
represented by node g , and its input and output ports are connected to only one branch respectively. A g is a matrix whose elements are the absolute value of those in A g (the port-branch incidence matrix of node g before piecewise linearization). If the number of segments in piecewise linearization is s , then:
11× s 01× s 01× s 11× s
A g = Ag 0
V =
(38)
W=
(40)
vgdisch
D g
T
(44)
C
D
(45)
1]
(46)
g
The input incidence matrix X is an m × B (B is the length of V ) matrix that relates the energy inputs to the branch energy flows. Similarly, the output incidence matrix Y is an n × B matrix that relates the energy outputs to the branch energy flows. We define: (47)
X=[ X 0 ], Y=[ Y 0 ]
where the number of columns of 0 matrix equals to the length of V . Then,
(41)
(48)
Vin = XV , Vout = YV
where and denote the charging and discharging efficiencies of component g , respectively. If the charging and discharging efficiencies of energy storage component g vary with its charging and discharging powers respectively, its efficiency curves can be also piecewise linearized in a similar way to the energy converters. In this case, secondary variables vgchar and vgdisch ,k ,k (k = 1, 2, s ) are introduced to represent the amounts of charging and discharging powers that fall in each segment and then Eg is: C g
Eg
H. Comprehensive energy flow equations for EH
The modeling of energy storage components is similar to that of SISO energy converters where the input port models charging and the output port discharging. The only difference is that energy can be stored in energy storage components and state of charge (SOC) is used to describe the ratio of stored energy to storage capacity. Let vgchar and vgdisch denote the charging and discharging powers of energy storage component g respectively, the rate of change of the SOC of this component is: C g
(43)
1
The formulation of the splitter/concentrator characteristic matrices of energy storage components is the same as that of energy converters.
G. Formulating matrices for energy storage components
Eg = vgchar ·
and augmented branch vector V :
0
T
Zg = H g A
(39)
WV = 0
g
where C = [ gC,1, gC,2, , gC, s ] and D = [1 gD,1, 1 gD,2, , 1 gD, s]. The nodal energy balance matrix for energy storage component g can be calculated as follows:
T T WG
V
Hg = [
The splitter/concentrator characteristic matrix for general situations can be formulated in a similar way. Suppose that EH has G nonlinear components in the EH. The system splitter/concentrator characteristic matrix W combines the nodal splitter/concentrator characteristic matrix of all nonlinear components: T W2
D g
Based on the above, the comprehensive energy flow equations for EH are: T
T T [ X T Y T Z T W T ]T V = [ Vin Vout 0T 0T ]
v2 vin ,1
v1
v3
CHP
ηe ηh
(49)
Except for the difference in formulation of the input/output incidence matrix and energy balance matrix, equations (49) have extra variables for energy flows in secondary branches and equations for
vout ,1
Gas
(42)
where A g can be established in a similar way to a SISO energy converter. Considering the energy balance of energy storage component g , i.e., (42), its component characteristic matrix is:
Then,
T W1
k=1
D g, k
D g, k
branch incidence matrix A
(37)
Wg [ V T V T ]T = Wg V = 0
vgdisch ,k
where and are the constant charging and discharging efficiencies of energy storage component g in the k -th segment. To model the SOC for energy storage component g under the proposed standardized matrix modeling framework, we introduce a virtual output port and a virtual branch Eg and construct augmented port-
(36)
A g Wg ]
C g, k
C g, k
where A g V denotes the energy flows in the primary branches related to node g . 11 × s (01 × s ) denotes a 1 × s vector whose elements are all ones (zeros). We define splitter/concentrator characteristic matrix Wg of node g as:
Wg = [
s
vgchar ,k · k=1
(35)
A g V = Wg V
Wg =
s
Eg =
v4
vout ,2
ηc WARG
Fig. 6. A CCHP-based EH. 7
v5
vout ,3
Electricity
Heat
Cooling
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
splitters/concentrators compared with (9). A simple example based on a combined cooling, heating and power (CCHP)-based EH shown in Fig. 6 is given to illustrate the proposed standardized matrix modeling approach of EH with variable efficiencies. The EH consists of a CHP unit and a water absorption refrigerator group (WARG). The input and output vectors of this CCHP-based EH are:
Vin = [vin,1]
(50)
Vout = [ vout ,1 vout ,2 vout ,3 ]T
(51)
and WARG can be calculated:
W2 = [ (52)
v 6 ]T
v5 v 1
are
A1 = [ A1 0 ] =
A2 =
0 1 0
03 × 5 I3 × 3 03 × 3 03 × 5 03 × 3 I3× 3
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
=
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
1 0 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0
e
0 0 0 0 0
H2 = [
c,3 × 3
h
I
]=
3×3
c,1
0 c,2
0 0
0
c,3
1 0 0 0 1 0 0 0 1
Gas
c,3
0 1 0
0 0 1
(59)
0 0 0 0 0 0
1 0
0 1 1 1 0 0 0 1 0 0 0 1 1 1
(60)
(61)
0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
(62)
0 0 1 0 0 1 0 0 0 0 0
0 0 0 0 0 1 0 0 0 1 0
0 0 0 1 0 0 0 0 0 0 1
0 0 0 0 0 0 c,1
0 0 0 0 0 0 0
0 0 1 0
c,2
0 0 0 0 0 0 0 0
0 1 0
1 0
c,3
0 0 0 0 0 0 1 0 0 0 1
0 0 0 0 0 0 0 1 0 0 1
0 0 0 0 0 0 0 0 1 0 1
v1 v2 v3 v4 v5 v1 = v2 v3 v4 v5 v6
vin,1 vout ,1 vout ,2 vout ,3 0 0 0 0 0 0 0
4. Standardized data structure This section presents a standardized data structure of the proposed modeling approach, which supports an automatic procedure for establishing the matrix model of arbitrary EH with variable energy efficiencies and related computations. The proposed data structure consists of two basic tables, i.e., node table and branch table, which act as the standardized input for MES information.
(57)
Given the port-branch incidence matrices and the component characteristic matrices, the nodal energy balance matrix for the CHP
vin ,1 v1
0
1 0 0
Eq. (63) provides a matrix model for EH with variable efficiencies. The proposed matrix modeling approach can be highly automated by computers and facilitates a highly efficient optimization in operation and planning.
(56)
0 0
0 0
(63)
0 0 0 0 0 1
1 0 0 1
A2 W2 ] =
0 1 0 0 1 0 0 0 0 0 0
h
(55)
e
c ,2
The input incidence matrix and output incidence matrix of the EH
1 0 0 0
The component characteristic matrices for the CHP and WARG are:
H1 = H1 =
0
c,1
0 0
(58)
Based on the above, the comprehensive energy flow equations for the EH are:
(54)
0 0 0 0 0 0 1 0 0 1 0 0
0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Y = [ Y 0] =
(53)
0 0 1
0 1
X = [ X 0 ] = [ 1 0 0 0 0 0 0 0 0 0 0]
Assume that the CHP is Node #1 and the WARG is Node #2, the port-branch incidence matrices for the CHP and WARG are:
1 0 0
1 0
The splitter/concentrator characteristic matrix of the WARG is:
The CHP is assumed to operate in backpressure mode (the production of electricity and heat are both proportional to the gas input) and produce electricity and heat with fixed efficiencies e and h , respectively. The WARG is assumed to convert its heat input into cooling with a variable efficiency c (v4) . Using the piecewise linearization method described above, the domain of function c (v4) can be divided into 3 segments (taking 3 segments as an example here) with fixed efficiencies c,1 , c,2 , and c,3 . As a result, the SISO WARG is converted to a 3-input 3output converter, and one splitter, one concentrator and 6 secondary branches vk (k = 1, 2, , 6) should be introduced. The EH after piecewise linearization is shown in Fig. 7. The branch vector V should be expanded:
V = [ V V ]T = [ v1
h
Z2 = H2 A2 =
The set of branches (energy flows) can be written as:
V = [ v1 v2 v3 v4 v5 ]T
e
Z1 = H1A1 =
(1) Node table
vout ,1
v2 CHP
v3
ηe ηh
vout ,2
v4 v1′ ~v3′
v4′ ~v6′
WARG
Splitter
Concentrator
ηc ,1 , ηc,2 , ηc,3 Fig. 7. A CCHP-based EH after piecewise linearization. 8
v5
vout ,3
Electricity
Heat
Cooling
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
In the node table, each row describes a single energy conversion/ storage component or an input/output of the EH, including:
Table 3 Parameters of the energy components in Fig. 8.
• Node No.: Numbering all nodes in EH, including energy conversion/storage components and inputs/outputs of EH, in order and starting from 1. • Node type: Identifying the type and operating mode of each node. • • •
1 2
3 1
3 4 5 6
−1 0 0 0
e A c
0 0 0 0
h B c
0 0 0 0
0 C c
0 0 0 0
Node Capacity
Segment Number
CCHP CWARG
0 3
0 0 0 0
0 0 0 0
Branch Type
Source Node No.
Sink Node No.
1 2 3 4 5
4 1 3 3 2
3 1 1 1 2
1 4 5 2 6
400
Out = 0.00003041· In3 + 0.01901·In2 + 0.2593· In * Out = 3·In
400 El: 300 Therm:420 900 800 (3.2 MWh)
A 2 B c = c ·v4 + c · v4 + WARG, respectively.
El: Out = 0.0001150·In2 + 0.2305·In
Therm: Out = 0.0001611·In2 + 0.3228·In ** Out = 0.8·In Charging Efficiency = 0.00005 In + 0.93 Discharging Efficiency = 0.00005 Out + 0.93
C c ).
CCHP and CWARG are capacities of CHP and
(2) Branch table In the branch table, each row describes one energy flow (only for primary branches) in the EH, including:
• Branch No.: Numbering all energy flows in the EH in order and starting from 1. • Branch Type: Identifying the type of each energy flow. For example, • •
the electricity, cooling, heat and gas flows are numbered as 1, 2, 3 and 4, respectively. Source Node No.: Node No. of the source node of the branch. Sink Node No.: Node No. of the sink node of the branch.
Taking the EH shown in Fig. 6 as example, the branch table is shown as Table 2. The component characteristic matrix Hg for each energy conversion or storage component is established based on the node table. The portbranch incidence matrix A g for each component, input incidence matrix X and output incidence matrix Y are established based on the branch table. Splitter/concentrator characteristic matrix Wg is established based on both the node table and branch table.
Table 2 Branch table of the EH in Fig. 6. Branch No.
CERG
* Part-load performance of the existing chiller with specific settings of condensing pressure and outdoor temperature; ** Part-load performance of natural gas-powered microturbine with a heat exchanger.
Table 1 Node table of the EH in Fig. 6. Energy Efficiency Parameter
Performance [5,26,31,33]
AB HS
The inputs/outputs of EH do not have data of Energy Efficiency Parameter, Node Capacity and Segment Number. Taking the EH shown in Fig. 6 as example, the node table is shown as Table 1. The energy conversion characteristics of CHP and WARG are described as v2 = e ·v1, (v3 + v4) = h · v1, v5 = cA ·v43 + cB ·v42 + cC · v4 (i.e.,
Node Type
Capacity (kW)
HP CHP
For example, −1 for input, 0 for output, 1 for WARG, 2 for CHP unit in extraction condensing operation mode, 3 for CHP unit in back pressure operation mode and so on. A component library can be constructed for index. Energy Efficiency Parameter: Providing information on energy conversion characteristics of energy conversion components and energy charging/discharging characteristics of energy storage components. The number of parameters depends on the type of node. Node Capacity: Providing information on capacity of each energy conversion or storage component, which is used when modeling the operating constraints of the EH. Segment Number: Determining the number of segments for piecewise linearization. 0 is for energy conversion and storage components with constant efficiencies.
Node No.
Converter
5. Case study This section illustrates the application of the proposed standardized matrix modeling approach in optimal operation of EH. A. System description Fig. 8 illustrates an EH which consists of a CHP in backpressure
Fig. 8. An EH with three nonlinear energy components. 9
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
operation mode (heat-to-power ratio is a constant), a compression electric refrigerator group (CERG), an electric heat pump (HP), an auxiliary boiler (AB) and a heat storage component (HS). The energy inputs are electricity and gas, while the outputs consist of cooling, electricity and heat. The CHP, CERG and HS are considered to exhibit considerable efficiency variations with their input or output powers changing. The electric efficiency of the CHP in partial loads can degrade more than 20% compared to the full-load efficiency [30]. The specific type of CERG in Fig. 8 is an air-cooled centrifugal chiller whose highest coefficient of performance (COP) occurs at a part-load ratio of 0.71–0.84 rather than at full load [31]. The efficiencies of the HP and AB do not differ significantly under part-load conditions [26,32]; therefore, their efficiencies are assumed to be constant. Table 3 lists the parameters of the energy components of this system. The hourly electricity, heat and cooling demand patterns refer to a hospital site in a Mediterranean area [5] and are indicated in Fig. 9. The gas price is set at 20 euro/MWh and is considered constant during the daily analysis period [5]. The variation in hourly electricity prices for the selected day corresponding to a real case in Italy [5] is also shown in Fig. 9. Using the proposed modeling approach described above, the piecewise linearized EH with splitters/concentrators is shown in Fig. 10 (taking 2 segments as an example here). Given the topology of EH, parameters of energy components and the pre-defined number of segments for piecewise linearization (taking 2 segments as an example here), the computer can automatically formulate the comprehensive energy flow equations of EH as shown in (64).
given time horizon. By replacing the comprehensive energy flow equations (9) with (49) and introducing piecewise linearization constraints (e.g., (20)–(21)), a revised optimal operation model based on the model in [22] is used. The decision variables of this problem are the continuous variables V , the continuous variables Vin , the continuous variables E (state of charge of energy storage components) and the binary variables Ug in each time period, making the problem a mixed integer linear programming (MILP) problem. C. Optimization results and sensitivity analysis The objective of the EH operational optimization in this case study is to minimize the daily operating cost of the whole system. Three cases are considered: (1) Case 1: The exact performances of CHP, CERG and HS are considered, making the optimal operation problem a NP problem; (2) Case 2: The efficiencies of the CHP, CERG and HS are approximated as constants (e.g., efficiencies of converters are calculated by dividing maximum output by maximum input), making the optimal operation problem a linear programming problem; (3) Case 3: The proposed standardized matrix modeling approach is used. The CHP, CERG and HS are piecewise linearized, making the optimal operation problem an MILP problem. In addition, subcases of Case 3 with different numbers of segments in piecewise linearization are considered. When the number of segments is large en-
(64)
B. Multi-period optimal operation model for EH
ough (e.g., 300), the optimization results are considered accurate and taken as the reference compared with the other optimization results.
The optimal energy flow problem of a single EH determines how energy flows should be dispatched among different components to minimize the cost or maximize the profit of the whole system over a
Case studies are executed in MATLAB with YALMIP [34] on a PC 10
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
Fig. 9. Hourly load patterns and electricity prices of the selected day.
with an Intel Core i7-7500U 2.70 GHz CPU and 16 GB RAM. Cases 2 and 3 are solved using GUROBI. Case 1 is solved using fmincon with the sequential quadratic programming (SQP) method. The maximum number of iterations is set at 3000. The optimization results of Case 2 are used as initial values for Case 1. The optimal operating costs, computation times and relative errors of optimal values for these three cases are summarized in Table 4 and Fig. 11, in which s denotes the number of segments in piecewise linearization. For Case 2, the relative error of the optimal operating cost is over 13%, which indicates that the constant efficiency approximation does not adequately model the EH. With an increase in the number of nonlinear components, nonlinearity of components and complexity of EH, the constant efficiency model of EH results in larger errors due to approximations, making its optimization results impractical. With the proposed standardized matrix modeling approach, the relative error of Case 3 is below 1.0% and 0.1% when s is equal to 12 and 36, respectively. With no evident calculation time increase (less than 0.2 s) compared to Case 2, the proposed approach reduces the approximation error from 13.7% to 0.1% in Case 3. Using the results provided by Case 2 as initial values, Case 1 obtains much more accurate results than Case 2. Nevertheless, the NP problem cannot guarantee a global optimal solution, and its solution is time-consuming. The proposed standardized
v2
Electricity
Gas
Splitter
vin ,1
v1′, v2′
CERG
v7′ , v8′
Table 4 Optimal operating costs, computation times and relative errors for the three given cases. Case Case Case Case Case Case Case Case Case Case Case Case Case Case Case Case
(s (s (s (s (s (s (s (s (s (s (s (s (s
= = = = = = = = = = = = =
2) 4) 6) 8) 10 ) 12 ) 20 ) 36 ) 50 ) 76 ) 100 ) 200 ) 300 )
Computation Time (s)
Relative Error* (%)
313.628 264.440 280.842 289.797 297.182 301.089 303.020 304.222 305.683 306.230 306.325 306.372 306.393 306.416 306.420
292.73 0.43 0.47 0.61 0.60 0.50 0.50 0.49 0.53 0.58 0.67 0.88 1.29 2.73 3.39
2.35 13.70 8.35 5.43 3.01 1.74 1.11 0.72 0.24 0.06 0.03 0.02 0.01 0.001 –
* Relative errors for each case/subcase are calculated according to the optimization results of Case 3 (s = 300 ).
Concentrator
vout ,1
v6
v1 v3
vin ,2
1 2 3 3 3 3 3 3 3 3 3 3 3 3 3
Optimal Value (Euro)
v4
HP
′ v9′ , v10
Concentrator
CHP
Splitter
v3′ , v4′
′ , v12 ′ v11
Concentrator
v11 v7
v8
Splitter
v5′ , v6′
AB
v10
Fig. 10. The piecewise linearized EH with splitters/concentrators. 11
HS
Electricity
′ , v14 ′ v13 Concentrator
v9
v5
vout ,2
Cooling
v12
vout ,3
Heat
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
Fig. 11. Computation times and relative errors for different numbers of segments in Case 3.
matrix modeling approach makes the optimal operation problem an MILP problem that can be efficiently solved by many existing commercial software (e.g., GUROBI, CPLEX, MOSEK). It should also be noted that the inefficiency of the operation is better approximated when the number of segments in piecewise linearization increases, which leads to a higher operating cost. As shown in Fig. 11, relative error decreases rapidly when s increases from 2 to 20, and computation time increases relatively slowly when s increases from 30 to 70; consequently, the optimal number of s is between 30 and 70 considering the balance between nonlinear approximation accuracy and calculation efficiency. Figs. 12 and 13 show the operating states of each energy component in the EH over a 24-hour period obtained in Case 1, Case 2 and Case 3
(s = 100 ). The electricity and heat demands are provided by the CHP, HS and AB during the daytime because the price of electricity is high. Load demands and electricity price decrease during the night; therefore, it is inefficient for the CHP to work at a low load level, and it is more cost-effective to purchase electricity from the distribution system to meet the electricity demand and provide power for the HP to produce heat. HS contributes to reducing the heating cost by storing heat when electricity price is low during the night and discharging it in the daytime. The AB is used as a supplemental source for heat supply. The constant efficiency approximation in Case 2 leads to a wrong and uneconomical operation strategy as shown in Fig. 12: (1) The CHP still works at a relatively low load level during the night without regard to its part-load performance. (2) The electricity input of the CERG
CERG: Case 3 HP: Case 3 CHP: Case 3 AB: Case 3 HS: Case 3
CERG: Case 2 HP: Case 2 CHP: Case 2 AB: Case 2 HS: Case 2
Fig. 12. Operating states of the energy components in the EH of Case 2 and Case 3 (s = 100 ). 12
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
CERG: Case 3 HP: Case 3 CHP: Case 3 AB: Case 3 HS: Case 3
CERG: Case 1 HP: Case 1 CHP: Case 1 AB: Case 1 HS: Case 1
Fig. 13. Operating states of the energy components in the EH of Case 1 and Case 3 (s = 100 ).
decreases because the efficiency is supposed to remain high at any load level. (3) The HS discharges heat with higher power between 10:00 and 14:00 without considering the increasing power loss with increasing discharging power. The NP problem of Case 1 results in a sub-optimal operation strategy by obtaining the local optimal solution, which is highly impacted by the initial values of decision variables. Comparing Fig. 12 and Fig. 13, it can be found that the operation strategy obtained by Case 1 is very similar to that of Case 2, which is because the optimization results of Case 2 are used as initial values for Case 1. Though Case 1 considers the exact performances of CHP, CERG and HS, its optimization results are misled by Case 2 which ignores the part-load performances. In addition, between 8:00 and 22:00, the input powers of CHP and CERG in Case 1 are almost the same as those in Case 3 (s = 100 ), but different from those in Case 2. This indicates that the part-load performances of components are well modeled by the proposed approach used in Case 3 (s = 100 ).
condition for the operation model and it can be obtained in advance by the weather forecast. Different efficiency curves can be set in advance for different temperature differences during different time periods in the optimal operation problem. In other words, the efficiency variation caused by external factor, e.g., environmental factor, can be obtained as conditional efficiency in advance based on the forecast of the factor. In addition, for the optimal operation problem, the external factor does not depend on the operation strategy. Therefore, the operation model under a certain external factor, e.g., temperature difference, still has a linear formulation. By contrary, the efficiency variation and nonlinearity caused by partial load is the key difficulty in the optimal operation problem, which is handled using the proposed modeling approach. The piecewise linearization process can be further improved using self-adaptive technique, i.e., the widths and number of piecewise linearization segments can be unevenly chosen according to the nonlinearity of efficiency curves. For example, the widths of different segments for one component can be different. Narrower segments are used to approximate the part with higher nonlinearity on the curve, which further improves the approximation accuracy. In addition, different energy conversion and storage components have different nonlinearities and they can use different numbers of segments in piecewise linearization to reduce the number of segment variables without sacrificing modeling accuracy. The proposed approach can effectively model efficiency curves with all kinds of shapes, irrespective of whether they are convex or nonconvex. E.g., linear, quadratic, cubic, higher order polynomial and exponential efficiency curves can all be easily handled by the proposed piecewise linearization based modeling approach.
6. Discussion The proposed matrix modeling approach of EH can be effectively applied to a more complex case because the proposed approach can be highly automated by computers. In practical MES, more linear and nonlinear components, more complicated system structure and a longer time horizon would lead to an increase in the size of the matrices and the number of continuous and binary variables in the operation model. The matrices obtained, though very large in scale, can be easily handled by computers using effective techniques such as sparse technique. Binary variables in the proposed approach are tightly constrained by (20)–(21). Once the value of one binary variable for one segment is determined, the values of the other binary variables for the other segments can be immediately determined. So the number of binary variables does not significantly impact the calculation time which is proved in the case study (e.g., Fig. 11). The energy efficiencies of components may depend on several parameters besides input or output power, e.g., temperature difference; however, it does not introduce nonlinearity to the operation model. For example, the efficiency of heat pump varies with the temperature difference, e.g., the heating efficiency is high when temperature difference is small during the daytime while it is low when temperature difference is large at night. However, the outdoor temperature is a boundary
7. Conclusion The variable energy efficiencies of the energy conversion and storage components in MES complicate the modeling of EH in two ways: (1) The nonlinearity caused by variable efficiencies complicates the analysis and optimization of EH; (2) The nonlinear energy components can hardly be automatically modeled by computers. To jointly address these issues, we propose a standardized matrix modeling approach to automatically model EH with variable energy efficiencies. Using piecewise linearization, a nonlinear energy component is converted into a “splitter, linear energy component, concentrator” three-layer structure. 13
Electrical Power and Energy Systems 119 (2020) 105876
W. Huang, et al.
The nonlinear energy conversion and storage relationship in EH can thus be further modeled under a linear modeling framework using matrices. Such modeling approach facilitates an automatic modeling of an arbitrary EH with various kinds of nonlinear energy components. Using the proposed modeling approach, the optimal operation of a fivecomponent EH with energy storage shows a great improvement on the accuracy and calculation efficiency compared to the exiting constant efficiency model and nonlinear model, respectively. Future work includes extending the modeling framework to consider the dynamics of energy conversion in different time scales.
[8] Tian X, Zhao R. Energy network flow model and optimization based on energy hub for big harbor industrial park. J Coastal Res 2015;73:298–303. [9] Krause T, Kienzle F, Liu Y, Andersson G. Modeling interconnected national energy systems using an energy hub approach. In: Proc IEEE powertech, 2011, Trondheim, Norway; Jun. 2011. p. 1–7. [10] Guler B, Çelebi E, Nathwani J. A ‘Regional Energy Hub’for achieving a low-carbon energy transition. Energy Policy 2018;113:376–85. [11] Mohammadi M, Noorollahi Y, Mohammadi-ivatloo B, Hosseinzadeh M, Yousefi H, Khorasani ST. Optimal management of energy hubs and smart energy hubs–a review. Renew Sustain Energy Rev 2018;89:33–50. [12] Wang Y, Zhang N, Zhuo Z, Kang C, Kirschen DS. Mixed-integer linear programmingbased optimal configuration planning for energy hub: Starting from scratch. Appl Energy 2018;210:1141–50. [13] Sadeghi H, Rashidinejad M, Moeini-Aghtaie M, Abdollahi A. The energy hub: an extensive survey on the state-of-the-art. Appl Therm Eng 2019:114071. [14] Geidl M, Andersson G. Optimal power flow of multiple energy carriers. IEEE Trans Power Syst 2007;22(1):145–55. [15] Biglia A, Caredda FV, Fabrizio E, Filippi M, Mandas N. Technical-economic feasibility of CHP systems in large hospitals through the Energy Hub method: the case of Cagliari AOB. Energy Build 2017;147:101–12. [16] Cheng Y, Zhang N, Wang Y, Yang J, Kang C, Xia Q. Modeling carbon emission flow in multiple energy systems. IEEE Trans Smart Grid July 2019;10(4):3562–74. [17] Waibel C, Evins R, Carmeliet J. Co-simulation and optimization of building geometry and multi-energy systems: Interdependencies in energy supply, energy demand and solar potentials. Appl Energy 2019;242:1661–82. [18] Geidl M, Andersson G. Optimal coupling of energy infrastructures. Proc IEEE Lausanne Power Tech 2007:1398–403. [19] Almassalkhi M, Hiskens I. Optimization framework for the analysis of large-scale networks of energy hubs. In: Proc 17th power syst comput conf; 2011. p. 1–7. [20] Beccuti G, Demiray T, Batic M, Tomasevic N, Vranes S. Energy hub modelling and optimisation: an analytical case-study. In: PowerTech 2015 IEEE Eindhoven; 2015. p. 1–6. [21] Wang Y, Cheng J, Zhang N, Kang C. Automatic and linearized modeling of energy hub and its flexibility analysis. Appl Energy Feb. 2018;211:705–14. [22] Wang Y, Zhang N, Kang C, Kirschen DS, Yang J, Xia Q. Standardized matrix modeling of multiple energy systems. IEEE Trans Smart Grid 2019;10(1):257–70. [23] Fabrizio E, Filippi M, Virgone J. An hourly modelling framework for the assessment of energy sources exploitation and energy converters selection and sizing in buildings. Energy Build 2009;41(10):1037–50. [24] Evins R, Orehounig K, Dorer V, Carmeliet J. New formulations of the ‘energy hub’ model to address operational constraints. Energy 2014;73:387–98. [25] Sakti A, Gallagher KG, Sepulveda N, Uckun C, Vergara C, de Sisternes FJ, et al. Enhanced representations of lithium-ion batteries in power systems models and their effect on the valuation of energy arbitrage applications. J Power Sources 2017;342:279–91. [26] Milan C, Stadler M, Cardoso G, Mashayekh S. Modeling of non-linear CHP efficiency curves in distributed energy systems. Appl Energy 2015;148:334–47. [27] Almassalkhi MR, Towle A. Enabling city-scale multi-energy optimal dispatch with energy hubs. In: Presented at 2016 power syst comput conf (PSCC), Genoa; 2016. p. 1–7. [28] Holjevac N, Capuder T, Zhang N, Kuzle I, Kang C. Corrective receding horizon scheduling of flexible distributed multi-energy microgrids. Appl Energy 2017;207:176–94. [29] Shabanpour-Haghighi A, Seifi AR. Energy flow optimization in multi-carrier systems. IEEE Trans. Ind. Informat. 2015;11(5):1067–77. [30] Bianchi M, Pascale AD, Melino F, Peretto A. Performance prediction of micro-CHP systems using simple virtual operating cycles. Appl Therm Eng 2014;71(2):771–9. [31] Yu FW, Chan KT. Part load performance of air-cooled centrifugal chillers with variable speed condenser fan control. Build Environ 2007;42(11):3816–29. [32] Anonymous. Boiler efficiency guide: facts about firetube boilers and boiler efficiency. CleaverBroocks; Mar. 2010. [Online]. Available: http://cleaverbrooks.com/ reference-center/insights/Boiler%20Efficiency%20Guide.pdf. [33] Zhou Z, Liu P, Li Z, Ni W. An engineering approach to the optimal design of distributed energy systems in China. Appl Therm Eng 2013;53:387–96. [34] Lofberg J. YALMIP: a toolbox for modeling and optimization in MATLAB. In: presented at proc IEEE int'l symp comput aided control syst des, Taipei; 2004. p. 284–289.
CRediT authorship contribution statement Wujing Huang: Methodology, Software, Validation, Investigation, Data curation, Writing - original draft, Visualization. Ning Zhang: Conceptualization, Writing - review & editing, Supervision. Yi Wang: Writing - review & editing. Tomislav Capuder: Writing - review & editing. Igor Kuzle: Writing - review & editing. Chongqing Kang: Writing - review & editing. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgement This work was supported in part by China-Croatia InterGovernmental S&T Cooperation Project Flexible Integrated Demand Response in Multi-Energy System (FIRE), National Natural Science Foundation of China (No. 71961137004), Tsinghua University Initiative Scientific Research Program (20193080026), and EU-CHINA Research and Innovation Partnership Project Instigation of Research and Innovation Partnership on Renewable Energy, Energy Efficiency and Sustainable Energy Solutions for Cities (IRES-8), contract number ICI+/2014/347-910. References [1] Mancarella P. MES (multi-energy systems): an overview of concepts and evaluation models. Energy 2014;65:1–17. [2] Holjevac N, Capuder T, Kuzle I. Adaptive control for evaluation of flexibility benefits in microgrid systems. Energy 2015;92(4):487–504. [3] Ceseña EAM, Capuder T, Mancarella P. Flexible distributed multienergy generation system expansion planning under uncertainty. IEEE Trans Smart Grid 2016;7(1):348–57. [4] Geidl M, Koeppel G, Favre-Perrod P, Klckl B, Andersson G, Frhlich K. Energy hubs for the future. IEEE Power Energ Mag 2007;5(1):24–30. [5] Chicco G, Mancarella P. Matrix modelling of small-scale trigeneration systems and application to operational optimization. Energy Mar. 2009;34(3):261–73. [6] Rastegar M, Fotuhi-Firuzabad M, Lehtonen M. Home load management in a residential energy hub. Elect Power Syst Res 2015;119:322–8. [7] Bahrami S, Safe F. A financial approach to evaluate an optimized combined cooling heat and power system. Energy Power Eng J 2013;5(5):352–62.
14