Model predictive control based robust scheduling of community integrated energy system with operational flexibility

Applied Energy 243 (2019) 250–265

Contents lists available at ScienceDirect

Applied Energy journal homepage: www.elsevier.com/locate/apenergy

Model predictive control based robust scheduling of community integrated energy system with operational flexibility

T

Chaoxian Lva, Hao Yua, Peng Lia, , Chengshan Wanga, Xiandong Xub, Shuquan Lic, Jianzhong Wub ⁎

a

Key Laboratory of Smart Grid of Ministry of Education, Tianjin University, Tianjin 300072, China Institute of Energy, School of Engineering, Cardiff University, Cardiff CF24 3AA, UK c State Grid Customer Service Centre, Tianjin 300300, China b

HIGHLIGHTS

flexibilities of multiple links are utilized with economic performance improved. • The adaption to uncertainties is enhanced by the MPC-based rolling optimization. • The • The scheduling strategy with reservation margins is robust for real-time errors. ARTICLE INFO

ABSTRACT

Keywords: Community integrated energy system (CIES) Flexibility Source-network-load Energy storage Model predictive control (MPC)

Community integrated energy system (CIES) couples multiple energy forms and links, it is important to cope with the inadequacy of flexibility for realizing efficient utilization and reliable supply. Based on the detailed modeling of central energy station, district heating network (DHN) and building loads, a unified scheduling model is established by fully considering the coupling characteristics of the sources, the networks, and the loads. Meanwhile, the flexibilities from adjustable capacity of thermal storage device and the indoor temperature on the demand side, as well as the energy storage of DHN, are utilized to improve the economy. Furthermore, a model predictive control (MPC) based robust scheduling strategy is proposed to maintain and utilize CIES flexibility for enhancing the uncertainty adaptability. The scheduling framework consists of rolling optimization and robust constraint generation. The rolling optimization schedules the system with a better adaptation to flexibility demands, and the robust constraints generate adjustable margin for building demands by modifying the upper/lower and ramping limits. After the linearization of nonlinear terms in the model, case studies are conducted based on the data of a typical day in summer. The results show that the unified source-network-load scheduling strategy can give full play to the flexibilities of different energy links with a better operation economy. Additionally, the MPC-based optimization can well adapt to the rolling forecasting uncertainties, and the schedule generated by advance reserving means has stronger robustness for real-time errors.

1. Introduction With the primary energy depletion and increasing environmental pressures all around the world [1], improving overall energy efficiency and reducing environmental pollution have become urgent tasks for realizing the sustainable development [2]. Community integrated energy system (CIES) deeply integrates various energy forms and supplies flexible energy demands such as cooling, heating and electricity to end users [3]; making it a promising means to build an efficient and green energy system via unified planning and operation strategies [4,5]. CIES contains multiple energy forms and the supply, delivery and

⁎

consumption are highly coupled, making the operation behavior more complex [6]. It’s difficult to maintain the supply/demand balance at all times; moreover, the variability and uncertainties from renewable energy systems (RESs), loads and some other links pose a growing challenge to reliable operation [7,8]. The coordination to meet the real-time energy balance and the ability to cope with uncertainties are usually regarded as the flexibility of CIES [9]; in addition, flexibility is helpful to promote the economy and utilization efficiency [10]. Different energy links in CIES have their own flexibilities, such as the combined heat and power (CHP) units/storage devices in the source, the heating/ gas network in the delivery link and flexible loads on the demand side

Corresponding author. E-mail address: [email protected] (P. Li).

https://doi.org/10.1016/j.apenergy.2019.03.205 Received 26 December 2018; Received in revised form 26 March 2019; Accepted 31 March 2019 0306-2619/ © 2019 Elsevier Ltd. All rights reserved.

Applied Energy 243 (2019) 250–265

C. Lv, et al.

QtDC,C / QtDC,I cooling/ice-making power of the i th DC at time t (kW) ,i ,i DC,C Ut , i / UtDC,I cooling/ice-making mode of i th DC at time t ,i mass flow rate of the i th pipe at time t (kg/s) mi , t Tis/r,start / Tis/r,end temperature at the start and end of the i th supply/ ,t ,t return pipe at time t (°C) indoor temperature of building i at time t (°C) Ti, t out Qiin injecting and dissipation power of building i at time t (kW) , t /Qi, t Qtces cooling power of CES at time t (kW) Ttces,s/Ttces,r supply/return temperature of CES at time t (°C) B,out TiB,in inlet and outlet water temperature of building i at time t , t /Ti, t (°C)

Nomenclature Abbreviations HP WC IS/IT DC EP/CWP CES CIES DHN MPC

ground source heat pump conventional water-cooled chiller ice-storage system/ice-storage tank double-duty chiller ethylene glycol pump/chilled water pump central energy station community integrated energy system district heating network model predictive control

Parameters

Sets

Q _ i

HP / WC / DC set p,s/r p,s/r i,in / i,out set of

p CES /

p B, i

Q

of HP/WC/DC units the supply/return pipelines flowing in/out of

HP

/ Q iHP lower/upper power limit of the HP (kW)

WC

/ QiWC lower/upper power limit of the WC (kW)

_ i

W IT/ Q IT upper limit of the stored energy/cooling power of the IT (kWh/kW) Q DC,C / Q iDC,C lower/upper cooling power limit of the DC (kW)

nodei set of pipelines connected to CES/buildingi

_ i

Variables

Q

DC,I

_ i EP

/ QiDC,I lower/upper ice-making power limit of the DC (kW)

Q / Q IS,CWP upper power limit of the EP/CWP collateral to the IS (kW) lower/upper limits of indoor temperature (°C) T /T

WC QtHP cooling power of the i th HP/WC unit at time t (kW) , i / Qt , i HP Qt /QtWC cooling power of the HPs/WCs at time t (kW) WC DC UtHP ON/OFF state of the i th HP/WC/DC at time t , i / Ut , i /Ut , i PtHP/ PtWC/ PtIS consumed power of the HP/WC/IS subsystem at time t (kW) WtIT cooling energy stored in IT at time t (kWh) QtIT /QtIS cooling power of the IT/IS at time t (kW) UtIS,CWP cooling state of the i th CWP collateral to the IS at time t ,i QtDC,C/ QtDC,I cooling/ice-making power of the DCs at time t (kW)

_

T ave,min /T ave,max lower/upper limits of average indoor temperature for each building (°C) T ces,min / T ces,max lower/upper limit of supply temperature of CES (°C) T ces,r,max upper limit of CES return temperature (°C)

[11]. The highly coupling characteristics of source-network-load in CIES bring new challenges to traditional operation strategies, which consider the flexibility of different links independently [12]. To improve the overall operation flexibility, the generation, distribution and consumption process of the system need to be reasonably coordinated by a unified scheduling framework based on flexibility modelling of different links. At present, many studies on the supply side flexibility of CIES have been conducted. In [13], a CHP unit, ice-storage system and other units worked together to supply space cooling/space heating/hot water for a building system. This system had a higher energy utilization rate and economy due to the flexibility of the energy storage device and CHP unit. The day-ahead scheduling of a cooling/electricity energy system including a cold water tank and CHP unit was studied in [14], and the introduction of the thermal storage device decreased the operation cost under the mechanism of dynamic electricity price. In [15], a bi-level economic dispatch was carried out for an integrated natural gas and electricity system considering wind power and power-to-gas (P2G) process, and the study proved that P2G process is effective on improving the wind power utilization rate. Ref. [16] presented an optimization model to achieve the potential increase in performance and flexibility by coordinating the operation of heat pumps and thermal energy storage of a heating system. And the system became cost competitive under TOU tariff. Refs. [17] and [18] conducted the flexibility evaluation of the integration of thermal storage devices from multiple dimensions, including time, energy, power, and cost. And the operational flexibilities from energy supply/storage devices are quantified for the improvement of economic performance. Besides the flexibility utilized from device level, energy networks can also introduce more options for providing additional flexibility in system operation, as well as more realistic scheduling strategies. Ref. [19] introduced a day-ahead

scheduling of integrated urban energy system and coordinated with the network reconfiguration of distribution network (DN). In [20], an analytical framework for quantifying the flexibility of active distribution networks was proposed, including the quantification of node flexibility and network flexibility. Also, soft open points and energy storage systems were used to improve the indexes of flexibility including magnitude, frequency and intensity dimensions. In [21], a dispatch strategy that coordinated the heat and electrical power of a multi-energy system was proposed, which can improve the flexibility and economic performance by the couple of district heating network (DHN) and DN. The authors of [22] had found that the amount of flexibility of the district heating system is determined by the delayed or forced operation mode. The authors in [23] proposed the dynamic temperature distribution model of DHN for integrated heat and power dispatch. The CHP flexibility was increased by the delay characteristic of DHN. In [24], a heat/electricity system was optimized to minimize the operation cost and the wind curtailment penalty with considering the energy storage capacity of DHN, and better economy and wind utilization rate were achieved. In [25], the dynamic characteristics of gas networks were modelled, and it was proved the line pack storage can provides major option space for gas system. For an integrated electricity and gas system, the authors of [26] considered different heating scenarios and presented a method to quantify the flexibility that the gas network can provide. In addition, some literature considered demand-side management (DSM), which can change the load profiles by multiple means [27,28], to achieve a flexible operation. For an integrated energy system with energy storage, Ref. [29] conducted the optimal operation by utilizing the complementarity of PV energy time-shift and demand load-shift, and achieved a better economy. In [30], a charging pricing strategy of electric vehicle (EV) stations was proposed, which consists of the simulation model of EV mobility and a 251

