Journal of Modern Power Systems and Clean Energy

ISSN 2196-5625 CN 32-1884/TK

网刊加载中。。。

使用Chrome浏览器效果最佳,继续浏览,你可能不会看到最佳的展示效果,

确定继续浏览么?

复制成功,请在其他浏览器进行阅读

A Hybrid State Estimation Approach for Integrated Heat and Electricity Networks Considering Time-scale Characteristics  PDF

  • Tongtian Sheng
  • Guanxiong Yin
  • Qinglai Guo
  • Hongbin Sun
  • Zhaoguang Pan
Department of Electrical Engineering, Tsinghua University, Beijing 100084, China

Updated:2020-07-16

DOI:10.35833/MPCE.2019.000230

  • Full Text
  • Figs & Tabs
  • References
  • Authors
  • About
OUTLINE

Abstract

State estimation (SE) usually serves as the basic function of the energy management system (EMS). In this paper, the time-scale characteristics of the integrated heat and electricity networks are studied and an SE model is established. Then, a two-stage iterative algorithm is proposed to estimate the time delay of heat power transportation in the pipeline. Meanwhile, to accommodate the measuring resolutions of the integrated network, a hybrid SE approach is developed based on the two-stage iterative algorithm. Results show that, in both steady and dynamic processes, the two-stage estimator has good accuracy and convergence. The hybrid estimator has good performance on tracking the variation of the states in the heating network, even when the available measurements are limited.

I. Introduction

ENERGY internet, also named the multi-energy system or the integrated energy system (IES), is proposed to fully release the flexibility of various energy networks such as heating, cooling, natural gas, and even transportation networks, with its advantages in improving energy efficiency, reducing carbon emissions and increasing renewable energy penetration [

1]-[3]. Among various energy networks, the integrated heat and electricity network proves to be promising [4]-[8]. Reference [4] considered the thermal inertia of the heating network in improving the wind power integration. Reference [5] studied the economic dispatching scheme considering the time-scale characteristics of the networks and buildings. Reference [6] proposed a decentralized solution of the coordinated dispatch of electric power and district heating (DH) networks. As the coupling between the electricity network and the DH network became tighter, the marketing model [7] and demand-side management model [8] had been further studied, to make the best of the potential flexibility.

Reference [

9] proposed the architecture of the integrated energy management system (IEMS), which contains the functions including modeling and assessing IES, and dispatch of the multiple energy flows. Reference [9] presented that the online dispatch and control strategies of IES rely on the overall perception and reliable assessment of the integrated system. Reference [10] studied the combined analysis method of integrated heat and electricity networks (IHENs), in which the grid-connected mode and the islanded mode were considered. However, the method was developed on the steady-state flow model. Reference [11] studied the interactions of distributed heat and electricity networks considering the discrepancy of time-scale characteristics. However, the dynamics of the network were simplified and only the quasi-steady energy flow was calculated. In analyzing the dynamics of the DH network, simulation methods were usually used [12]-[15]. The dynamic models of the DH network included transient models [12] and pseudo-transient models [13]-[15]. However, the simulation methods were not practical for real-time monitoring and analysis.

Electricity networks are monitored by remote terminal units (RTUs), phasor measurement units (PMUs), and other intelligent electronic devices (IEDs), the automation level of which is high. However, the monitoring and communication infrastructures of the DH network are not so advanced as the electricity networks due to technical and economic reasons. In recent years, more attention has been paid to improve the automation level of the DH network. In [

16], the concept of the smart thermal grid was proposed which claimed that more information technologies are necessary to be integrated. According to the directive report of European Union, independent heat metering devices had been utilized to obtain the consumption of both heating and domestic hot water in multi-apartment buildings by the end of 2016 [17]. Owing to the heating market reform in China, more automated calorimeters were also equipped in many DH systems [18]. Moreover, supervisory control and data acquisition (SCADA) systems had already been installed in some advanced DH systems to guarantee proper daily operations as well [19].

However, errors caused by the measurements are inevitable. To get a complete description of the integrated system state and detect the bad data, it is necessary to make state estimation (SE). SE has been widely used in electric systems [

20], [23] and water distribution systems [22]-[24]. As for DH networks, SE is usually used to evaluate heat loss and improve energy efficiency [25], [26]. As for the IHEN, [27] proposed an SE approach for the IHEN and proved that the combined SE approach helped to improve the accuracy with higher redundancy of the measurement data. However, the steady model of the DH network was adopted and the dynamics were neglected. In [28] and [29], two novel SE approaches were proposed with the dynamics in pipelines. However, since the measurement data of the DH network were still used for billing purpose in most circumstances [26], the sampling resolutions of the DH network (minute-hour) were relatively lower than those of the electricity network (second-minute), which was not considered in the previous studies.

In this paper, a two-stage SE (TSE) approach and a hybrid estimation algorithm for IHEN is proposed, which considers both the dynamics of the pipeline and different sampling resolutions of two systems.

The remainder of this paper is organized as follows: Section II introduces the model of IHEN. Section III describes the SE model, the TSE approach and the hybrid estimation algorithm. Section IV introduces and discusses the results of numerical simulation and case study. Section V presents the conclusions and directions for future work.

II. Model Description

The model of IHEN consists of an electricity network model, a heating network model, and coupling component models. The heating network includes the hydraulic network and thermal network. In the integrated network, three types of energy are considered: the electric energy in the electricity network, the pressure energy in the hydraulic network, and the thermal energy in the thermal network.