Applied Energy 243 (2019) 250–265

C. Lv, et al.

double-layer optimization mode. And the strategy fully exploited EV’s flexibility by price scheme to improve the voltage profiles of DN. Ref. [31] considered the thermal inertia of buildings and the thermal comfort of consumers, and the DSM of a heating/electricity system was conducted. Thus, the system could accommodate more wind power. Ref. [32] carried out optimal scheduling with the virtual energy storage by regarding buildings as flexible and controllable units, and the operation economy was promoted. When facing a system with flexible energy devices, heating/gas network and adjustable demand side resources, it is of necessity to utilize the flexibility of different links for a better operation performance. Besides better economy and higher RES utilization rates, flexibility offers an efficient way to improve the adaption to uncertainties. During day-ahead scheduling, the optimal strategy assumes an accurate forecast for RESs, environmental factors and loads. However, the forecast information may deviate from the actual values in practical operation, causing a lack of flexibility with increasing time scale. Model predictive control (MPC) based on rolling optimization and feedback correction can adapt well to the problem of uncertain forecasting, and it has advantages compared with multi-time advance scheduling such as stochastic planning [33] and interval planning [34]. In [35], an MPC method was adopted for coordinating the operation of energy supply and storage devices in a cooling/electricity system. In [36], after constructing typical scenarios of electricity price, a stochastic MPC method was adopted to conduct energy management of a microgrid to reduce the negative impact of uncertain information. For a grid-connected microgrid, Ref. [37] modeled the variabilities of RES and loads, and implemented a rolling method based on MPC approach to capture the interaction behaviour with DN. Thus, RES and load forecast errors were well balanced. Although the above literatures have enhanced adaptability to inaccurate forecasting information by using MPC method, the errors of forecasting and actual operation still exist and are not taken into account. The authors in [38] carried out rolling optimization with building storage, but the operation should be adjusted to maintain indoor temperatures consistent with the scheme when facing real-time fluctuations. For a combined cooling, heating, and power (CCHP) microgrid, the multiple time-scale scheduling method was used to reduce the impact of real-time errors in [39]. The methods of [38,39] belong to post-checkout scheduling with no uncertainty quantification for flexibility supplements, and they have high requirements for equipment responsiveness, data transmission timeliness and accuracy and have difficulties in engineering applications. Therefore, the demand side flexibility should be further utilized to adapt to the real-time uncertainties by advance scheduling method. In this paper, we propose an MPC-based robust scheduling method considering multiple flexibilities from source-network-load to adapt the uncertainties from the forecasting data. First, the operation model of central energy station (CES), dynamic transfer model of DHN and the storage model of buildings are built; then a unified source-network-load scheduling is carried out to increase flexibility for adaption to the uncertainties by using the robust framework, which consists of MPC-based rolling optimization and robust constraint generation of the control domain. Finally, the correctness and effectiveness of the system model and scheduling strategy are verified through case analysis. The main contributions of this paper are summarized as follows:

complementary operation of the supply and storage devices in the CES, the adjustable indoor temperature of building loads, as well as the dynamic transfer characteristics of DHN. Thus, the operation cost is decreased by the flexible operation of different links under TOU mechanism. (3) The proposed MPC-based robust strategy can adapt well to the uncertainties of the environment, RESs and loads. With rolling optimization for better adaption to the flexibility required in multitime scheduling and robust constraint generation for the uncertainty in the control domain, robust scheduling can be realized, achieving better demand-side management with the economy guaranteed. The rest of this paper is organized as follows. Section 2 describes the system structure and configuration of the CIES. Section 3 establishes the operation model of the central energy station, heating network and building loads. Section 4 proposes the MPC-based robust scheduling model integrated with source-network-load. Section 5 selects a typical day to conduct the case analysis. Finally, Section 6 presents the conclusions of this paper. 2. System structure and configuration This paper selects the CIES at State Grid North Customer Service Centre in Tianjin as the study case. The project has a total construction area of 143,000 m2, including one monitoring building, one public service building, two call buildings and five staff dormitory buildings. In the system, electricity is the core form of energy conversion and the system driver of the CIES. The system has multiple energy demands including space heating/space cooling and electricity, and different types of renewable energy sources, such as photovoltaic (PV) and ground source heat pumps (HPs). Additionally, the CIES integrates multiple energy processes, such as energy supply (central energy station), energy delivery (district heating network) and energy consumption (space heating/space cooling for buildings), as well as the energy storage process of thermal storage devices in the CES. We mainly focus on the scheduling problem for cooling and electricity supplement in the cooling period. In this period, the external grid and PV system meet the electricity demand together. The CES produces chilled water for air conditioning. The chilled water is delivered through the supply pipelines to each building and the fan coil systems of the buildings work to supply space cooling. The CES can be divided into three energy subsystems: the ground source heat pump system (three ground source heat pump units), conventional water-cooled chiller system (two conventional water-cooled chillers (WCs)), and icestorage system (IS) which contains two double-duty chillers (DCs) and one ice-storage tank (IT). There are alternative patterns of energy supply in the system. The cooling demand can be directly supplied by electrically driven chillers, as well as by energy storage devices [40]. The configuration of the CES and corresponding energy flow are shown in Fig. 1, and the structure of DHN is shown in Fig. 2. 3. Modelling of the CIES operation 3.1. Modelling of central energy station

(1) A unified source-network-load scheduling strategy is proposed to improve the overall performance of the CIES operation. Based on the detailed modeling of CES, DHN and building loads, the coupling characteristics of different energy links are integrated into the strategy. And the coordinated operation of sources, networks and loads can be achieved by unified dispatching, which is essential for actual operation. (2) The flexibility potentials of source-network-load are fully utilized with the economy increased. The strategy can give full play to the

3.1.1. Ground source heat pump The constraints on the ground source heat pumps’ cooling power are shown in Eqs. (1)–(3). The cooling powers of the ground source heat pumps in operation are kept equivalent, as shown in Eq. (1). The relationship between the cooling powers of ground source heat pumps and each unit is shown in Eq. (2). Eq. (3) describes the lower/upper limits of the cooling power of the units. N HP is the number of ground source heat pump units. 252

Applied Energy 243 (2019) 250–265

C. Lv, et al.

Pt PV

PV

HP subsystem

LEt Ground source heat pump 1 Ground source heat pump 2

Pt HP

QtHP

Ground source heat pump 3

External grid

Pt TL

WC subsystem

Conventional water-cooled chiller 1

Pt WC

Conventional water-cooled chiller 2

Double-duty chiller 1

QtDC,C DC,I

Double-duty chiller 2 / Qt

Wt

IT

Ice-storage

QtWC

PHE

Ice-storage subsystem

Pt IS

Electricity

ces t

Q

IS t

Q

IT t

Q

tank

Cooling

QtDC,I water distributor

Central energy station

Elctricity

Water pump

Chilled water

Ethylene glycol pump

Cooling water

Plate heat exchanger

Ethylene glycol

PHE

Fig. 1. Configuration and energy flow of the central energy station [40]. N HP

QtHP ,i

HP HP UtHP , i = Ut , i Qt , i=1

i

UtHP ,i Q HP

(1)

QtHP ,i i=1

QtHP ,i

HP UtHP , i Qi ,

i

HP

(3)

The operation mode/on-off state of the chillers and pumps is 1 when it is in the execution/on state and 0 when in the non-execution/off state. Eq. (4) restricts the startup/shutdown orders of the ground source heat pumps.

N HP

QtHP =

HP

_ i

(2)

Fig. 2. Structure of DHN. 253

Applied Energy 243 (2019) 250–265

C. Lv, et al.

UtHP ,i

UtHP , i + 1,

i = 1, 2,

, N HP

(4)

1

UtDC,C Q ,i

Eq. (5) describes the power consumption of the ground source heat pumps, where COP HP is the coefficient of performance (COP) of the i th i ground source heat pump; P HP,CWP and P HP,CP are the rated powers of the chilled water pump (CWP) and the cooling water pump (CP) interlocked with the ground source heat pump, respectively. N

PtHP =

HP

UtDC,I Q ,i

HP,CP + PHP,CWP ) UtHP , i (P

(5)

=

WC UtWC , , i Qt

i

WC

(6)

i=1

N WC

QtWC =

QtWC ,i

(7)

i=1

UtWC ,i Q

WC

QtWC ,i

UtWC ,i

UtWC , i + 1,

_ i

WC UtWC , , i Qi

i = 1, 2,

i

WC

, N WC

(8)

The power consumption of the conventional water-cooled chillers is depicted in Eq. (10), where COPiWC is the COP of the i th conventional water-cooled chiller and PWC,CWP, PWC,CP and PWC,CT are the rated powers of the CWP, CP and cooling tower (CT) interlocked with conventional water-cooled chiller, respectively.

PtWC =

N WC i=1

WC QtWC + UtWC , i /COPi ,i

(PWC,CWP + PWC,CP + PWC,CT)

i

DC

DC

(16) (17)

(18)

1

UtDC,C ,i

UtDC,C ,i +1 ,

UtDC,I ,i

UtDC,I , i + 1,

, N DC

i = 1, 2,

, N DC

i = 1, 2,

UtDC ,i

UtDC,C , ,i

i

UtDC ,i

UtDC,I ,i ,

i

UtDC ,i

UtDC,C + UtDC,I ,i ,i ,

1

(19) (20)

1

(21)

DC

(22)

DC

i

(23)

DC

The cooling and ice-making modes of the double-duty chillers cannot operate simultaneously, as depicted in Eq. (18). Eqs. (19) and (20) constrain the startup/shutdown orders of the double-duty chillers. Eqs. (21)–(23) describe the relationship between the ON/OFF state of each double-duty chiller and its cooling/ice-making mode. A doubleduty chiller is on when operating in the cooling or ice-making mode and off otherwise. Eqs. (24) and (25) describe the relationship between the stored cooling energy and charge/discharge power of the ice-storage tank, and restrict the total cooling energy stored in the ice-storage tank within its limits, where IT is the heat loss rate of the ice-storage tank. From Eq. (26), we can see that the cooling power of the ice-storage tank should be lower than its upper limit and that the ice-storage tank cannot release cooling when the double-duty chillers are operating in ice-making mode.

(9)

1

DC,I UtDC,I , , i Qi

QtDC,I ,i

UtDC,C + UtDC,I ,1 ,1

N WC

UtWC ,i

DC,I

i

Eqs. (12) and (13) ensure that the power outputs of the double-duty chillers in cooling mode and ice-making mode are the same. The relationships between the cooling/ice-making power of double-duty chillers and each unit are shown in Eqs. (14) and (15). Eqs. (16) and (17) restrict QtDC,C /QtDC,I within their lower and upper bounds when ,i ,i operating in the cooling/ice-making mode.

3.1.2. Conventional water-cooled chiller Eq. (6) ensures that the cooling powers of conventional watercooled chillers in operation remain consistent. The relationship between the cooling powers of conventional water-cooled chillers and each unit is shown in Eq. (7). Eqs. (8) and (9) describe the upper/lower power limits and startup/ shutdown order of the units, respectively. N WC is the number of conventional water-cooled chillers.

QtWC ,i

UtDC,C QiDC,C, ,i

QtDC,C ,i

_ i

HP QtHP , i /COPi +

i=1

DC,C

_ i

(10)

3.1.3. Ice-storage system The double-duty chiller includes two main operating modes: cooling mode and ice-making mode. The ice-storage tank makes use of the water/ice phase change to achieve the charge and discharge of cooling energy.

N DC

WtIT = (1

IT ) W IT t 1

QtDC,I t ,i

+ i=1

WtIT

W IT _

W IT

QtIT t

(24) (25)

NDC

QtDC,C + QtIT = QtIS ,i

0

(11)

i=1

N DC

UtDC,C = UtDC,C QtDC,C, ,j ,i

i

DC

j=1

N IS,CWP

DC,I UtDC,I = UtDC,I , ,j , i Qt j=1

i

DC

(12)

i=1 N IS,CWP

(13)

i=1

N DC

QtDC,C =

QtDC,C ,i i=1

QtDC,I ,i i=1

(26)

N DC

UtDC,C ,i i=1

UtIS,CWP Q IS,CWP ,i

QtIS

(27)

(28)

Eq. (29) describes the power consumption of the ice-storage system, where COP DC,C and COP DC,I are the COP of the i th double-duty chiller in i i the cooling mode and ice-making mode, respectively, and P DC,CP , P DC,CT , P EP and P IS,CWP are the rated powers of the CP, opening CT, ethylene glycol pump (EP) and CWP in the ice-storage system, respectively.

(14)

NDC

QtDC,I =

IT UtDC,I ,1 ) Q

UtIS,CWP ,i

N DC

QtDC,I ,i

(1

Eq. (27) restricts the running number of CWPs in the ice-storage system according to the principle of “one chiller-one pump”, and Eq. (28) guarantees that the CWPs in operation can satisfy the cooling supply requirement of the ice-storage system.

The sum of the cooling power of the double-duty chillers and the ice-storage tank is equal to that of the ice-storage system, which is depicted in Eq. (11), where N DC is the number of double-duty chillers.

QtDC,C ,i

QtIT

(15) 254

Applied Energy 243 (2019) 250–265

C. Lv, et al.

PtIS

=

in Eq. (33). Eq. (34) denotes that the starting temperature of the pipeline flowing out of the intersection node is equal to the mixing temperature.

DC,I /COPiDC,C + QtDC,I QtDC,C + ,i , i /COP i

N DC

(PDC,CP + P DC,CT)+ UtDC,C ,i DC,I Ut , i (PDC,CP + PEP + P DC,CT)

i=1

+

EP QtIS P / Q EP

N IS,CWP

UtIS,CWP P IS,CWP ,i

+ i=1

k

(29)

i

+ Ta

s. t.

Ri, t =

k

p,x i,in

k

p,x i,out

(33) (34)

ni, t k=0

Z

mi , t + k · t

Ai Li

(35)

mt + k · t

ni, t 1 k=0

0,

(36)

mt + k · t ,

if ni, t

if ni, t < 1

1 (37)

Since the temperatures of both sides of the pipeline are adjusted once per schedule interval, the number of delay intervals is usually not an integer. Different from the rounding method with large errors [42], the weighted delay method is adopted in this paper. The weighted delay coefficient qi, t and parameter i1, t , i2, t are defined as follows:

(2) Flow balance: The incoming and outgoing mass flows of each node are the same, as shown in Eq. (32). In the actual energy supply system, the mass flows of the supply and return pipelines with a consistent index are the same.

mk, t

ni, t k=0

Si, t =

where x {s, r} ; s and r represent the supply and return pipelines, respectively; i is the heat transfer coefficient of pipeline i ; Li is the length of pipeline i ; Cw is the specific heat capacity of water; Ta is the ambient temperature of the pipelines; di is the internal diameter of pipeline i ; and µ is the heat loss coefficient of the pipelines.

mk, t =

mk, t

p,x i,out

k

min ni, t , ni, t

(30) (31)

= di µ

p,x i,out

In Fig. 3, mi, t t is the water quality injecting to the pipeline i at time t; , Ai , and Li are the water density, cross-sectional area and length of pipeline i , respectively. We assume that the water mass flowing into pipe i at time t flows out for the first time at time t + ni, t t , where ni, t can be obtained by Eq. (35). Parameter Si, t denotes the water quality flowing into pipeline i between t and t + ni, t t ; similarly, Ri, t denotes the water quality flowing into pipeline i between t and t + (ni, t 1) t , as shown in Eqs. (36) and (37).

(1) Energy loss: In the process of transmission, there are differences between the medium (cold/hot water) and ambient temperature of the pipeline, and there is energy exchange through the tube wall of the pipeline. The energy loss is represented by the decrease (heating) and increase (cooling) of the heat medium temperature from the start to the end of the pipeline, as shown below: i Li Cw mi, t

k

(4) Transfer delay: The water flows in the pipeline at a certain speed, so the temperature change at the end of the pipeline has a delayed effect relative to that at the start. The transmission process can be viewed in Fig. 3.

The DHN transfers the energy generated by heating/cooling sources to the demand side through pipelines and heat media (hot water/steam, cold water) [41], and then meets the user's energy demand through the work of the radiator or fan coil unit. This paper mainly focuses on the modelling of hot/cold water networks.

Ta ) e

mk, t Tkx,end = Tix,mix ,t ,t

Tkx,start = Tix,mix , ,t ,t

3.2. Dynamic transfer modelling of DHN