The spread speed of thermal energy in the thermal network is much slower than that of electric energy in the electricity network and the pressure energy in the hydraulic network [

28], which may cause a large time delay in the spread of nodal temperature in IHEN. This delay is mainly caused by the flow of mass through long pipelines. Therefore, in this paper, we utilize a partial-dynamic model to represent the thermal network, in which the dynamic characteristic of the pipeline is considered, while the steady-state models of the electricity network, the hydraulic network and the rest of the thermal network (including heat sources and heat consumers) are retained.

A. Electricity Network

The electricity network model is established based on the Kirchhoff’s current law (KCL), Kirchhoff’s voltage law (KVL), and Ohm’s law. In this paper, the AC power flow model for transmission network [

22] is used. One slack bus bar is needed and other bus bars are given as PQ or photovoltaic (PV) bus bars. The electro-mechanical transient model and the electromagnetic transient model are not considered.

The steady-state power flow equations of electricity network reflect the relationship between electric power injection and bus voltage, which could be represented with the following algebraic equations:

Pi=UijiUj(Gijcos(θi-θj)+Bijsin(θi-θj)) (1)
Qi=UijiUj(Gijsin(θi-θj)-Bijcos(θi-θj)) (2)

where Pi, Qi are the active and reactive power injections of bus i, respectively; Ui, θi are the voltage magnitude and phase angle at bus i, respectively; and Gij+jBij is the ijth element of the complex bus admittance matrix.

B. DH Network

A typical DH network is a dual-pipe system where heat is delivered to customers as heated water. The supply pipe takes the hot water to customers and the return pipe brings the cooled water back to the heat source. The water is circulated with pumps which are usually located at the heat source.

When analyzing a DH network, the entire network is usually separated into a hydraulic network and a thermal network [

10]. The mathematical model is set up based on the following assumptions:

1) Water flow in the pipeline is one-dimensional, incompressible and fully developed.

2) The properties of water such as density, kinematic viscosity, specific heat and so on do not change with the pressure and temperature.

3) Convective heat transfer from the fluid to the surroundings is in the radial direction only.

1) Hydraulic Network

The relationship between the head pressure and mass flow rate is considered in the model of the hydraulic network. Mass flowing through pipe and heat exchanger will cause pressure losses, which could be written as [

13]:

Δppi,i=Kpi,impi,i2i=1,2,...,Npi (3)
Δphex,i=Khex,imhex,i2i=1,2,...,Nhex (4)

where Δppi,i, Δphex,i are the pressure losses of the pipe and heat exchanger, respectively; mpi,i, mhex,i are the mass flow rates of the pipe and heat exchanger, respectively; Kpi,i, Khex,i are the relative resistance coefficients of the pipe and heat exchanger, respectively; and Npi, Nhex are the numbers of the pipes and heat exchangers, respectively.

As for the circulation pump, variable speed pumps are considered in this paper, the relationship between the mass flow rate and head pressure generated by the pump can be expressed as [

30]:

Δppu,i=Kpu,impu,i2i=1,2,...,Npu (5)

where Δppu,i is the pressure generation of the circulation pump; mpu,i is the mass flow rate of the pump; Npu is the number of the pumps; and Kpu,i is the relative coefficient.

Equations(3),(5) can be expressed with a unified matrix form:

Δp=Km2 (6)

where Δp is the pressure loss vector; K is the relative resistance coefficient vector; and m is the mass flow rate vector.

Considering the hydraulic network operates in the steady state, the conservation of mass and pressure energy is also satisfied, respectively [

10], which is described by the flow continuity equation in (7) and the loop pressure equation in (8).

Am=mq (7)
B(Δp)=0 (8)

where A is the node-branch incidence matrix; B is the loop-branch incidence matrix of the DH network; and mq is the nodal flow injection vector.

Note that in [

10], MQ stands for the mass flow rate through each node injected from a supply or discharged to a load. However, in this paper, a heat supply (including adjacent circulation pump) or a heat load is also modeled as a branch in the DH network. With this model, the entire DH network becomes a closed-loop network. If the leak of the water is neglected, we have:

mq=0 (9)

After mq is specified, the hydraulic model could be performed independently [

10]. The mass flow rate m of each pipe is determined by (6)-(8). The nonlinear algebraic equations can be solved using the Hardy-Cross method or Newton-Raphson method. Each loop is considered independently in the Hardy-Cross method and all loops are considered simultaneously in the Newton-Raphson method [31].

2) Thermal Network

The model of the thermal network describes the temperature distribution at each node and thermal power flow in the network. In the thermal network, the components that could cause the change of temperature, e.g., the mixer, the heat exchanger at the source or the loads, and the pipeline, are considered, while circulation pumps and valves are assumed to be lossless of heat. Equation (10) describes the model of a mixer, which serves as a node in the thermal network. The process of fluid mixing is assumed to be instant and homogeneous.

mnd,ioutTnd,iout=(mnd,iinTnd,iin)i=1,2,...,Nnd (10)

where Tnd,iout is the mixture temperature of a node; mnd,iout is the mass flow rate within a pipe leaving the node; Tnd,iin is the temperature of the flow at the end of an incoming pipe; mnd,iin is the mass flow rate within a pipe coming into the node; and Nnd is the number of the nodes.

Equations(11) and (12) represent the thermal model of the heat sources and the loads, respectively, and the dynamics are not considered.

Φsrc,i=Cpmsrc,i(Tsrc,iout-Tsrc,iin)i=1,2,...,Nsrc (11)
Φld,i=Cpmld,i(Tld,iin-Tld,iout)i=1,2,...,Nld (12)

where Φsrc,i, msrc,i, Tsrc,iin, and Tsrc,iout are the heat power, mass flow rate, inlet temperature, and outlet temperature at sources, respectively; Φld,i, mld,i, Tld,iin, and Tld,iout are the heat power, mass flow rate, inlet temperature, and outlet temperature at loads, respectively; Cp is the specific heat capacity of the water; and Nsrc, Nld are the numbers of the heat sources and the heat loads, respectively.

The dynamics of the pipeline are introduced by using a pseudo-transient model [

14], [15]. In this model, the steady-state lumped mass model is combined with a variable transport delay model to enable the progress of the fluid to be tracked in time as the flow varies. Based on the energy balance across a differential pipe segment, one-dimensional partial differential equation (PDE) model is established as [32]:

Tt+mρSTx+λCpρS(T-Ta)=0 (13)

where S is the cross-section area of the pipeline; m is the mass of flow; T is the temperature; λ is the overall heat transfer coefficient per unit length of the pipeline; ρ is the density of water; Ta is the ambient temperature; and t, x are the time point and space point, respectively.

The characteristic lines of (13) are defined by:

dxdt=m(t)ρS (14)

Considering the temperature along the characteristic line, the PDE along the line can be written as (15) according to (13) and (14).

dTdt=Tt+mρSTx=λCpρS(Ta-T) (15)

Figure 1 shows the scheme of the characteristic line in the pipe i, where Lpi,i is the total length; Tpi,iout(t) is the outlet temperature at the time t; and Tpi,iin(t-spi,i) is the inlet temperature at the time t-spi,i, and spi,i is the time delay of mass flow from the inlet of the pipeline to the outlet along the characteristic line, which is calculated by (16) according to (14).

t-spi,itmpi,i(τ)dτ=ρSpi,iLpi,ii=1,2,...,Npi (16)

Fig. 1 Scheme of characteristic line of (13).

where mpi,i is the mass flow rate of the pipe; and Spi,i is the cross-section area of the pipe.

Thus the integral equation along the characteristic line of pipe i is calculated by (15) and (17), and the relationship of the inlet and outlet temperature of the pipe is represented by (18).

Tpi,iint-sTpi,iouttdTT-Ta=-t-stλpi,iCpρSpi,idt (17)
Tpi,iout(t)=(Tpi,iin(t-spi,i)-Ta)e-λpi,iCpρSpi,ispi,i+Ta (18)

However, since the measurements have specific sampling resolution, and the SE only considers the state at sampling points, the integral equation (16) cannot be applied directly and the discretization is needed. Given the current time point t and the measuring period Δt, there is an integer γ satisfying (19) for each pipe:

k=0γ-1mt-kρSΔtLk=0γmt-kρSΔt (19)

where mt-k is the mass flow rate at time t-kΔt and can be obtained from the historically estimated results.

Equation (19) can be transformed into:

k=0γ-1m̂t-kSΔtρLk=0γm̂t-kSΔt=V* (20)

where V* is the virtual column of the pipe, and the relationship can be illustrated as Fig. 2.

Fig. 2 Schematic illustration of estimating time delay in a pipeline.

And γ can be calculated by (21) based on (19):

γ=minn ns.t.k=0nmt-kΔt/ρSLn0, n (21)

Then, the time delay s of each pipe can be estimated with the virtual volume V*

s=γΔt-V*-SLmt-γ/ρ (22)

For pipe i, the historically estimated inlet temperature which is closest to time t-s is selected as Tpi,iin(t-spi,i) in (23):

 Tpi,ioutt=(Tpi,iin(t-γpi,iΔt)-Ta)e-λpi,iCpρSpi,ispi,i+Ta  (23)

According to (21)-(23), the outlet temperature of a pipeline could be determined by the time delay of mass flow and historical inlet temperature:

Tpi,iout=Θi(spi,i,Tpi,iin(spi,i))i=1,2,...,Npi (24)

Usually, in the thermal network, the supply temperature at each source and the outlet temperature at each load are specified [

10]. Supposing that:

1) The mass flow rate m of each pipe is determined by the hydraulic model.

2) The historical inlet temperature Tpi,iin of each pipeline is known.

3) The delay of mass flow in each pipeline could be determined by (22).

Then, the thermal model could be performed independently by solving (10), (11) and (24).

C. Coupling Component

The electricity network and the DH network are coupled with coupling components at the boundary area, as shown in Fig. 3.

Fig. 3 Illustration of coupling area and coupling components of integrated network.

In this paper, only a few coupling components are modeled [

10]. The combined heat and power (CHP) units generate both electric power and heat power simultaneously. Gas turbines and internal combustion reciprocating engines can be modeled as in (25). When working in the condensing mode, extraction steam turbines can be modeled using (26).

ZCHP,igt=ΦCHP,igtPCHP,igti=1,2,...,NCHPgt (25)
ZCHP,ist=ΦCHP,istPCon,ist-PCHP,isti=1,2,...,NCHPst (26)

where ΦCHP,igt, PCHP,igt are the thermal power, the electric power output of a gas turbine of internal combustion reciprocating engine, respectively; ΦCHP,ist, PCHP,ist are the thermal power, the electric power output of an extraction steam turbine, respectively; PCon,ist is the electric power output of the extraction steam turbine in full condensing mode; NCHPgt, NCHPst are the numbers of the gas turbine of internal combustion reciprocating engine and the extraction steam turbine, respectively.