Tix,end = (Tix,start ,t ,t

p,x i,in

qi, t =

(32)

(3) Temperature mixing: According to the first law of thermodynamics, the energy flowing in and out of each node is conserved. The mixing temperature at the intersection node of different pipelines is shown

1 i, t

=

1 i, t 1 i, t

+

2 i, t

Ai Li Ri, t , mi, t + n 1

if ni, t

Ai Li

if ni, t < 1

Ri , t

mi, t + n

Fig. 3. Pipeline transfer diagram. 255

(38)

,

1 (39)

Applied Energy 243 (2019) 250–265

C. Lv, et al.

2 i, t

=

Si, t

Ai Li mi , t + n

Qtces = Cw

(40)

where i1, t and i2, t indicate the time errors between the delay interval using ni, t , ni, t + 1 and the standard delay respectively. Therefore, the temperature of the ending node can be represented by the weighted form of two adjacent delay intervals, and the energy loss equation can be modified as follows:

qi, t · T

x,end i, t + ni, t + 1

x,start qi, t )· Tix,end , t + ni, t = (Ti, t

+ (1

Ta ) e

i Li Cw mi, t

+ Ta

Ttces,r = Tir,mix ,i ,t

(41)

1

=

Qiout , t = (Tout, t

(Tout, t

Ti, t 1) Ki Fi

T

T

T

T

T ave,min

Ti, t

Ti, t

1

NT

Ti, t / NT

T

T ave,max

t=1

Ti,ini = Ti,ter

T ces,max

T ces,min

Ttces,r

T ces,r,max

TiB,out ,t

(52) (53)

=

T jr,start , ,t

TiB,in ,t , j

p B,i

p B, i

j

p B, i

(54) (55) (56)

3.5. Power balance The cooling balance of CES is shown in Eq. (57), and the electricity balance of the system is described in Eq. (58).

(43)

NHP

(44)

N WC

QtHP ,i + i=1

IS ces QtWC , i + Qt = Qt i=1

(57)

PtTL + PtPV = LtE + PtHP + PtWC + PtIS

(58)

PtTL

(59)

PtTL,max

where is the electricity load except for the electricity needed in the CES at time t (including the power consumed by the pumps of DHN), and PtPV is the power output of the PV system at time t. The purchasing power PtTL should not exceed the maximum allowed power of the tieline P TL,max (see Eq. (59)).

LtE

(45)

i, t

Ttces,s

TiB,in = T js,end ,t ,t , j

where Tout, t is the outdoor temperature at time t, and Ki and Fi are the average heat dissipation coefficient and surface area of building i , respectively. To ensure better user experiences, the indoor temperature of the building should remain within the comfort temperature range, and the indoor temperature change rate at adjacent times and the average temperature of a complete scheduling period should be limited, as shown in Eqs. (45)–(47). Eq. (48) ensures that the terminal temperature is the same as the initial temperature in a scheduling cycle. _

T ces,min

B,out Qiin , t = Cw mj, t Ti, t

Qiin ,t

air · Vi

(51)

The power injected to the can be calculated by Eq. (54). Eqs. (55) and (56) describe that the ending temperature of the supply pipelines and the starting temperature of the return pipelines connected with the building are the same with the inlet and outlet temperatures, respectively.

(42)

Ti, t 1)· Ki ·Fi Cair ·

N col

building Qiin ,t

where Cair and air are the specific heat capacity and density of air, respectively. To meet the scheduling requirement, the differential equation needs to be discretized. The dynamic temperature equation for buildings after discretization can be expressed as:

Ti, t t

(50)

is the set of collector-distributor nodes, and the supply/rewhere turn temperature Ttces,s /Ttces,r can also be regarded as the media temperature at the water distributor/water collector. To ensure the quality of the energy supply, the supply and return water temperatures of the CES should be less than the upper temperature limit; considering the refrigeration characteristics of the watercooled equipment, the water supply temperature should be greater than the lower limit, as shown in Eqs. (52) and (53).

The indoor temperature of buildings can change within a certain range. At this point, buildings can be considered to be demand response resources with energy storage ability. Taking the cooling season for example, the heat balance equation of a building is obtained according to the law of energy conservation [43], as shown in Eq. (42):

Ti, t

(49)

N col

3.3. Demand response modelling of buildings

Cair

p CES

Ttces, s, i

p CES

Ttces,s = Tis,start ,i ,t

The above weighted delay method, which can avoid the problem of large delay mismatch caused by the accumulation of traditional rounding errors, is used to simulate the dynamic temperature distribution of the nodes in the network, and this method is more suitable for practical applications.

dTi, t = Qiout Qiin ,t ,t air Vi dt

mi, t Ttces, r i

4. MPC-based robust scheduling with source-network-load

(46)

In order to adapt to the predictive information inaccuracies of environment factors, RESs and loads in multiple periods and the scheduling robustness requirement for real-time uncertain information, this paper adopts flexibility-oriented MPC optimization method for system scheduling, including MPC-based rolling optimization and robust scheduling for the execution/control domain.

(47) (48)

where NT is the number of scheduling intervals for a complete scheduling cycle; T is the ramping limit; Ti,ini and Ti,ter are the initial and terminal temperature of building i .

4.1. MPC-based rolling optimization In the stage of rolling optimization, the energy management system first updates the predictive data of outdoor temperature, solar radiation and loads based on the weather and history information. Then, multiinterval equipment scheduling and demand management plans are generated based on the predicted values, including the ON-OFF states, operating conditions, supply power of chillers, operation conditions and power of the thermal storage devices, supply water temperature of the CES and the injected powers of buildings. However, the system only

3.4. Coupling of the sources, networks and loads In the cooling season, central energy station chills the high-temperature water in the water collector, and then the chilled water is shunted to each pipeline by the water distributor. Thus, the energy transfer is realized. The coupling relationship between the CES and the pipelines is shown as follows: 256

Applied Energy 243 (2019) 250–265

C. Lv, et al.

executes the schedule of the first interval. In each rolling period, multiple devices are coordinated to minimize the operation cost and unit startup/shutdown penalty. NT

minF = t = tS

CtP PtTL t +

NT i S

j

Uti, j (1

Uti

1, j ) Ei

Incorporating the rolling optimization method, we formulate the robust scheduling plan for adapting to the outdoor temperature uncertainty, through adjustments of the lower/upper limits and the ramping rates of the indoor temperature on the basis of uncertain information in the control domain. Each rolling optimization should determine not only the predictive data of the prediction domain, but also the outdoor temperature deviation amplitude of the control domain relative to the rolling prediction value. The deviation amplitude of the outdoor temperature in the . Therefore, the indoor temcontrol domain can be denoted as Mtout S perature’s deviation amplitude of building i can be calculated by Eq. (61).

(60)

i t = tS

where tS is the first interval of the rolling optimal scheduling; denotes the purchase price of electricity at time t; PtTL denotes the tie-line power (purchasing power from the external grid) at time t; Ei denotes the startup penalty of different types of units; and NT represents the total interval of the rolling period. Because all devices in the system are electrically driven, the operation cost is the cost of purchasing electricity from the external grid. The startup/shutdown cost is the total penalty of each unit in each subsystem, where S is {HP, WC, DC}. The purpose of considering the startup penalty of each unit in the objective function is to avoid frequent equipment start/stop. Here, we ignore the shutdown cost due to its lower value. For better adaption to the cyclical changes of TOU tariff, the rolling method with prediction horizon reducing and control horizon unchanged is selected to minimize the cost of each scheduling day. The scheduling plan has the optimality throughout the whole scheduling day. And the overall economy for long-time operation can be realized through the reasonable consideration of different scheduling days. The schematic of the rolling optimization based on MPC is shown in Fig. 4. At the time t 0 + m t , the system scheduling plan of the prediction domain is obtained by the optimization solution based on the prediction domain information (t 0 + m t to t 0 + n t ), and only the scheduling plan in the control domain is executed. Then the prediction domain is compressed to the last scheduling interval, and the control domain is continuously moved backward. This rolling process repeats until the execution scheduling is generated for all scheduling intervals.

CtP

Mtout Ki Fi t S

Mi, t S =

Cair

(61)

air Vi

The indoor temperature upper and lower limits of building i in the control domain can be modified as follows:

Ti, t S = T

T

_ i, tS

(62)

Mi, t S

= T + Mi, tS

(63)

_

Similarly, the indoor temperature ramping limits of building i in the control domain can be modified as follows:

Ti, tS,upper = T

Mi, tS

(64)

Ti, tS,lower = T

Mi, tS

(65)

The adjustment schematics of the lower/upper and the ramping limits of the indoor temperature are presented in Fig. 5. The compact form of the robust scheduling model based on MPC can be written as:

min F (60) s. t. (1) (29), (31) (41) (43) (59), (61) (65)

4.2. Robust scheduling of the control domain In the actual operation process, due to the existence of prediction error, the outdoor temperature, solar radiation, and electrical load will change in every scheduling interval. For electricity, the system can automatically balance the fluctuations of solar radiation and electrical load. However, deviation of the outdoor temperature will directly affect the indoor temperature of buildings, causing problems that seriously affect the energy supply quality, such as too high/low temperatures, and large temperature fluctuations between adjacent times.

(66)

where the constraints include the operational constraints of the CES, (1)–(29); the constraints of DHN, (31)–(41); the constraints of demand response of buildings, (43)–(48); the coupling constraints of the CES and pipelines, (49)–(53); the coupling constraints of the pipelines and the buildings, (54)–(56); the power balance of cooling/electricity, (57)–(59); and the robust constraints of the indoor temperature in the control domain, (61)–(65).

Prediction domain

t0 t0 +∆t

t0 + m∆t

t0 + n∆t

Control domain

t0 t0 +∆t

Prediction domain

Time-window moves back

Control domain

t0 + n∆t

t0 + m∆t Control domain

Prediction domain t0 + n∆t Control domain

t0 + (m + 1)∆t

Prediction domain

t0 t0 +∆t

t0 + (n − 1)∆t t0 + n∆t

t0 + (m + 1)∆t Fig. 4. Rolling optimization based on MPC. 257

Repeat rolling

t0 t0 +∆t

Applied Energy 243 (2019) 250–265

C. Lv, et al.

T

Ti ,ts M i ,t

s

Temperature range

Tset

T i, ts M i ,t

T

4.3. Linearization of the optimization model

M i ,ts

∆T

∆Ti ,ts ,upper

Ti ,ts −1

∆Ti ,ts ,lower

∆T

In the above optimization model, constraints (1), (6) and (12)–(13) contain product terms of a binary variable and a continuous variable, and the objective function (60) contains the product terms of two binary variables when expanded. To reduce the solving difficulty, we linearize the nonlinear terms by the method in [40]. After linearization, the original optimization problem can be transformed into a mixed-integer linear programming (MILP) problem. This MILP model is solved by using the OPTI optimization toolbox for linking ILOG’s CPLEX 12.8 solver [44]. The computation is performed on a PC with an Intel Xeon CPU E5-2603 v3 @1.60 GHz processor and 8 GB RAM.

Ramping range

M i ,ts

s

Fig. 5. Adjustment of lower/upper ramping limits of indoor temperature.

5. Case study This section applies the proposed method to a typical day in the cooling period. The scheduling period is from 23:00 to 22:00 of next day, and the scheduling interval is 30 min, i.e., a whole scheduling cycle contains 48 scheduling intervals. Table 1 shows the time of use (TOU) tariff of the relevant region. The TOU tariff presents the “peakvalley-flat” characteristics and the first 16 intervals are in the valley period of the electricity price. The detailed equipment parameters and energy supply conditions of centralized power station can be found in [40]. The startup penalties of the ground source heat pump, conventional water-cooled chiller and double-duty chiller units are set as 40.0 RMB, 120.0 RMB and 120.0 RMB per time, respectively. The upper and lower limits of the indoor temperature are 25 °C and 19 °C, respectively, and the standard room temperature is 22 °C; the upper limit of indoor temperature variation between adjacent dispatch intervals is 2 °C; the upper and lower limits of the average indoor temperature are 23 °C and 21 °C, respectively [12,32]. The specific heat capacity and density of air are 1.007 kJ/(kg °C) and 1.2 kg/m3, respectively. The specific heat capacity and density of water are 4.2 kJ/(kg °C) and 1000 kg/m3, respectively. The upper/lower limits of the supplied temperature are 8 °C and 5 °C. The heat loss coefficient of the pipelines is set as 10 kW/(m2 °C). The structure of the system’s heating network is shown in Fig. 7. The network contains 13 pipelines, 9 load nodes (buildings) and 5 intersection nodes. The control modes of heating networks can be quality regulation (the mass flow rates are constant; the water temperatures are variables) or quality-quantity regulation (the mass flow rates and the water temperatures are both variables) [45], and quality regulation is widely used in northern China. Likewise, the CIES adopts the quality regulation mode. The parameters of the heating network and the buildings are shown in Tables A1 and A2 of the Appendix, respectively. In the rolling optimization stage, the multi-period prediction information of outdoor temperature, solar radiation, electricity load, and the deviation amplitude of outdoor temperature relative to the rolling prediction value in the control domain can be obtained by various methods [46]. This paper focuses on the effectiveness of the scheduling strategy, and the prediction method is not discussed. It is assumed that the real-time values of the solar radiation and electricity load are obtained by superimposing the errors of normal distribution on the rolling prediction information (both expectations are 0; the standard deviations of the radiation and load are set as 6% and 3% respectively). The deviation amplitude of the outdoor temperature relative to the rolling prediction value is 1.5%; the real-time data is simulated by superimposing the errors of uniform distribution on the rolling prediction value. The curves of electricity load, solar radiation, and outdoor temperature in a typical day are shown in Fig. A1 of the Appendix, including the day-ahead forecast, rolling forecast and actual values. To verify the correctness and effectiveness of the proposed operation model and the robust scheduling strategy, the rolling optimization scheduling strategy without robust scheduling model (i.e., without the consideration of real-time errors) and the rolling optimization

Fig. 6. Coordinated control process of the MPC-based robust scheduling.

Table 1 TOU tariff of the CIES. Category

Period

Electricity price/(RMB/kWh)

Peak Valley Flat

8:00–11:00, 18:00–23:00 00:00–7:00, 23:00–00:00 7:00–8:00, 11:00–18:00

1.35 0.47 0.89

In engineering practice, the application and coordination of the proposed scheduling strategy are depicted in Fig. 6. Based on the statues and forecasting data of CIES, energy management system delivers the corresponding schedules of control domain to each device for execution. After responding to the actual RES, loads and environment information, the operation statuses of the physical system will be acquired for the next rolling optimization, which plays a role in feedback correction. The proposed scheduling framework provides a novel operation flexibility of CIES by unified source-network-load scheduling and feedback strategy. Based on rolling optimization and advance reserving means, the strategy regulates multiple flexible and controllable links to facilitate the system operation, including the economic performance and adaptability to the uncertainties. 258

Applied Energy 243 (2019) 250–265

C. Lv, et al.

scheduling strategy including the robust scheduling model are adopted for analysis.

N4

5.1. MPC-based rolling optimization

P3 180/273

In this scenario, it is assumed that the rolling forecasting information is accurate, so there is no need to perform robust scheduling of the execution interval. We use the rolling forecasting values of electricity load, solar radiation and outdoor temperature shown in Fig. A1 for calculation.

P2 95/273 N2

N3

N5

Building 2

Building 1

Building 3

Building 4 CES N1

P1 280/426

5.1.1. Optimal operation considering source-network-load characteristics Fig. 8 shows the cooling power allocation of the CES. In the valley period (23:00–07:00), the cooling demand is mainly supplied by the ground source heat pumps (double-duty chillers are in ice-making mode). In the non-valley period, the ice-storage tank and ground source heat pumps are the main cooling suppliers, and the conventional watercooled chillers serve as the complementary suppliers due to the lower COP. Throughout the whole period, the electricity demand of the CIES is satisfied by the PV and the external grid. When the PV output is fully absorbed, the insufficient electricity demand will be balanced by the external grid. Fig. 9 shows the curve of the cooling energy stored in the ice-storage tank. The double-duty chillers are in ice-making mode during most of the valley period. During this period, the energy stored in the cooling storage equipment increases gradually. In the non-valley period, the energy stored in the ice-storage tank is discharged completely to meet the cooling demand and there is no energy redundancy. The chargedischarge characteristics of the thermal storage device and the “peakvalley-flat” characteristic of the TOU tariff are well incorporated by the scheduling strategy for the goal of less operation cost. The supply/return water temperature of the CES is shown in Fig. 10. The supply water temperature varies between 5 and 8 °C according to the schedule, and the return water temperature rises due to the heat exchange of the pipelines with ambient soil and the power consumption of buildings on the demand side. Meanwhile, due to the influence of the pipeline transfer delay, the return temperature curve has an overall rightward trend relative to the supply curve. Combining Fig. 10 with Fig. 8, it can be seen that the cooling power of the CES is directly proportional to the difference between the supply and return water temperatures. Fig. 11 depicts the variation of the indoor temperature in each building. The temperature change trends of multiple buildings are similar overall, and there are slight shifts on the time axis duo to delay differences. In the valley period, the overall building temperatures are relatively low, and plenty of cooling energy is stored in the building storage system. In the non-valley period, the indoor temperatures are higher than those in the valley period. The cooling energy stored during the valley period is released and the overall stored cooling energy decreases. The DSM method with the cooling energy transfer between the different periods adjusts the indoor temperature in the comfort range (19–25 °C) rather than remaining it at the standard temperature (22 °C) and can obtain a better economy incorporating the TOU tariff. The power stored in the pipelines and buildings can be calculated by the output power difference of the CES between scenario 4 and scenario 1 (the scenario adopting the scheduling strategy of this paper and the scenario adopting the conventional economy scheduling strategy, respectively; see Section 5.1.2), as shown in Fig. 12. It can be seen that the pipelines and buildings are in the charge/discharge state in the valley/non-valley period. Combined with the TOU tariff, the operation cost can be significantly reduced by utilizing the characteristics of “valley to charge and non-valley to discharge” in the pipelines and buildings.

P4 600/325

Building 9

P13 100/219

P5 250/325

N7

N8

N6

P7 100/159

P6 110/219

P8 60/273

P9 100/159

N9

N10

Building 5

Building 6

P10 90/273

Building 8

N12

N11

P11 110/219

P12 130/159 N13

Pipeline

Intersection node

Building 7 Load node

Fig. 7. Structure of the heating network in the CIES (P denotes pipeline; N denotes node. The numbers attached to the pipelines are lengths (meter) and diameters (millimeter)).

4000 3000 2000 1000 0

23:00 00:00 01:00 02:00 03:00 04:00 05:00 06:00 07:00 08:00 09:00 10:00 11:00 12:00 13:00 14:00 15:00 16:00 17:00 18:00 19:00 20:00 21:00 22:00

Power (kW)

5000

Time (h) Ground source heat pump

Conventional water-cooled chiller

Double-duty chiller

Ice-storage tank

Central energy station

Fig. 8. Cooling power allocation of the CES.

40000 35000

4000

0 -2000

25000 23:00 0:00 1:00 2:00 3:00 4:00 5:00 6:00 7:00 8:00 9:00 10:00 11:00 12:00 13:00 14:00 15:00 16:00 17:00 18:00 19:00 20:00 21:00 22:00

Power (kW)

30000 2000

20000 15000 10000

-4000 -6000 Ice-making power of DCs

5.1.2. Analysis of different scenarios To analyse the influences of the dynamic transfer characteristic of

5000

Time (h) Releasing power of IT

0 Stored energy in IT

Fig. 9. Curve of the cooling energy stored in the ice-storage tank. 259

Stored energy in IT (kWh)

6000

Applied Energy 243 (2019) 250–265

C. Lv, et al.

DHN and the energy storage characteristic of buildings on the scheduling of the CIES, four scenarios are adopted for comparative analysis:

14

Temperature ( )

12 10

(1) Scenario 1: conventional economy scheduling, i.e., without the consideration of the pipelines’ dynamic characteristic and the buildings’ storage characteristic; (2) Scenario 2: the scheduling model only considers the buildings’ storage characteristic; (3) Scenario 3: the scheduling model only considers the pipelines’ dynamic characteristic; (4) Scenario 4: the scheduling model simultaneously considers both the pipelines’ dynamic characteristic and the buildings’ storage characteristic.

8 6 4

0

23:00 00:00 01:00 02:00 03:00 04:00 05:00 06:00 07:00 08:00 09:00 10:00 11:00 12:00 13:00 14:00 15:00 16:00 17:00 18:00 19:00 20:00 21:00 22:00

2

Time (h) Supply temperature

Return temperature

The output of the CES and the load power from different scenarios are shown in Figs. 13 and 14. When the energy storage characteristics of the building are not taken into account (Scenarios 1 and 3), the indoor temperature is kept at standard temperature of 22 °C. The load powers of the buildings are related only to the outdoor temperature, and the load powers under the two scenarios are equal. When the storage characteristics are considered (Scenarios 2 and 4), the load demand is adjusted according to the environmental temperature and electricity price information. The buildings will charge in the valley period and discharge in the non-valley period for a better economy performance. The comparison of Figs. 13 and 14 shows that the output of the CES and the load power are equal when the transfer effect of the pipelines is not considered (Scenarios 1 and 2). After considering the transfer effect (Scenarios 3 and 4), the supply and demand are not balanced due to the influences of heat loss and delay; the output of the CES at a specific time passes through the pipelines with energy loss and is matched with the load power after a certain time delay. Table 2 compares the system operation cost of different scheduling scenarios. In contrast with Scenario 1, Scenario 2 has considered the building storage characteristic. Different from the conventional scheduling method of keeping the temperature at the standard value, the indoor temperature can be adjusted within the comfortable range to improve the energy system’s optimization space in Scenario 2. The operation cost is reduced by 4.9% with buildings charging in the valley period and discharging in non-valley period. The comparison between Scenarios 1 and 3 shows that the operation cost increases in Scenario 3 due to the heat loss of the heating network. However, the scheduling plan is more practical because it considers the dynamic transfer process which exists in actual operation. In Scenario 4, the network transfer characteristics and building storage characteristics are considered together, the operation cost decrease by 2.7%. The scheduling flexibility of Scenario 4 increases compared with that of Scenario 3 when the building storage characteristic is considered, and the operation cost is 2744.9 RMB lower. To analyse the influence of transfer delay of the pipelines, two scenarios are added:

27 25 23 21 19 17 15

23:00 00:00 01:00 02:00 03:00 04:00 05:00 06:00 07:00 08:00 09:00 10:00 11:00 12:00 13:00 14:00 15:00 16:00 17:00 18:00 19:00 20:00 21:00 22:00

Indoor temperature ( )

Fig. 10. Supply/return water temperature curve of the CES.

Time (h) Building 1

Building 2

Building 3

Building 4

Building 5

Building 6

Building 7

Building 8

Building 9

Fig. 11. Indoor temperature curves of the buildings.

6000 5000

3000 2000 1000 0 -1000

23:00 00:00 01:00 02:00 03:00 04:00 05:00 06:00 07:00 08:00 09:00 10:00 11:00 12:00 13:00 14:00 15:00 16:00 17:00 18:00 19:00 20:00 21:00 22:00

Power (kW)

4000

-2000 -3000

Time (h) CES output-scenario 1

CES output-scenario 4

Storage power

Fig. 12. Storage curves of the pipelines and buildings.

6000

Scenario 5: the scheduling model only considers the pipelines’ heat loss; Scenario 6: the scheduling model simultaneously considers both the pipelines’ heat loss and the buildings’ storage characteristic.

4000 3000 2000

Table 3 shows the results without considering the transfer delay of pipelines. When the influence of the time delay is not considered, the output of the CIES is equal to the sum of the network loss and the consumption of the building loads. In this case, the heating network only plays the role of energy transfer. When the time delay is considered, the supply and demand do not need to be synchronous and there is an imbalance between the injected power, power loss of the network and the building consumption. At this time, the network has a certain energy storage capacity that provides more flexibilities for the

1000 0

23:00 0:00 1:00 2:00 3:00 4:00 5:00 6:00 7:00 8:00 9:00 10:00 11:00 12:00 13:00 14:00 15:00 16:00 17:00 18:00 19:00 20:00 21:00 22:00

Power (kW)

5000

Time (h) Scenario 1

Scenario 2

Scenario 3

Scenario 4

Fig. 13. The output of the CES from different scenarios. 260

Applied Energy 243 (2019) 250–265

5000 4500 4000 3500 3000 2500 2000 1500 1000 500 0

operation optimization of the CES. Comparing the operation cost of Scenario 3 in Table 2 with that of Scenario 5, as well as Scenario 4 with Scenario 6, the operation costs decease by 352.9 RMB and 254.0 RMB, respectively. 5.2. Analysis of MPC-based robust scheduling In each control domain, the actual outdoor temperature may differ from the forecasting value, causing a series of comfort problems on indoor temperature. Therefore, it is necessary to consider the real-time forecast error for a better scheduling.

23:00 0:00 1:00 2:00 3:00 4:00 5:00 6:00 7:00 8:00 9:00 10:00 11:00 12:00 13:00 14:00 15:00 16:00 17:00 18:00 19:00 20:00 21:00 22:00

Time (h) Scenario 4

5.2.1. Results of MPC-based robust scheduling When the robust scheduling model is included, the scheduling strategy can generate robust scheduling schemes that can adapt to uncertain changes according to the deviation limit of outdoor temperature in the control domain. The outdoor temperature in the typical day is shown in Fig. 15, including the upper/lower limits, rolling forecast value, and actual value. Building 5 is chosen to illustrate the schemes’ adaptability to the uncertainty of the environment. The modified upper/lower and ramping limits of the indoor temperature are depicted in Figs. 16 and 17, respectively. In Fig. 16, the indoor temperature’s upper limit is lowered according to the deviation amplitude for adapting the uncertainty of upward deviation; the lower limit is adjusted upward according to the deviation amplitude for adapting the uncertainty of downward deviation. Likewise, the upper ramping limit decreases with upward uncertainty and the lower ramping limit increases with downward uncertainty incorporating the possible deviation range, as shown in Fig. 17. Fig. 18 depicts the actual indoor temperature curve of Building 5. With considering the real-time errors with forecasting values, the rolling temperature schemes vary in the modified temperature range,

Fig. 14. Load power from different scenarios.

Table 2 Operation cost comparison of different scenarios (√ - with — - without). Scenario

Building storage

Pipeline transfer

Operation cost/ RMB

Proportion of reduction

1 2 3 4

— √ — √

— — √ √

56769.8 53983.7 57973.6 55228.7

– 4.9% −2.1% 2.7%

Table 3 Cost without the transfer delay of the pipelines (√ - with — - without). Scenario

Heat loss

Transfer delay

Building storage

Operation cost/RMB

5 6

√ √

— —

— √

58326.5 55482.7

36

Ramping rate ( )

32 30 28 26 24

Rolling forecast

Time (h)

Upper limit

3

2

2

1

1

0

0

-1

-3

Lower limit

22

0

21

-1

20

-2

19

-3

18

-4

Indoor temperature ( )

1

Deviation limit ( )

2

23

23:00 00:00 01:00 02:00 03:00 04:00 05:00 06:00 07:00 08:00 09:00 10:00 11:00 12:00 13:00 14:00 15:00 16:00 17:00 18:00 19:00 20:00 21:00 22:00

3

24

Time (h) Lower limit

Upper deviation limit

Lower deviation limit

-3

Time (h) Upper ramping limit

Lower ramping limit

Upper deviation limit

Lower deviation limit

Fig. 17. The modified ramping limits of the indoor temperature in Building 5.

25

Upper limit

-1 -2

Actual value

Fig. 15. Curves of outdoor temperature in the typical day.

Indoor temperature ( )

3

-2

23:00 0:00 1:00 2:00 3:00 4:00 5:00 6:00 7:00 8:00 9:00 10:00 11:00 12:00 13:00 14:00 15:00 16:00 17:00 18:00 19:00 20:00 21:00 22:00

Temperature ( )

34

Deviation limit ( )

Scenario 3

27

3

25

2 1

23

0

21

-1

19

-2

17

-3

15

-4

Time (h)

Fig. 16. The modified upper/lower limits of the indoor temperature in Building 5.

Actual value

Upper limit

Rolling value

Deviation value

Lower limit

Fig. 18. Curve of actual indoor temperature of Building 5. 261

Deviation value ( )

Scenario 2

23:00 00:00 01:00 02:00 03:00 04:00 05:00 06:00 07:00 08:00 09:00 10:00 11:00 12:00 13:00 14:00 15:00 16:00 17:00 18:00 19:00 20:00 21:00 22:00

Scenario 1

23:00 00:00 01:00 02:00 03:00 04:00 05:00 06:00 07:00 08:00 09:00 10:00 11:00 12:00 13:00 14:00 15:00 16:00 17:00 18:00 19:00 20:00 21:00 22:00

Power (kW)

C. Lv, et al.

Applied Energy 243 (2019) 250–265

0

-1

-1 -2

-2 -3

-3

Time (h) Actual ramping rate

Upper raming limit

Rollling ramping rate

Deviation value

Indoor temperature ( )

Indoor temperature ( )

21 19 17 4

5

6

7

8

3

5

6

7

8

9

Building Minimum ramping rate

27

3

25

2 1

23

0

21

-1

19

-2

17

-3

15

-4

9

Time (h)

Building Minimum temperature

Fig. 20. The maximum/minimum indoor temperature without robust scheduling.

Actual value

Upper limit

Rolling value

Deviation value

Lower limit

Fig. 22. Actual indoor temperature curve of Building 5 without robust scheduling.

Ramping rate ( )

with lower values in the valley period and higher values in the nonvalley period for economy improvement. In actual operation, the indoor temperature will differ from the rolling value for responding to the deviation of outdoor temperature with the cooling power unchanged, and the variations are limited within the allowable margin generated by the robust constraints. Therefore, the indoor temperature always stays within the comfort range. The actual ramping rate of the indoor temperature of Building 5 is shown in Fig. 19. Similarly, after taking into account the uncertain change interval of outdoor temperature in the control domain, the rolling ramping rate of generated temperature is schemed to change within the modified ramping range. When responding to the upward/ downward change of outdoor temperature, the actual ramping rate will change upward/downward within the allowable ramping margin by superimposing the indoor temperature’s deviation values on the rolling values; and there is no rapid temperature fluctuations (greater than 2 °C) between adjacent intervals. Incorporating the result of Fig. 18, it is shown that the uncertainty of the outdoor temperature is well adapted by the advanced robust method, no out-of-limit appears under the uncertain situation and the consumer comfort is maintained.

3

3

2

2

1

1

0

0

-1

23:00 00:00 01:00 02:00 03:00 04:00 05:00 06:00 07:00 08:00 09:00 10:00 11:00 12:00 13:00 14:00 15:00 16:00 17:00 18:00 19:00 20:00 21:00 22:00

Maximum temperature

4

Fig. 21. The maximum/minimum ramping rate without robust scheduling.

23

3

2

Maximum ramping rate

25

2

1

-1

-3

27

1

0

-2

Lower ramiping limit

Fig. 19. Actual ramping rate of the indoor temperature of Building 5.

15

1

Deviation value ( )

0

2

-1 -2

-2 -3

Time (h) Actual ramping rate

Upper raming limit

Rollling ramping rate

Deviation value

Deviation value ( )

1

23:00 00:00 01:00 02:00 03:00 04:00 05:00 06:00 07:00 08:00 09:00 10:00 11:00 12:00 13:00 14:00 15:00 16:00 17:00 18:00 19:00 20:00 21:00 22:00

2

1

)