A circulation pump is propelled by a motor that consumes the electricity to create and maintain a pressure difference. The electricity consumption of the circulation pump is shown in (27).

Ppu,i=mpu,ig(Δppu,i)ηpu,ii=1,2,...,Npu (27)

where Ppu,i is the electric power consumption of the circulation pump; and ηpu,i is the efficiency of the pump.

Denoting (28) as the pressure energy supplied by a circulation pump [

30], (25)-(27) can be expressed in one equation as in (29).

Ψpu,i=mpu,ig(Δppu,i)i=1,2,...,Npu (28)
ΦCHP,igtΦCHP,jstΨpu,k=ZCHP,igt000-ZCHP,jst000ηpu,kPCHP,igtPCHP,jstPpu,k+0ZCHP,jstPCon,jst0  (29)

Denoting the power flow at the boundary area and the relative coefficient in matrix form, we can obtain (30)-(32), and represents the transition:

yE=PCHPgtPCHPstPpu (30)
yH=ΦCHPgtΦCHPstΨpu (31)
Γ1=ZCHPgt000-ZCHPst000ηpuΓ0=0ZCHPstPConst0 (32)

Then we can obtain:

yH=Γ1yE+Γ0 (33)

Equation (33) describes the coupling constraint of the boundary area of the combined network.

III. Approach of Hybrid SE (HSE)

A. Measurements and State Variables

SE of the combined heat and electricity network utilizes the measurements from both electricity network and DH networks. As for electricity network, the most commonly-used measurements include the magnitude of bus voltage denoted as U; active power and reactive power at both ends of a transmission line, denoted as Pb=[Pbf,Pbt] and Qb=[Qbf,Qbt]; and the electrical power injection at the bus bar, denoted as Pinj and Qinj. The state variables are the magnitude and phase angle of the bus voltage, denoted as U and θ, respectively. Note that one of the voltage phases in θ should be chosen as a reference phase angle. The measurement vector of the electricity network and the state vector are represented in (34) and (35), respectively.

zE=[UPbfPbtPinjQbfQbtQinj] (34)
xE=[Uθ] (35)

As for the DH network, the measurements are usually located at the heat source or consumer nodes as well as some important pipelines [

25]. In this paper, the measurement allocation problem of DH network is not discussed, thus we assume that the following four types of real-time measurements are available: ① nodal head pressure, denoted as p; ② nodal temperature including the temperature at the start and the end of each branch, denoted as T; ③ mass flow rate within the branch of the network, denoted as m; ④ heat power generated by a heat source or consumed by a heat consumer, denoted as Φ. This is a full measurement configuration.

In [

33], the nodal head pressure and injected flow are selected as the state variables in the water distribution system. However, they are not sufficient for thermal calculation in DHN as heat power transferred by mass cannot be determined by the nodal head pressure and mass flow rate in the pipelines. In this paper, m and T are selected as state variables of DHN. The measurement vector and state vector of DHN are represented in (36) and (37), respectively.

zH=[pTmΦ] (36)
xH=[mT] (37)

B. Formulation of SE

SE of IHEN is to estimate the solution of the following optimization problem with a given measurement set of the integrated network:

min J(xE,xH)=(rE)WErE+(rH)WHrH (38)

s.t.

rE=zE-fE(xE) (39)
rH=zH-fH(xH) (40)
Am=0 (41)
BKm2=0 (42)
A¯mTndout=A̲mTndin (43)
Tpiout=Θ(spi,Tpiin(spi)) (44)
yH(xH)=Γ1yE(xE)+Γ0 (45)

Equations(38),(45) describe the objective function of the SE problem using a weighted least square estimator [

20]; rE, rH are the vectors of measurement residual; and WE, WH are the weight matrices of measurement, which are the inverse of the variance matrices of the error of the meters. The variance matrices are given in advance, which are determined by the historical estimates and measurements. Equation (39) is the measurement function of the electricity network [20]; (40) is the measurement function of the DH network including (6), (11) and (12) [27]; (41) represents the flow continuity equation in (7), by setting mq=0; (42) represents the loop pressure equation, by substituting (6) into (8); (43) is the matrix form of (10), which represents the thermal constraints of water mixers, where A¯ and A̲ are the positive and negative incidence matrices which can be obtained from matrix A [34]; (44) is the matrix form of (24), which represents the thermal constraints of pipelines considering the transport delay; and (45) represents the constraints of the coupling components at the boundary area of the integrated network.

C. Two-stage Iteration Algorithm

To solve the optimization problem in (38)-(45), a moving horizon estimation approach with a two-stage iterative algorithm is adopted to solve the problem [

28] and the main steps of the algorithm are introduced as follows:

Step 1:   initialize all the state variables and parameters of the fluid.

Step 2:   with real-time measurement metered at time step t, add measurement functions in (39) and (40).

Step 3:   initialize the index of outer iteration, i = 1.

Step 4:   in the first stage, the time delay s of each pipeline is estimated using (22).

Step 5:   in the second stage, the entire SE problem in (38) is modeled and solved. Firstly, select the historically estimated inlet temperature of each pipeline which is closest to time ts, and add the dynamic thermal constraint of each pipeline in (44) into the model.

Step 6:   add other steady-state constraints into the model (37), including (41)-(43) of the DH network, and (45) of the coupling area.

Step 7:   solve the SE problem by Lagrange multipliers in (46) and get the estimates x̂t,i=x̂t,iE,x̂t,iH. An inner iteration algorithm is performed with the correct equation as in (47) [

35].

L(x,ω)=J(x)+ωTc(x) (46)
FTWFCTC0Δxω=FTWΔz-Δc (47)

where ω is the vector of Lagrange multipliers; W is the weight matrix of measurement; and F, C are the Jacob matrices of f(x) and c(x), respectively.

Step 8:   if the outer iteration is convergent, i.e., maxx̂t,i-x̂t,i-1<δ, or the index of the iteration i reaches the maximum number, go to the next time step, i.e., t=t+Δt, and repeat from Step 2; otherwise, update state variables with newly estimated results, set i=i+1, and repeat from the first stage in Step 4.

D. Integration of Different Measuring Resolutions

The discrepant measurement resolutions of the electricity network and the DH network should be considered when performing the combined SE. As for the electricity system, the monitoring resolution of SCADA system is usually seconds to minutes. Since the fast dynamics of the electricity system may cause nasty accidents, data are sampled with high frequency and transmitted with a special wireline network. However, in the DH network, the measurement data are still used for billing purposes in most circumstances, with the resolution of minutes to hours [

26]. In this paper, we assume that the measuring period of the electricity system, i.e., zE, is Δt, and that of the DH network, i.e., zH, is Δτ. Without the loss of generality, we assume that Δτ is integer multiple of Δt, i.e., Δτ=nΔt, where n is a positive integer.

During the measuring period of zH, the electricity network can perform an independent SE algorithm when zE arrives only, while the DH network does not perform SE until zH arrives. However, this method will lose the variation detail of the DH network in Δτ, and lose the accuracy of estimating the time delay s. That is because in (22), the discrete form of (16) assumes that the mass flow rate m does not change in Δt. When Δt is smaller, the result of the discrete form is closer to that of the integral form.

Although there is no available measurement, from the analysis in Section II-A, we found that with some specified conditions, the hydraulic network and thermal network can be solved independently. Therefore, we propose a method that uses only the measurements of the electricity network to estimate the states and dynamics of the DH network. The steps are introduced as follows:

Step 1:   perform the independent SE of the electricity system [

20] and obtain xE.

Step 2:   calculate the power flow of the boundary area of the electricity network as yE = yE(xE).

Step 3:   calculate the power flow of the boundary area of the DH network by using coupling constraints in (33), and obtain yH=ΦCHPgt,ΦCHPst,Ψpu.

Step 4:   estimate the mass flow rate of each circulation pump, i.e., mpu, by using (5), (28), and the power consumption of each circulation pump in yH, i.e., Ψpu.

Step 5:   solve the hydraulic network by using (6)-(8), and the mass flow rate of each circulation pump, i.e., mpu, and then obtain the mass flow rate of each pipe, source, and consumer.

Step 6:   estimate the time delay s of each pipe by using the mass flow rate of each pipe calculated in Step 5, i.e., mpi.

Step 7:   supposing that the outlet temperatures of all heat consumers, i.e., TL,spr , are specified, calculate the nodal temperature of the return network by using (10) and (23) as well as the time delay sr of each pipe of the return network estimated in Step 6.

Step 8:   estimate the supply temperature of each heat source TCHPr by using (11), the mass flow rate estimated in Step 5, the power injection of each source in yH, i.e., ΦCHP, and the return temperature of each source calculated in Step 7.

Step 9:   calculate the nodal temperature of the supply network by using (10) and (23), the supply temperature of each heat source TSs, and the time delay ss of each pipe of the supply network estimated in Step 6.

Step 10:   estimate the heat power of each consumer ΦL by using (12), the mass flow rate estimated in Step 5, the supply temperature estimated in Step 9 and the specified outlet temperature.

When ZE and ZH arrive at the same time, the integrated SE algorithm is performed for the entire network. The Algorithm 1 shows the HSE approach.

Algorithm 1  : HSE approach

Input: network topology, parameters

Output: estimates of state variables

1

Initialize state variables and time indices τ=0, t=0

2

while true do

3

τ=τ+Δτ

4

ifτ=t+Δtthen

5

  t=t+Δt

6

   Receive the measurements ztE, ztH

7

   Perform the two-stage iterative algorithm illustrated in Fig. 2       and get x̂t=x̂tE,x̂tH

8

else

9

   Receive the measurements zτE

10

   Perform the independent SE of the electricity system and get x̂τE

11

   Calculate yτE=yτE(x̂τE)

12

   Calculate yτH=Γ1yτE+Γ2

13

   Estimate mpu,τ of each circulation pump

14

   Solve the hydraulic network and get mτ

15

   Estimate the time delay s of each pipe

16

   Calculate the nodal temperature of the return network Tτr

17

   Estimate the supply temperature of each heat source TCHP,τs

18

   Calculate the nodal temperature of the supply network Tτs

19

   Estimate the heat power of each consumer ΦL,τ

20

End

21

End

IV. Numerical Tests

In this section, we validate the effectiveness of the proposed TSE approach and the HSE approach via numerical simulation on different cases. Firstly, the TSE approach is tested on a simple integrated network compared with the steady SE (SSE) approach proposed in [

27]. In this part, the measuring resolution differences are not considered. Then, the HSE approach is tested on a complicatedly integrated network with discrepant measurement resolutions of the electricity network and DH network.