2

3

Ramping rate(

3

Deviation value ( )

3

23:00 00:00 01:00 02:00 03:00 04:00 05:00 06:00 07:00 08:00 09:00 10:00 11:00 12:00 13:00 14:00 15:00 16:00 17:00 18:00 19:00 20:00 21:00 22:00

Ramping rate( )

C. Lv, et al.

-3

Lower ramiping limit

Fig. 23. Actual ramping rate of the indoor temperature of Building 5 without robust scheduling.

problems of too high/low temperature and ramping rate, which will heavily affect the comfort level. Building 5 is chosen for specific illustration, and the indoor temperature and ramping rate are shown in Figs. 22 and 23, respectively. As shown in Fig. 22, the times at which the real-time temperature goes below the lower limit often occur in the valley period. In this period, the temperature is kept lower for storing cooling energy and no adjustable margins are reserved for the downward forecasting errors; thus, the indoor temperature will exceed the lower limit if the temperature deviates downward. Similarly, the indoor temperature is maintained at a high value during the non-valley period and no adjustable margins are reserved for the upward forecasting errors, and this situation will cause the temperature to go beyond the upper limit if the actual value deviates upward. For the ramping rate, when the values of the rolling optimization are close to the upper/lower limit with no redundancy for adapting to the uncertainty, the temperature will rise/fall

5.2.2. Actual operation without robust scheduling The actual operation without considering the real-time forecasting errors is applied to further demonstrate the effectiveness of the MPCbased robust scheduling strategy. Under uncertain conditions, in which the system performs rolling optimization without robust strategy, the maximum and minimum values and ramping rates of the buildings’ indoor temperatures during a scheduling cycle are shown in Figs. 20 and 21, respectively. The simplex rolling optimization does not carry out the correction of the temperature upper/lower and ramping limits for enhancing the robustness to the forecasting errors, and the indoor temperature may deviate from the comfort range and fluctuate rapidly after responding to the change of outdoor temperature. It appears that almost every building has the 262

Applied Energy 243 (2019) 250–265

C. Lv, et al.

too quickly if the uncertain values deviate upward/downward, which can be seen in Fig. 23. By comparing the results of the robust scheduling model included or not, it is obvious that the robust scheduling strategy based on MPC can effectively adapt to the forecasting errors of outdoor temperature by rolling optimization and robust constraints generated for adjustable margins, achieving better DSM and satisfying consumers' demand comfortably and reliably.

of source-network-load and the forecasting information of outdoor temperature, solar radiation and electricity load in each rolling step, effectively reducing the influence of inaccurate forecasting. The robust scheduling scheme considering the uncertain interval of real-time outdoor temperature in the control domain has stronger robustness and can adapt well to uncertain changes of environmental factor, with maintaining the indoor temperature in the comfortable range and preventing fast fluctuations. In conclusion, the model predictive control based robust strategy can achieve optimal scheduling and coordination of multiple energy links, as well as demand management on the consumer side. Furthermore, the uncertainties of solar radiation, electricity load, and environment factor can be well adapted. For a whole scheduling cycle, the numbers of total variables and integer variables are 9545 and 1577, and the solving time is about 560 s and can be suitable for the scheduling with half an hour’s interval. Note that the computational efficiency will be challenged with the increasing system scale and the scheduling’s timeliness for real application will not be maintained. Thus, an efficient solution for large-scale communities should be developed to reduce problem complexity and number of variables, such as the methods based on benders decomposition [47], decentralized coordination [48] and load equivalent [49].

6. Conclusions This paper presents a model predictive control based robust scheduling strategy that addresses the inaccuracy of forecasts for community integrated system. On the basis of building the dynamic transfer model of pipelines and storage model of buildings, a rolling robust scheduling strategy based on model predictive control is proposed for adapting to uncertainties, which takes into account the flexibilities of source-network-load. Then the operation data of a typical cooling day is chosen to conduct a case study. The proposed scheduling strategy fully considers the coupling of the central energy station, district heating network and building loads. The power of multiple energy equipment is reasonably allocated with the energy storage device charging in the valley period and discharging in the non-valley period. The operation flexibility is enhanced by considering the pipeline delay and adjustable indoor temperature of buildings, which make the operation more economical than the methods that consider single or no flexible characteristics of the district heating network and buildings. The model predictive control based rolling robust scheduling strategy can generate more reasonable schemes by updating the states

Acknowledgments This work was financially supported by the National Key Research and Development Program of China (No. 2018YFB0905000) and conducted in cooperation of APPLIED ENERGY UNiLAB-DEM: Distributed Energy & Microgrid. UNiLAB is an international virtual lab of collective intelligence in Applied Energy.

Appendix A See Tables A1 and A2 and Fig. A1.

Table A1 Parameters of the pipelines. Pipeline

Length/m

Inner diameter/mm

Mass flow/(kg/s)

Pipeline

Length/m

Inner diameter/mm

Mass flow/(kg/s)

1 2 3 4 5 6 7

280 95 180 600 250 110 100

426 273 273 325 325 219 159

77 21 56 42 53.2 17.5 9.1

8 9 10 11 12 13

60 100 90 110 130 100

273 159 273 219 159 219

26.6 9.1 17.5 11.9 5.6 37.1

Table A2 Parameters of the buildings. Building

Surface area/m2

Volume/m3

Thermal conductivity/(W/(m2 °C))

Building

Surface area/m2

Volume/m3

Thermal conductivity /(W/(m2 °C))

1 2 3 4 5

20,000 80,000 40,000 50,000 20,000

40,000 63,000 70,000 135,000 25,500

2 2 2 2 1.2

6 7 8 9

15,000 10,000 20,000 25,000

24,000 22,400 21,500 35,000

1.2 1.2 1.2 1.2

263

Applied Energy 243 (2019) 250–265

C. Lv, et al. 3000

Power (kW)

2500 2000 1500 1000 Day-ahead forecast

Rolling forecast

Actual value

500

23:00 00:00 01:00 02:00 03:00 04:00 05:00 06:00 07:00 08:00 09:00 10:00 11:00 12:00 13:00 14:00 15:00 16:00 17:00 18:00 19:00 20:00 21:00 22:00

0

Time (h) 0.9

Solar radiation/(W/m2)

0.8

Day-ahead forecast

Rolling forecast

Actual value

0.7 0.6 0.5 0.4 0.3 0.2 0.1 23:00 00:00 01:00 02:00 03:00 04:00 05:00 06:00 07:00 08:00 09:00 10:00 11:00 12:00 13:00 14:00 15:00 16:00 17:00 18:00 19:00 20:00 21:00 22:00

0

Time (h)

40

Day-ahead forecast

Rolling forecast

Actual value

Temperature ( )

35 30 25 20

23:00 00:00 01:00 02:00 03:00 04:00 05:00 06:00 07:00 08:00 09:00 10:00 11:00 12:00 13:00 14:00 15:00 16:00 17:00 18:00 19:00 20:00 21:00 22:00

15

Time (h)

Fig. A1. Curves of the load, solar radiation and outdoor temperature in the typical day.

References [10]

[1] Rifkin J. The third industrial revolution: how lateral power is transforming energy, the economy, and the world. Macmillan; 2011. [2] Shen J, Luo C. Overall review of renewable energy subsidy policies in China – contradictions of intentions and effects. Renew Sustain Energy Rev 2015;41:1478–88. [3] Cao Y, Li Q, Tan Y, Li Y. A comprehensive review of Energy Internet: basic concept, operation and planning methods, and research prospects. J Mod Power Syst Clean Energy 2018;6(3):1–13. [4] Guo L, Liu W, Cai J, Hong B, Wang C. A two-stage optimal planning and design method for combined cooling, heat and power microgrid system. Energy Convers Manage 2013;74(10):433–45. [5] Yu H, Zhang C, Deng Z, Bian H. Economic optimization for configuration and sizing of micro integrated energy system. J Mod Power Syst Clean Energy 2017;6(2):1–12. [6] Wu Q, Zheng J, Jing Z. Coordinated scheduling of energy resources for distributed DHCs in an integrated energy grid. CSEE J Power Energy Syst 2015;1(1):95–103. [7] Cho H, Smith AD, Mago P. Combined cooling, heating and power: a review of performance improvement and optimization. Appl Energy 2014;136:168–85. [8] Koirala BP, Koliou E, Friege J, Hakvoort RA, Herder PM. Energetic communities for community energy: a review of key issues and trends shaping integrated community energy systems. Renew Sustain Energy Rev 2016;56:722–44. [9] Zhang Y, Campana PE, Yang Y, Stridh B, Lundblad A, Yan J. Energy flexibility from

[11] [12] [13] [14] [15] [16] [17]

264

the consumer: integrating local electricity and heat supplies in a building. Appl Energy 2018;223:430–42. Li G, Zhang R, Jiang T, Chen H, Bai L, Cui H, et al. Optimal dispatch strategy for integrated energy systems with CCHP and wind power. Appl Energy 2017;192:408–19. Zhou Y, Hu W, Min Y, Dai Y. Integrated power and heat dispatch considering available reverse of combined heat and power units. IEEE Trans Sustain Energy 2018. Cai H, Ziras C, You S, Li R, Honore K, Bindner HW. Demand side management in urban district heating networks. Appl Energy 2018;230:506–18. Ruan Y, Liu Q, Li Z, Wu J. Optimization and analysis of Building Combined Cooling, Heating and Power (BCHP) plants with chilled ice thermal storage system. Appl Energy 2016;179:738–54. Lu Y, Wang S, Sun Y, Yan C. Optimal scheduling of buildings with energy generation and thermal energy storage under dynamic electricity pricing using mixedinteger nonlinear programming. Appl Energy 2015;147:49–58. Li G, Zhang R, Jiang T, Chen H, Bai L, Li X. Security-constrained bi-level economic dispatch model for integrated natural gas and electricity systems considering wind power and power-to-gas process. Appl Energy 2017;194:696–704. Renaldi R, Kiprakis A, Friedrich D. An optimization framework for thermal energy storage integration in a residential heat pump heating system. Appl Energy 2017;186:520–9. Oluleye G, Allison J, Hawker G, Kelly N, Hawkes AD. A two-step optimization model for quantifying the flexibility potential of power-to-heat systems in

Applied Energy 243 (2019) 250–265

C. Lv, et al. dwellings. Appl Energy 2018;228:215–28. [18] Stinner S, Huchtemann K, Müller D. Quantifying the operational flexibility of building energy systems with thermal energy storages. Appl Energy 2016;181:140–54. [19] Jin X, Mu Y, Jia H, Wu J, Xu X, Yu X. Optimal day-ahead scheduling of integrated urban energy systems. Appl Energy 2016;180:1–13. [20] Ji H, Wang C, Li P, Song G, Yu H, Wu J. Quantified analysis method for operational flexibility of active distribution networks with high penetration of distributed generators. Appl Energy 2019;239:706–14. [21] Lu S, Gu W, Zhou J, Zhang X, Wu C. Coordinated dispatch of multi-energy system with district heating network: modeling and solution strategy. Energy 2018;152:358–70. [22] Nuytten T, Claessens B, Paredis K, Bael JV, Six D. Flexibility of a combined heat and power system with thermal energy storage for district heating. Appl Energy 2013;104:583–91. [23] Zheng J, Zhou Z, Zhao J, Wang J. Integrated heat and power dispatch truly utilizing thermal inertia of district heating network for wind power integration. Appl Energy 2018;211:865–74. [24] Li Z, Wu W, Shahidehpour M, Wang J. Combined heat and power dispatch considering pipeline energy storage of district heating network. IEEE Trans Sustain Energy 2016;7(1):12–22. [25] Fang J, Zeng Q, Ai X, Chen Z. Dynamic optimal energy flow in the integrated natural gas and electrical power systems. IEEE Trans Sustain Energy 2018;9(1):188–97. [26] Clegg S, Mancarella P. Integrated electrical and gas network flexibility assessment in low-carbon multi-energy systems. IEEE Trans Sustain Energy 2016;7(2):718–31. [27] Xue X, Wang S, Sun Y, Xiao F. An interactive building power demand management strategy for facilitating smart grid optimization. Appl Energy 2014;116(3):297–310. [28] Esther BP, Kumar KS. A survey on residential demand side management architecture, approaches, optimization models and methods. Renew Sustain Energy Rev 2016;59:342–51. [29] Parra D, Norman SA, Walker GS, Gillott M. Optimum community energy storage for renewable energy and demand load management. Appl Energy 2017;200:358–69. [30] Dong X, Mu Y, Xu X, Jia H, Wu J, Yu X. A charging pricing strategy of electric vehicle fast charging stations for the voltage control of electricity distribution networks. Appl Energy 2018;225:857–68. [31] Yang Y, Wu K, Long H, Gao J, Yan X, Kato T, et al. Integrated electricity and heating demand-side management for wind power integration in China. Energy 2014;78:235–46. [32] Jin X, Mu Y, Jia H, Wu J, Jiang T, Yu X. Dynamic economic dispatch of a hybrid energy microgrid considering building based virtual energy storage system. Appl Energy 2016;194:386–98. [33] Najafi A, Arsalan H, Contreras J, Ramezani M. Medium-term energy hub management subject to electricity price and wind uncertainty. Appl Energy

2016;168:418–33. [34] Bai L, Li F, Cui H, Jiang T, Sun H, Zhu J. Interval optimization based operating strategy for gas-electricity integrated energy systems considering demand response and wind uncertainty. Appl Energy 2016;167:270–9. [35] Zhao Y, Lu Y, Yan C, Wang S. MPC-based optimal scheduling of grid-connected low energy buildings with thermal energy storages. Energy Build 2015;86:415–26. [36] Zhang Y, Meng F, Wang R, Zhu W, Zeng X. A stochastic MPC based approach to integrated energy management in microgrids. Sustain Cities Soc 2018;41:349–62. [37] Holjevac N, Capuder T, Kuzle I. Adaptive control for evaluation of flexibility benefits in microgrid systems. Energy 2015;92:487–504. [38] Zhang F, Jin X, Mu Y, Jia H, Yu X, Liu D. Model predictive scheduling method for a building microgrid considering virtual storage system. Proc CSEE 2018;38(15):4420–8. [39] Luo Z, Wu Z, Li Z, Cai H, Li B, Gu W. A two-stage optimization and control for CCHP microgrid energy management. Appl Therm Eng 2017;125:513–22. [40] Wang C, Lv C, Li P, Song G, Li S, Xu X, et al. Modeling and optimal operation of community integrated energy systems: a case study from China. Appl Energy 2018;230:1242–54. [41] Yokoyama R, Kitano H, Wakui T. Optimal operation of heat supply systems with piping network. Energy 2017;137:888–97. [42] Gu W, Wang J, Lu S, Luo Z, Wu C. Optimal operation for integrated energy system considering thermal inertia of district heating network and buildings. Appl Energy 2017;199:234–46. [43] Cristofari C, Norvaišienė R, Canaletti JL, Notton G. Innovative alternative solar thermal solutions for housing in conservation-area sites listed as national heritage assets. Energy Build 2015;89:123–31. [44] IBM. IBM ILOG CPLEX optimization studio. Academic Initiative ed. Available: < https://www.ibm.com/products/Ilog-cplex-optimizationstudio > [assessed 20.12.18]. [45] Li Z, Wu W, Wang J, Zhang B. Transmission-constrained unit commitment considering combined electricity and district heating networks. IEEE Trans Sustain Energy 2016;7(2):480–92. [46] Raza MQ, Khosravi A. A review on artificial intelligence based load demand forecasting techniques for smart grid and buildings. Renew Sustain Energy Rev 2015;50:1352–72. [47] Siddiqui S, Gabriel SA, Azarm S. Solving mixed-integer robust optimization problems with interval uncertainty using Benders decomposition. J Operat Res Soc 2015;66(4):664–73. [48] Guo Y, Pan M, Fang Y, Khargonekar PP. Decentralized coordination of energy utilization for residential households in the smart grid. IEEE Trans Smart Grid 2013;4(3):1341–50. [49] Wang J, Peng J, Duan J, Liu Z, Yuan S, Chao Q. Control strategy for suppressing power fluctuation of equivalent load in microgrids based on demand-side reservation and response. Autom Electr Power Syst 2017;41(4):69–77.

265