In all of the test cases, some conditions in the network are specified, then true values of other states of the heating network are obtained by 12-hour simulation conducted by Bentley sisHYD®, and those of the electricity network are constructed by performing power flow calculation by Matpower under corresponding boundary conditions. At the boundary of the two networks, the coupling constraints are satisfied. Practical programs of different methods are developed in MATLAB, and the largest normalized residual test method for bad data identification is used in each program. The standard deviations for all the measurements are set equally as 0.01 p.u.. The convergence criterion is set as δ=0.0001.

Following indexes are used to measure the accuracy of the SE:

SM=1Ntt=1Nt1Nz(zt-f(x̂t))W(zt-f(x̂t)) (48)
SH=1Nt1Nz(f(x˜t)-f(x̂t))W(f(x˜t)-f(x̂t)) (49)

where Nt is the number of time points; Nz is the number of measurements; zt is the measurement; x˜t is the true value of the estimates; and x̂t is the estimates. Use SH/SM to evaluate the effects of SE, and a small value of the index means a good estimation.

A. Verification of TSE Approach

In this part, the proposed TSE approach is verified via numerical test on a simple integrated network. The integrated network consists of a 6-bus electricity network as shown in Fig. 4 and a simple heating network as shown in Fig. 5. In the heating network, there is only one heat source and one heat load connected by a supply pipe and a return pipe. The heat load is a heat exchanger. The CHP unit is a gas turbine unit, with ZCHP=1.3. The CHP unit is connected to bus bar 6, and the circulation pump is connected to bus bar 5 in the electricity network. Ts,in and Ts,out are the inlet and outlet temperatures of the supply pipe, respectively. Tr,in and Tr,out are the inlet and outlet temperatures of the return pipe, respectively. Parameters of the heating network in Fig. 5 are given in Table I [

28].

Fig. 4 One-line diagram of 6-bus electricity network.

Fig. 5 Structure of combined network.

TABLE I Parameters of Heating Network
PipelineL (m)S (m2)λ (W/(m·K))
Supply pipe 1000 4.1187×10-4 0.14
Return pipe 1000 4.1187×10-4 0.14

The SSE approach proposed in [

27] is used as a contrast. The measurements of the electricity network and the DH network are sampled with the same resolutions as 10 min, i.e., Δτ=Δt=10 min.

TABLE II SH and SM of TSE and SSE
ApproachSHSMSH/SM
SSE 1.30330 0.98625 1.32140
TSE 0.18712 0.98625 0.18972

Table II gives the estimate result of both approaches. The SH/SM of the estimates from the TSE approach is lower than the SH/SM from the SSE approach, which confirms that the TSE approach is more accurate than the SSE approach.

SH at each time point of both approaches are calculated (NT=1), and the results are depicted in Fig. 6. The true values of inlet and outlet temperature are shown in Fig. 7. Figures6 and 7 show that the error of the SSE approach is mainly from the adjusting period compared with the TSE approach, e.g., the 2nd hour and the 7th hour. The application of SSE will cause significant estimation error, especially during the dynamic process. The proposed TSE approach can accurately estimate system state, despite the large delay caused by water flow, which is more than 1 hour as shown in Fig. 7. Besides, the average computation time of the TSE approach is 0.5634 s, which proves the efficiency of the TSE approach.

Fig. 6 SH of estimated results at different time points.

Fig. 7 Measurement and true value of Ts,in and Ts,out.

B. Tests of HSE Approach

In this part, we illustrate the HSE approach with a 6-bus and 5-node IHEN. The 6-bus electricity network used here is the same one as shown in Fig. 4. The topology of the 5-node DH network used here is referred to [

36], which contains 5 pairs of pipes and 3 heat consumers. The electricity network and DH network are coupled with 2 CHP units. CHP1 is a gas turbine unit with Z1=1.3; CHP2 is an extraction steam turbine unit with Z2=8.1 and PCon=0.6. CHP1 and CHP2 are connected to bus bar 2 and bus bar 6 in the electricity network, and node 5 and node 4 in the DH network. The efficiencies of the circulation pumps are set equally as 0.85. The outlet temperatures of all the heat consumers are controlled locally and fixed at Tr=50  , and the ambient temperature is set as Ta=0  . The supply temperature of the heat sources and the heat power consumption of the heat loads are adjustable in the following cases. Measurements of the electricity network are sampled per 3 min, i.e., Δτ=3 min; while the DH network is measured per 30 min, i.e., Δt=30 min. Parameters of the pipelines are given in Table III. In this case, the heat power consumptions of the three consumers are time-invariant and equal to 0.3 MW. The heat power supply of CHP2 is fixed as 0.45 MW. The supply temperatures of the two heat sources are regulated with the following schemes:

TABLE III Parameters of Pipeline
PipelineL (m)S (m2)λ (W/(m·K))
Pipe 1 500 0.0177 0.14
Pipe 2 500 0.0177 0.14
Pipe 3 500 0.0177 0.14
Pipe 4 500 0.0177 0.14
Pipe 5 500 0.0177 0.14

1) Heat source 1: from 0-2 hours, 80 ℃; from 2-12 hours, 100 ℃.

2) Heat source 2: from 0-3 hours, 80 ℃; from 3-12 hours, 100 ℃.

The HSE approach introduced in Algorithm 1 is performed. The average computation time is 1.2216 s, and the SH/SM equals 0.085994, which shows the computation efficiency and the accuracy of the HSE approach. More detailed results are shown in Fig. 8.

Fig. 8 Test results of true value, measurements and estimates using HSE. (a) Power output and supply temperature of CHP1. (b) Power consumption and mass flow rate of pump 1. (c) Power output and supply temperature of CHP2. (d) Power consumption and mass flow rate of pump 2. (e) Mass flow rate of pipes 1, 3, 4, and 5. (f) Supply temperature of heat consumers 1, 2, and 3.

The supply temperatures of CHP units are regulated following curves (the third sub-figure of Fig. 8(a)) (CHP1) and (the third sub-figure of Fig. 8(c)) (CHP2). Since the thermal power output of CHP2 is fixed at 0.45 MW, CHP1 adjusts its thermal power output to satisfy the temperature regulation. Therefore, its electrical power output is changed. Besides, the circulation pumps are regulated co-ordinately to meet the demands of heat consumers at any time. The results are shown in Fig. 8. Note that measurements of the electricity network are sampled with high resolution (per 3 min); while in the DH network, the sampling period is 30 min. During a sampling period, the measurement cannot be updated. As shown in the results, the low sampling resolution will lead to the loss of perception of rapidly changing processes, especially when the regulations frequently occur. In the second sub-figure of Fig. 8(a), (b), and (e), the peak values or rapid variations are missed by the supervisor system of the DH network. The estimated results of using the hybrid estimation algorithm are represented with pink curves. Results show that this approach can utilize measurements with high resolution of the electricity network to estimate the dynamics of the DH network. The designed estimator and algorithm have desirable estimation performance.

The other advantage of HSE is to precisely estimate the time delay of heat transport in the pipeline by using the measurements in the electricity network. Here, the performance of estimating the time delay of each pipe is illustrated in Fig. 9. The method that only uses TSE illustrated in Section III-C is utilized as a contrast. The start-up period of TSE is equal to the measuring period of the DH network, i.e., Δt=30 min. Results show that HSE has better performance than TSE since it can make better use of the information provided by the measurements of the electricity network during the sampling period of the DH network.

Fig. 9 Estimates of time delay of each pipe using HSE, comparing with TSE. (a) Pipe 1. (b) Pipe 3. (c) Pipe 4. (d) Pipe 5.

V. Conclusion

In this paper, we firstly study the time-scale characteristics of IHEN. Then a TSE approach is proposed considering time delay caused by the process of heat power transportation in the pipeline. After that, a HSE approach is developed to integrate measurements of different networks with discrepant resolutions. Results show that, in both steady and dynamic processes, the two-stage estimator has good convergence. The hybrid estimator has good performance on tracking the variation of states in the heat network, even when the available measurements are very limited. When the heat network becomes active, the hybrid estimator shows obvious advantages.

The proposed estimation approach could be easily expanded when the dynamics of other components of the network should be considered. In our future work, dynamic properties of other components in the heat network shall be considered such as heat exchangers and buildings. Besides, the bad data identification with the dynamic model should also be considered.

References

1

P. Mancarella, “MES (multi-energy systems): an overview of concepts and evaluation models,” Energy, vol. 65, pp. 1-17, Feb.2014. [百度学术

2

H. Sun, Q. Guo, Z. Pan et al., “Energy internet: driving force, review and outlook,” Power System Technology, vol. 39, no. 11, pp. 3005-3013, Nov. 5, 2015. [百度学术

3

J. Wu, J. Yan , H. Jia et al., “Integrated energy systems,” Applied Energy, vol. 167, pp. 155-157, Apr.2016. [百度学术

4

Z. Li, W. Wu, M. Shahidehpour et al., “Combined heat and power dispatch considering pipeline energy storage of district heating network,”IEEE Transactions on Sustainable Energy, vol. 7, no. 1, pp. 12-22, Jan.2016. [百度学术

5

Z. Pan, Q. Guo, and H. Sun, “Feasible region method based integrated heat and electricity dispatch considering building thermal inertia,” Applied Energy, vol. 192, pp. 395-407, Apr.2017. [百度学术

6

J. Huang, Z. Li, and Q. H. Wu, “Coordinated dispatch of electric power and district heating networks: a decentralized solution using optimality condition decomposition,” Applied Energy, vol. 204, pp.1508-1522, Nov.2017. [百度学术

7

M. Ito, A. Takano, T. Shinji et al., “Electricity adjustment for capacity market auction by a district heating and cooling system,” Applied Energy, vol. 206, pp. 623-633, Nov.2017. [百度学术

8

L. Brange, J. Englund, and P. Lauenburg, “Prosumers in district heating networks–a Swedish case study,” Applied Energy, vol. 164, pp. 492-500, Feb.2016. [百度学术

9

H. Sun, Z. Pan, Q. Guo et al., “Energy management for multi-energy flow challenges and prospects,” Automation of Electric Power Systems, vol. 40, no. 15, pp. 1-8, Aug.2016. [百度学术

10

X. Liu, J. Wu, N. Jenkins et al., “Combined analysis of electricity and heat networks,” Applied Energy, vol. 162, pp. 1238-1250, Jan.2016. [百度学术

11

Z. Pan, Q. Guo, and H. Sun. “Interactions of district electricity and heating systems considering time-scale characteristics based on quasi-steady multi-energy flow,” Applied Energy, vol. 167, pp. 230-243, Apr.2016. [百度学术

12

I. Gabrielaitienė, R. Kačianauskas, and B. Sunden. “Thermo-hydraulic finite element modelling of district heating network by the uncoupled approach,” Journal of Civil Engineering and Management, vol. 9, no. 3, pp. 153-162, Jan.2003. [百度学术

13

A. Yan, J. Zhao, Q. An et al., “Hydraulic performance of a new district heating system with distributed variable speed pumps,” Applied Energy, vol. 112, pp. 876-885, Dec.2013. [百度学术

14

R. Kicsiny, J. Nagy, and C. Szalóki, “Extended ordinary differential equation models for solar heating systems with pipes,” Applied Energy, vol. 129, pp. 166-176, Sept. 2014. [百度学术

15

J. Duquette, A. Rowe, and P. Wild, “Thermal performance of a steady state physical pipe model for simulating district heating grids with variable flow,” Applied Energy, vol. 178, pp. 383-393, Sept. 2016. [百度学术

16

H. Lund, S. Werner, R. Wiltshire et al., “4th generation district heating (4GDH). integrating smart thermal grids into future sustainable energy systems,” Energy, vol. 68, pp. 1-11, Apr.2014. [百度学术

17

European Parliament. Directive 2012/27/EU of the European parlia ment and of the council. [Online]. Available: https:// eur-lex.europa.eu/LexUriServ/LexUriServ.do?uri=OJ:L:2012:315:0001:0056:EN:PDF [百度学术

18

H. Wei. “Central heating thermal stations in Beijing energy consumption status and energy saving effect analysis,” M.S. thesis, Beijing University of Civil Engineering and Architecture, Beijing, China, 2016. [百度学术

19

S. Yuan and W. Xu. “Establishment and validation of a sustainable evaluation model for heat metering technology in China,” Energy and Building, vol. 99, pp. 153-61, Jul.2015. [百度学术

20

A. Abur and A. G. Exposito, Power System State Estimation: Theory and Implementation. Boca Raton, CRC Press, 2004. [百度学术

21

F. C. Schweppe and J. Wildes, “Power system static state estimation,” IEEE Transactions on Power Apparatus Systems, vol. 89, pp. 120-135, 1970. [百度学术

22

A. Preis, A. J. Whittle, A. Ostfeld et al., “Efficient hydraulic state estimation technique using reduced models of urban water networks,” Journal of Water Resources Planning and Management, vol. 137, no. 4, pp. 343-351, Jul.2011. [百度学术

23

J. H. Andersen and R. S. Powell, “Implicit state-estimation technique for water network monitoring,” Urban Water, vol. 2, no. 2, pp.123-130, Jun.2000. [百度学术

24

C. T. C. Arsene and B. Gabrys. “Mixed simulation-state estimation of water distribution systems based on a least squares loop flows state estimator,” Applied Mathematical Modelling, vol. 38, no. 2, pp. 599-619, Jan.2014. [百度学术

25

T. Fang and R. Lahdelma, “State estimation of district heating network based on customer measurements,” Applied Thermal Engineering, vol. 73, no. 1, pp. 1211-1221, Dec.2014. [百度学术

26

H. Wang, H. Meng, and T. Zhu, “New model for onsite heat loss state estimation of general district heating network with hourly measurements,” Energy Conversion and Management, vol. 157, pp. 71-85, Feb.2018. [百度学术

27

J. Dong, Q. Guo, H. Sun et al., “Research on state estimation for combined heat and power networks,” in Proceedings of2016 IEEE PES General Meeting (PESGM), Boston, USA, Jul.2016, pp. 1-5. [百度学术

28

T. Sheng, Q. Guo, H. Sun et al., “Two-stage state estimation approach for combined heat and electric networks considering the dynamic property of pipelines,” Energy Procedia, vol. 142, pp. 3014-3019, Dec. 2017 [百度学术

29

T. Zhang, Z. Li, Q. H. Wu et al., “Decentralized state estimation of combined heat and power systems using the asynchronous alternating direction method of multipliers,” Applied Energy, vol. 248, pp. 600-613, Aug.2019. [百度学术

30

CRANE, Flow of Fluids Through Valves, Fittings and Pipe. New York: Metric Edition, 1942. [百度学术

31

H. Zhao, “Analysis, modelling and operational optimization of district heating systems,” Ph.D. dissertation, Technical University of Denmark, Denmark, 1995. [百度学术

32

L. Dobos, J. Jäschke, J. Abonyi et al., “Dynamic model and control of heat exchanger networks for district heating,” Hungarian Journal of Industrial Chemistry Veszprem, vol. 37, pp. 37-49, Jan.2009. [百度学术

33

M. J. H. Sterling and A. Bargiela, “Minimum norm state estimation for computer control of water distribution systems,” IEE Proceedings of Control Theory and Applications, vol. 131, no. 2, pp. 57-63, Mar.1984. [百度学术

34

J. Wang, Z. Zhou, and J. Zhao, “A method for the steady-state thermal simulation of district heating systems and model parameters calibration,” Energy Conversion and Management, vol. 120, pp. 294-305, Jul.2016. [百度学术

35

W. Lin and J. Teng, “State estimation for distribution systems with zero-injection constraints,” IEEE Transactions on Power Systems, vol. 11, no. 1, pp. 518-524, Feb.1996. [百度学术

36

X. Qin, H. Sun, X. Shen et al., “A generalized quasi-dynamic model for electric-heat coupling integrated energy system with distributed energy resources,” Applied Energy, vol. 251, no.1, Oct.2019. [百度学术