Multi-energy-storage energy management with the robust method for distribution networks

Electrical Power and Energy Systems 118 (2020) 105779

Contents lists available at ScienceDirect

Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes

Multi-energy-storage energy management with the robust method for distribution networks

T

Quan Sui, Fanrong Wei, Xiangning Lin, Chuantao Wu, Zhixun Wang, Zhengtian Li

⁎

State Key Laboratory of Advanced Electromagnetic Engineering and Technology, Huazhong University of Science and Technology, Wuhan, Hubei 430074 China

ARTICLE INFO

ABSTRACT

Keywords: Distribution network Mobile energy storage Thermal energy storage Energy management system Iterative method

The randomness, volatility and anti-peaking characteristic from distributed renewable energy generation rise great challenges for the safe and economic operation of the distribution networks (DNs). To address this problem, this paper proposes a novel multi-energy-storage energy management system (EMS) to co-optimize the electricity-driven mobile energy storage (MES) and inverter air-conditioning (AC)-based thermal energy storage (TES). To facilitate the energy management of the DN, the MES that considers the delay factors and the TES that regulates reactive power have been developed into a unified analytic model capable of charging and discharging. In addition, considering the impact of the forecasting uncertainties and the risk-aversion of the dispatcher, a novel robust optimization method is proposed to obtain more accurate “worst scenario”. The dispatching model is then converted into a mixed integer second-order cone programming problem (MI-SOCP) and a mixed integer linear programming problem (MILP), and linearized techniques and an iteration method are used to efficiently solve these problems. Simulation studies on a 41-node DN in Ontario indicate that the operational cost and power loss of the DN can be reduced by no less than 1% and 8% using the proposed EMS, respectively, while a safer voltage level with a voltage deviation of 5% can be obtained. The results confirm the effectiveness of the MES and TES for peak shaving, valley filling and voltage supporting.

1. Introduction Energy crises and environmental requirements have stimulated the increasing level of penetration of renewable energy sources (RESs) in the distribution network (DN). Meanwhile, the randomness, volatility and anti-peaking characteristic of the RES cause great challenges for the safe and economic operation of the DN, such as worsening the node voltage quality and power flow distribution [1]. Although these problems can be solved well by multiple stationary electric energy storage systems, the application of these multiple stationary electric energy storage systems is limited because of the high operation and maintenance costs [2]. The numerous fixed-frequency air conditioner (AC) integrated buildings can act as this thermal energy storage (TES), but the TES can only be used to carry out peaking cutting [3]. Therefore, it is crucial to find some economical and practical solutions to optimize the economic and reliable operation of the DN. Research shows that compared with stationary electric energy storage systems and fixed-frequency ACs, mobile energy storage systems (MESs) and inverter ACs have more flexible regulation capability [4–12,13–21]. For example, the MES can be utilized to eliminate transmission line congestion [4–6] and accelerate load recovery [7–10]. ⁎

Meanwhile, the inverter AC can be used to conduct harmonic current mitigation [15–16], peaking cutting [19] and frequency modulation [20]. However, essentially, as an electrochemical energy storage system, the MES with a high operational cost can be utilized to carry out energy arbitrage only when the electricity price difference is large. Meanwhile, with the small thermal storage capacity and significant heat loss, the TES cannot even store the energy for tens of minutes. Focusing on separately scheduling the MES or TES, the existing operations must endure their shortcomings, resulting in the unsatisfactory regulation effect. To resolve the above problems, a novel multi-energy storage EMS, together with the modified models of the MES and TES, is proposed in this paper to exploit their complementary advantages and evaluate the co-optimization performance. Several studies have considered MES modelling. Refs. [4,5] propose the MES model integrated time-space network to depict its location and charging/discharging state. Considering the uncertainties of the load and wind [6], presents a scenario-based stochastic model of the MESintegrated power system scheduling. Ref. [7] depicts a set of methodologies to assess the reliability of a power network with the MES. A temporal and spatial MES model that is related to both transportable networks and DNs is presented in [8] to illustrate a post-disaster joint

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

https://doi.org/10.1016/j.ijepes.2019.105779 Received 27 August 2019; Received in revised form 8 December 2019; Accepted 11 December 2019 0142-0615/ © 2019 Elsevier Ltd. All rights reserved.

Electrical Power and Energy Systems 118 (2020) 105779

Q. Sui, et al.

Nomenclature

d UTES,ik, t HTES, i, t Pi, t Qi, t

Abbreviations EMS DN MES AC TES MI-SOCP MILP RES SOC OLTC

energy management system distribution network mobile energy storage air conditioning thermal energy storage mixed integer second-order cone programming mixed integer linear programming renewable energy source state of charge on-load tap changer

i, t

kijs Pij, i, t Qij, i, t Pij,j, t Qij, j, t Lij, t Vj, t ij, t

j, t

P l oss, ij, t Ppurchase, t F ,t u1 , t u2 , t UT,s j, t v js, t

Indices t i, j A Ai + Ai MES AC TES o d max min c d b r l n s T k pdc iter

index for time index for grid node or station set of routes set of routes which start from station i set of routes which end at station i index for MES index for AC Index for TES index for outdoor environment index for indoor environment index for maximum index for minimum index for charging index for discharging index for benchmark index for renewable energy generator index for load index for RES output, load, electricity price and transit delay index for uncertain domain bound index for OLTC’s voltage level index for OLTC index for building index for prediction index for iteration number

Parameters

Tij, t Tij, t T c Pmax d Pmax SMES Rik Ca, ik k1 ,k2 ,l1,l2 To, t Pv Emax MES Emax TES,ik Tmax Tmin P max TES,ik STES,ik Ki P l, i, t Pr, i, t Ql, i, t Qr, i, t rij x ij ¯ij

Variables

Iij, t Uij, t PMES,i,t QMES,i, t c UMES,i, t d UMES,i, t c PMES,i, t d PMES,i, t b PMES, t EMES, t socMES, t WTES,ik,t PTES,ik,t Td,ik,t E TES,ik,t b PTES, ik ,t c PTES,ik, t d PTES,ik, t QTES,ik,t c UTES,ik, t

discharging flag of the TES k at node i at time t variables of the aggregated TESs at node i at time t active power injected into node i reactive power injected into node i voltage of node i transformer ratio of voltages at node i and at node j active power of branch ij from node i reactive power of branch ij from node i active power of branch ij into node i reactive power of branch ij into node i current of branch ij at time t voltage of node j at time t second norm current of branch ij at time t second norm voltage of node j at time t power loss of branch ij at time t power purchased by the DN from the external grid actual RES output, load, electricity price and transit delay flag of lower deviation from prediction flag of upper deviation from prediction flag of voltage level voltage at level s

MES existence flag on route ij at time t value switching flag of Iij, t active power of the MES at station i at time t reactive power of the MES at station i at time t charging flag of the MES at station i at time t discharging flag of the MES at station i at time t charging power of the MES at station i at time t discharging power of the MES at station i at time t driving power of the MES at time t energy storage of the MES at time t SOC of the MES at time t refrigeration power of the TES k at node i at time t electric power of the TES k at node i at time t indoor temperature of building k at node i at time t energy storage of the TES k at node i at time t electric power maintaining indoor temperature at time t charging power of the TES k at node i at time t discharging power of the TES k at node i at time t reactive power of the TES k at node i at time t charging flag of the TES k at node i at time t

j

¯j p1, t Clabor p2, t F m, t in

F max ,t F pdc ,t M U

2

duration of the MES transiting route ij traffic congestion duration of route ij dispatching cycle maximum charging power of the MES maximum discharging power of the MES inverter capacity of the MES thermal resistance of building k at node i thermal capacity of building k at node i Parameters of AC outdoor temperature rated driving power of the MES charging/discharging efficiency energy storage capacity of the MES energy storage capacity of the TES k at node i upper limit of the comfortable temperature range lower limit of the comfortable temperature range maximum active power of the TES k at node i inverter capacity of the TES k at node i at time t number of the TESs at node i active power of common loads at node i at time t active power of RES at node i at time t reactive power of common loads at node i at time t reactive power of the RES at node i at time t resistance of branch ij reactance of branch ij current limit of branch ij voltage lower limit of branch ij voltage upper limit of branch ij electricity purchase price daily labour cost for the MES driver charging/discharging cost of the MES minimum predicted value of the RES output, load, electricity price and transit delay maximum predicted value of the RES output, load, electricity price and transit delay predicted value of the RES output, load, electricity price and transit delay sufficiently large positive value uncertain domain budget of uncertainty

Electrical Power and Energy Systems 118 (2020) 105779

Q. Sui, et al.

restoration scheme. A novel co-optimization model is formulated in [9] to route repair crews and MESs in the transportation network, schedule them in the DN. In addition [10], implements resilient routing and scheduling of MESs via a two-stage framework where MESs are prepositioned to enable rapid pre-restoration before the event and dynamically dispatched to coordinate with conventional restoration efforts after the event. Ref. [11] discusses a two-stage optimization dispatching model for the MES, where the MES designated stations are optimized in the first stage while the transit times are calculated in the second stage. Ref. [12] shows a two-stage optimization model that optimizes investments in MES in the first stage and can re-route the installed MES in the second stage. We can learn a great deal from the above research studies. However, the existing studies neglect the impact of the traffic congestion risk on the MES’s transit time and energy consumption, where dispatching strategies may not actually be supported. Additionally, in the future, the transportable batteries moved by a fuel-driven truck in [4–12] will inevitably be replaced by the electricity-driven MES, where traffic constraints (such as state of charge, SOC) cannot be regarded as irrelevant limits for the DN, but will be reflected directly in the EMS. As a further study, the increasing inverter ACs with multiple functions and service potentials have also recently received considerable attention [13]. The control model of the compressor speed is presented in [14] to achieve optimum inverter AC using the concept of the voltage control at the proportional dc link. A pulse width modulation converter inverter model for building an AC compressor is presented in [15,16] to mitigate the harmonic current generation. Ref. [17] proposes a dynamic model of an experimental room to estimate the effect of the direct load control application to the inverter AC on its indoor air temperature. Ref. [18] develops a data-driven dynamic model of the inverter AC through real-time experimental studies with a time horizon ranging from seconds to hours. Ref. [19] models an inverter AC system as a thermal battery by storing the thermal energy into the buildings to reduce the costs of electricity to the users. Ref. [20] introduces a cluster aggregation model of AC to carry out the demand response such as frequency modulation. Additionally, a consensus control strategy is designed in [21] to solve the inverter AC group operation model. However, with the AC-DC-AC structure, the inverter ACs are capable of reactive power regulation, which has never been discussed in the existing studies. Additionally, the inverter AC integrated buildings in [19] are formulated as a TES model with some exponential equations to ensure the result accuracy, resulting in a huge computation burden. Finally, only the power consumption of the TES in [13–21] is optimized according to the real-time electricity price, where dispatching strategies may not be applied in the DN with the grid constraints (such as the power flow). To ensure the performance of the MES and TES under the uncertainties of the RES output, etc., a novel robust optimization method is proposed. The robust optimization method using the uncertainty sets to represent the uncertainties has been a hot research topic, which is often used to analyse power system reserve dispatch considering the RES uncertainty [22,23]. Furthermore [24], proposes that the number of RES uncertain parameters that reach the bounds in the robust optimization model can reflect the risk-aversion level of the microgrid developer. In [25], the number of the price deviation period is used to express the conservative, moderate and aggressive attitudes of the dispatcher. However, unlike the existing dispatching system, in this paper, both the power network and traffic network will be significantly affected by the uncertain factors, resulting in the dispatching difficulty for the EMS. Additionally, for the nonconvex optimization model, “the worst scenario” integrating the risk-preference of the dispatcher may not always occur at the bounds, but within the uncertainty domain. With the observations mentioned above, this paper proposes a novel multi-energy-storage EMS for a DN with high renewable energy penetration. The MES is driven by the electricity to transit between stations. Meanwhile, the inverter AC integrating the building is developed into the TES with reactive and active regulation abilities. These abilities are

co-optimized by the EMS to improve the operational economy and voltage deviation of the DN. Three main obstacles exist for solving the EMS problem of the DN. First, compared with [4–21] focusing on either the MES or TES, this paper is devoted to formulating a unified model for the electricity-driven MES and inverter-AC-based TES and exploiting their complementary advantages. Second, unlike the existing robust methods in [22–25], where the uncertain factors at the bounds only worsen the power flow, in this paper, the forecasting errors will have a profound impact on the power network and traffic network, and “the worst scenario” may occur within the uncertainty domain. Third, in contrast to the existing studies in [4–12], where the optimization model is a classic mixed integer second-order cone programming (MI-SOCP) problem or a mixed integer linear programming (MILP) problem, the operational model of the DN in this paper incorporates complex nonlinear constraints, resulting in the difficulty in solution. To address the first two challenges, the MES is modelled considering the traffic delay factors and transit electricity consumption, while the inverter AC integrating the building is developed into the TES model with the charging and discharging functions. Considering the uncertain factors and risk-aversion of the dispatcher, a more accurate “worst scenario” is identified by discretizing uncertain domains into multiple virtual bounds and analysing their multidimensional impact. To address the third challenge, the optimization model is converted into a MISOCP problem and MILP problem, and a 3-stage iteration method is used to efficiently solve these problems. Hence, the highlighting features of the paper are summarized as follows: (1) A novel multi-energy-storage EMS for the DN is proposed for the first time. The electricity-driven MES is modelled considering both the power consumption and transit delay factor, while the inverterAC with the reactive power function is developed into a TES. These systems are formulated into a unified model with charging and discharging capabilities. (2) Considering the multidimensional impact of the uncertainties and risk-aversion of the dispatcher, a novel robust method oriented nonconvex optimization model is proposed. Benefiting from the proposed method, a more accurate “worst scenario” can be found by discretizing uncertain domains into multiple virtual bounds. (3) To efficiently solve the multi-energy-storage EMS, a linearized technique is presented to convert the optimization model into a MILP and MI-SOCP, and a three-stage iteration algorithm, i.e., the main-problem preprocessing stage, main-problem fine processing stage and sub-problem processing stage is used to accelerate the solution. (4) Using the realistic DN parameters, the effectiveness and benefits of the proposed multi-energy-storage EMS for the DN are illustrated. The rest of this paper is organized as follows: the structure of the multi-energy-storage EMS is introduced in Section 2. Section 3 introduces an optimization model conversion method and a two-stage iteration algorithm. Simulation studies are conducted in Section 4. Finally, the conclusions are drawn in Section 5. 2. Structure of multi-energy-storage energy management system Fig. 1 depicts the multi-energy-storage management framework for the DN. As observed from the figure, distributed RES generators, MES stations, smart buildings consisting of AC loads and conventional loads in the radial feeder exist. According to the social development trend, the MES will be driven by the electricity to move between stations, while the inverter AC will dominate the market and perform as the TES. In this paper, a multi-energy-storage EMS, together with the modified models of the MES and TES, is proposed to improve the operational economy and node voltage quality of the DN. The primary goal of the EMS is optimizing a predefined objective function, i.e.., minimizing the DN operational costs considering many 3

Electrical Power and Energy Systems 118 (2020) 105779

Q. Sui, et al.

Subsation

Photovoltaics Wind turbine

Smart house

N

F , t = F pdc ,t +

N

n=1 N

U: =

B1 Station 1

Station 2

F pdc ) u1n , t ,t

n (F min ,t

(u1n

n=1 T N

Station 3

t=1 n =1

,t

(u1n

+

+ u2n , t ) ,t

n (F max ,t

F pdc ) u2n , t ,t N

1

+ u2n , t )

(10)

where the first term expresses the possible occurrence value according to the prediction, the second term ensures the only deviation occurrence at time t with the binary variable u1 , t or u2 , t equaling 1, the third term expresses the total forecasting errors in scheduling cycle which reflects the risk aversion. Note that the actual deviation of parameters from their predicted values can be calculated according to the u1 , t or u2 , t . Additionally, when the “budget of uncertainty” is 0, the robust optimization model is completely equivalent to the deterministic optimization model.

Fig. 1. Multi-energy-storage management framework for the DN.

constraints and limitations associated with the DN, MES and TES. To mitigate the impact of forecasting errors of the RES output, load, electricity price and transit delay, a robust optimization method is used to evaluate the “worst scenario” according to the uncertainty interval in this paper. As a result, a multi-energy-storage energy management problem can be presented as follow:

max min f (x , y )

(1)

2.2. Modelling of the electricity-driven MES

s . t . g (x , y ) = 0

(2)

h (x , y)

(3)

In this paper, the location and charging/discharging power of the MES are embedded into the DN schedule. A traffic network with travel time is used to describe the transit process of the MES, formulated as the following:

u U

x, y

0

where f (x , y ) expresses the objective function, g (x , y ) denotes the equality constraint set, h (x , y ) represents the inequality constraint set, u presents binary variables related to the risk parameters, x consisting of x1, t and x2, t expresses the binary variables excluding the risk parameters, y denotes continuous variables. Accordingly, u , x and y can be descripted as follows:

]T

n T x1, t = [Iij, t , UT, j, t ]

c d c d T x2, t = [UMES, i, t , UMES, i, t , UTES , i, t , UTES, i, t ]

c PMES, i, t ,

y=

d PMES, i, t ,

ij Ai +

(6)

Iij, T -

max min U

+

T t=1

p2, t

(7)

i

d b + PMES, i, t ) - PMES, t

Iij,0 = 0

Iij, T

1

i

i

=0

t + Tij, t + Tij, t - 1

(14)

[Tij, t + Tij, t ] Uij, t

ij

t = 1, ..T (8)

Iij,1 Iij, t Iij, t

Uij, t =

(12) (13)

i

Iij, t

t =t

Tij, t

Tij, t + 1,

i

t = 1, i t = 2, ..T , i

1

(15) (16)

where (11) ensures the only MES location with the binary variable Iij, t equal to 1. (12)–(14) express the continuity of the MES trip. (15)–(16) represent the transmit duration constraints. The electricity-driven MES energy storage will change when travelling or exchanging power with power grids. The energy variation process can be described as follows:

0

p1, t Ppurchase, t + Clabor

c (PMES, i, t

t = 1, ..T - 1,

ij Ai

In this paper, the operational costs of the DN consist of the electricity purchase cost from the external grid, the daily MES driver cost and the life loss cost of the MES. The multi-energy-storage EMS integrated the robust method can be expressed as follows:

t=1

Iij, t = 0 ij Ai

Iij,1

(5)

2.1. Objective function

T

(11)

Iij, t + 1 -

Q MES, i, t , soc MES, t ,

c d PTES , i, t , PTES , i, t , QTES, i, t , socTES, i, t , ij, t , iji, t , Pi, t , Qi, t , Pij, i, t , Qij, i, t

t

ij Ai +

(4)

u = [u1n , t , u2n , t ] x = [x1, t , x2, t

Iij, t = 1 ij A

(9)

c PMES,i, t

c c Pmax UMES,i, t

d d Pmax UMES,i, t

d PMES,i, t

SMES Iii, t

Considering the errors between the prediction value of the RES output, load, electricity price and transit delay F pdc , t and actual value F , t , the existing robust method is generally used to evaluate the “worst scenario” according to the uncertainty interval F m, t in F , t F max and ,t reflect the risk aversion using the number of hours in which the actual values can reach the uncertainty set bounds. However, this paper presents a nonconvex optimization model (please refer to the MI-SOCP and MILP in Section 3), the “worst scenario” may occur within an uncertainty domain rather than at the uncertainty set bounds. To find a more accurate “worst scenario”, we finitely discretize the uncertain domain and expands virtual bounds into 2 N. The uncertain domain U can be expressed as follows:

Q MES,i, t

c d UMES ,i, t + UMES,i, t

b PMES, t =

(17)

0

(18)

SMES Iii, t

(19) (20)

Iii, t

Pv 1

Iii, t

(21)

i c d b PMES, i, t = PMES,i, t + PMES,i, t + PMES, t

(22)

2 PMES, i, t

(23)

+

2 QMES, i, t

EMES, t = EMES, t

1

S2MES c PMES, i, t

+ i

4

d b PMES, i, t + PMES, t

+ i

(24)

Electrical Power and Energy Systems 118 (2020) 105779

Q. Sui, et al.

socMES, t = EMES, t Emax MES

(25)

2 2 PTES , i, t + QTES, i, t

0

(26)

ETES , i, t = ETES, i, t

socMES, t

1

(27)

soc MES, T = soc MES,0

socTES, i, t =

where (17)–(23) denote the relationship between the charging power, discharging power, driving power, capacity and location, respectively. Eqs. (24)–(27) represent the energy storage constraints of the MES.

0

Ca,ik dTd,ik,t dt = (To, t

Td,ik,t ) Rik

ETES ,ik,t = Ca,ik (Tmax

Td,ik,t )

Emax TES,ik = Ca,ik (Tmax

0

ETES,ik,t

(28) (29)

WTES,ik,t

(30) (31)

Tmin )

max ETES ,ik

b PTES, ik ,t =

Td,ik,t ] = [To, t

k1 To, t Td, ik,t k2 R ik

Td,ik,t ] Rik

k1 l2

WTES,ik,t

k2 l1 k1

c PTES,ik, t

[Pmax TES,ik

b d PTES,ik, t UTES,ik,t

b c PTES,ik, t ] UTES,ik,t

d PTES,ik, t

c d UTES,ik, t + UTES,ik,t

0

(49)

1

Pi, t = P l, i, t + PTES, i, t + PMES, i, t Qi, t = Q l, i, t + QTES, i, t + Q MES, i, t

max

Pb

Pc

Pd

Uc

Ud

(50)

Pr, i, t

(51)

Qr, i, t

The on-load tap changer (OLTC) can regulate the voltage of node i and adjacent node j with the voltage ratio kijs . If without OLTC between node i and node j, the value of kijs equals 1. The electrical relationship between the two adjacent nodes can be expressed as (52). The second norm of current of the branch ij and the voltage of node j are represented in (53)–(54).

(33) (34)

vj, t (kijs ) 2 =

2(rij Pij, i, t + x ij Qij, i, t ) + (r 2ij + xij 2)

i, t

ij, t

(52)

(35)

ij, t

= Lij2, t

(53)

(36)

j, t

= V j2, t

(54)

The branch current can be calculated according to the apparent power and node voltage:

(37)

1

socTES, i, t

(48)

The RES output, load, charging and discharging power of the MES and TES determine the power injected into node i:

A TES model capable of ‘charging’, ‘discharging’ and regulating reactive power can be developed as:

0

(47)

max ETES, i, t ETES ,i

2.4. Modeling of the power flow

(32)

where (28) denotes the refrigeration power as determined by the electric power. Eq. (29) denotes the relationship between the refrigeration power, indoor temperature and outdoor temperature. Eqs. (30)–(32) denote the relationship between the actual energy storage and the maximum energy storage. With the small scheduling time unit, Eqs. (28)–(29) can be approximately processed as follows:

Ca,ik [Td,ik,t + 1

c d + PTES , i, t + PTES, i, t

-1

where H in (41) denotes E , , P , , , , , , Q and S. These new variables related to the aggregated TES have the same meaning as the original variables related to the AC. socTES, i, t in (48) denotes the SOC of the aggregated TES at node i. Considering the difference in the AC inverter capacity, (46) is treated conservatively according to (39). On the one hand, we can see that the TES model can be compatible with the MES model. The consistent operation models may reduce the management difficulty for the multi-energy storage EMS. On the other hand, when the users participate in the DN operation, they only need to c b d submit the electric power PTES,ik, t , Pmax and Pmax instead of the building parameters such as the heat capacity. Therefore, the proposed TES model can protect the user privacy to some extent.

The inverter AC is scheduled to participate in the DN operation in this paper. The energy storage property of the building integrating the inverter AC is represented as [17]:

k2l1) k1

(46)

E max

2.3. Modelling of the inverter AC-based TES

WTES,ik,t = k2 k1 × PTES,ik,t + (k1l2

2 STES ,i

= (Pij, i, t 2 + Qij, i, t 2 )

(55)

c d b PTES,ik,t = PTES,ik, t + PTES,ik,t + PTES,ik,t

(38)

2 2 PTES,ik, t + Q TES,ik,t

S2TES,ik

(39)

The relationship between the line inflow power and line outflow power can be expressed as follows:

c d + PTES,ik, t k2 k1 + PTES,ik,t k2 k1

(40)

Pij, j, t = Pij, i, t

E TES,ik,t = ETES,ik,t

1

ij, t

Eqs. (35)–(37) denote the constraints of charging power, discharging power and basic electric power maintaining indoor temperature at time t. Eq. (39) represents the relationship between the active power and reactive power. Eq. (40) expresses the energy storage variation limit. If the Ki AC loads exist at node i, the aggregated model of the TES can be expressed as:

Qij, j, t = Qij, i, t

Pi, t =

0

max [PTES ,i

b d PTES , i, t UTES, i, t

d PTES , i, t

c d UTES , i, t + UTES, i, t

PTES, i, t =

c PTES , i, t

b c PTES , i, t ] UTES, i, t

+

0

1 d PTES , i, t

+

b PTES , i, t

(57)

ij, t

Pji, j, t j

Pij, i, t j

Qi, t =

(41)

k=1 c PTES , i, t

xij

Qji, j, t j

HTES, ik, t

(56)

ij, t

The relationship between the power injected into node, line inflow power and line outflow power can be expressed as follows:

Ki

HTES, i, t =

rij

i, t

Qij, i, t j

(58) (59)

The power loss can be expressed as follows: (60)

(42)

P loss, ij, t = rij

(43)

The line current and nodal voltage have operational limits [26,27], as follows:

(44)

ij, t

(45)

j

5

¯ij j, t

ij, t

(61)

¯j

(62)

Electrical Power and Energy Systems 118 (2020) 105779

Q. Sui, et al.

3. Solution algorithm

(47)–(51), (56)–(62) and (63)–(64), the fifth term expresses (23), (46) and (68). (x , u) is the set of continuous variables y with given x and u . A1, A2, C, B2, B3 and B4 represent coefficient matrixes, b1, b2, b3, and b4 denote the constant vectors. To this end, the robust optimization model can be treated as the classic MI-SOCP and MILP. Note that the transit time of the MES affected by the traffic congestion in (15) may not be an integral multiple of an hour, while the linearization from (29) to (33) may result in large errors. Therefore, a smaller dispatching time unit needs to be adopted to improve the resulting accuracy.

In this section, the multi-energy-storage EMS with nonlinear terms such as (15), (42), (43), (52) and (55) is converted to the MI-SOCP and MILP, which are solved by a 3-stage iterative method. 3.1. Optimization model conversion To deal with the step-by-step voltage regulation of the OLTC in (52), a simpler model based on the pinch principle rather than the existing large M method is proposed in this paper: 5

v js, t =

i, t

- 2(rij Pij, i, t + x ij Qij, i, t ) + (rij2 + x ij 2)

ij, t

3.2. Iterative solution algorithm

(63)

s=1

Although the classic MI-SOCP, etc., can be directly solved using commercial solvers, the numerous optimized Boolean variables from the small dispatching time unit processing will reduce the solution efficiency. To solve the problem, a novel iterative solution framework including the main-problem preprocessing stage, main-problem fine processing stage and sub-problem solution stage is proposed in this paper. Step 1: Initialization The parameters including loads, RES output, DN data, MES parameter, road congestion delay, AC parameter, outdoor temperature, building parameter, etc, are input. In addition, initial iteration number iter = 1, maximum iteration number itermax, and error ( > 0) are set. Additionally, the initial value 0 is given. Step 2: Main-problem preprocessing stage We observe that the binary variables related to the OLTC and MES only switch their values several times. These variables can be regarded as slow-action variables. Therefore, in the preprocessing stage, the slow-action variables are optimized on a 1-hour time scale to accelerate the solution. The objective function and constraints are formulated as follows:

5

v js, t (kijs ) 2

vj, t =

(64)

s=1

v js, t

0

UT,s j, t

(65)

M, s = 1, 2, ...5

5

UT,s j, t = 1

(66)

s =1

where the first term expresses that the actual voltage of node j is the sum of the voltages at different levels, and the latter terms ensure the only voltage level with the binary variable UT,s j, t equaling 1, where the ratio of the original voltage to the secondary voltage of OLTC kijs belongs to {1, 1 - 2.5\% , 1 + 2.5\% , 1 - 5\% , 1 + 5\% } . Ref. [12] notes that (55) can be approximately converted into the quadratic constraint (67) in the general distributed system, and (67) can be equivalently treated as the second-order cone limit in (68). (67)

Pij, i, t 2 + Qij, i, t 2

ij, t i, t

2Pij, i, t 2Qij, i, t ij, t

ij, t i, t

+

i, t

(68)

2

min T xiter + v iter s. t . A1x iter b1

In addition, (42)–(43) can be converted into the mixed integer linear constraints [28]: c PTES , i, t

0

b PTES , i, t

max c PTES , i UTES, i, t

c MUTES , i, t bc PTES , i, t

0

(1

bd PTES , i, t b PTES , i, t

d MUTES , i, t bd PTES , i, t

0

(69)

0

bd PTES , i, t

b PTES , i, t

d UTES , i, t )M

(1

(70)

where represents the value of and expresses d d U . the value of PTES , i, t TES, i, t Through the above treatment, the objective function of the robust optimization model can be formulated as max min C T y . Considering c c PTES , i, t UTES , i, t ,

bc PTES , i, t

u U

bd PTES , i, t

min s. t .

x, y

min x

s. t .

+ max u U

y

min

(x , u)

A1 x b1 A2 x + B2 y b2, B3 y b3, B4 y 2 b4 y

(72)

b2)

T x iter + v iter A1x iter b1

1

iter T x iter = [x1,iter t , x 2, t ]

the solution difficulty due to the strong coupling constraints such as (15), a decomposition optimization model is proposed in this paper as follows: Tx

iter 1 (A x iter 2

where iter 1 and viter 1 are the multiplier vectors and cut constraints from the sub-problem, respectively. By solving the main problem, the slow-action variable results x1,iter t on the long-time scale can be obtained. Step 3: Main-problem fine processing stage On a short-time scale (20 min), the values of slow-action variables may switch at non-integer times near the integer hour. Therefore, in the fine processing stage, the slow-action variables are finely modified to improve the result accuracy.

b PTES ,i

c UTES , i, t )M

d PTES , i, t

iter T x iter = [x1,iter t x 2, t ]

v iter 1

bc PTES , i, t

bc PTES , i, t

1

, iter x1,iter t = x 1, t ,

v iter

CT y

1

t

iter 1 (A x iter 2

b2)

(73)

where Ψ is the time set related to the value switching. For example, if x1,4,iter equals 1 while x1,5,iter equals 0, only the values of x1iter at 4:20, 4:40, 5:00, 5:20 and 5:40 need to be optimized while the others (for example, 4:00, 3:40, etc.) can be directly determined by x1 , iter . (x iter , v iter 1) The results of can be obtained, and LBiter = T x iter + v iter 1 can be updated. Step 4: Sub-problem processing stage Based on the KKT condition and strong duality theory, the subproblem can be converted as follows:

(71)

where the first term expresses (9), the second term represents (11)–(16), (20), (44) and (66), the third term represents (17)–(19), (21), (65),(69),(70), the forth term denotes (22), (24)–(27), (41), (45), 6

Electrical Power and Energy Systems 118 (2020) 105779

Q. Sui, et al.

max

uiter U

s. t .

iter B

2

iter (A x iter 2 iter B

3

+ B4

T iter

iter iter ,

iter ,

iter b

b2)

2 iter ,

b4

iter

simulation results show that the road congestion factor has a significant impact on the transit time and energy consumption. Therefore, for the electricity-driven MES, the transmit energy consumption must be fully reflected and well managed by the EMS. To explain and illustrate the energy storage characteristics of the TES clearly, the authors simulate the TES SOC variation excluding the constraint (49) in the one-hour dispatching cycle by specifying the initial SOC value of 1 at 0:00, 1:00, 2:00, etc., respectively. The simulation result is shown in Fig. 5(b). As depicted in the figure, the energy storage loss of the TES is significantly affected by the ambient temperature in case of maintaining a comfortable indoor temperature, and the greater the temperature difference, the greater the thermal storage pressure. During 11:00–19:00, the SOC exceeds the lower limit, implying that the indoor temperature deviates from the comfort zone. Therefore, a reasonable scheduling strategy is important.

3

+C=0

iter iter

0

(74)

where uiter , iter , iter , iter and iter are variables. By solving the sub-problem, the result iter can be acquired, and iter b can be updated. UBiter = Tx iter + iterT (A2 x iter b2) 3 Step 5: Check the stopping criteria There are two stopping criteria in the solution procedure: Rule 1: iter > itermax Rule 2: UBiter LBiter If either criteria is satisfied, the iteration process will be terminated. Otherwise, update iter = iter + 1 and go to Step 2. Although the steps above have well described the solution process, we further draw a flow chart to illustrate it more clearly, as shown in Fig. 2.

4.2.2. Regulation property of the MES and TES Fig. 6 shows the power regulation property of the MES and TES according to the typical outdoor temperature in Fig. 5(b). As illustrated in Fig. 6(a), the power regulation region of the MES makes up a standard cylinder with the 24-hour high and 2-MW radius. However, affected by the outdoor temperature, the TES can carry out the power regulation within an irregular cylinder in Fig. 6(b). For instance, to maintain the indoor temperature of 24 °C, the maximum charging/ discharging power and adjustable reactive power of the TES are 0.73/ 2.20 MW and 2.79/1.86 MW, respectively at the outdoor temperature of 38 °C, while those parameters are 2.50/0.43 MW and 1.43/2.85 MW, respectively at the outdoor temperature of 26.7 °C. The results show that the MES can serve the power grid with a fixed power regulation, while the power regulation service of the TES is affected by the ambient temperature.

4. Numerical simulations 4.1. Test system A 41-node DN in Ontario shown in Fig. 3 is introduced to verify the effectiveness of the proposed strategy [28]; the rated voltage (i.e., 16 kV) and branch parameters including the length and resistance of this DN can be found in Table A1 in Appendix A [29]. The allowed voltage deviations are −5 to 5% and the maximum branch current is 0.71 kA. The electricity-driven MES is integrated with the lithium battery, whose rated energy storage capacity, charging/discharging power/efficiency, driving power and travel speed are set to 2.2 MWh, 2 MW/0.9, 20 kW and 40 km/h, respectively. The investment cost and life cycle times of the lithium battery are 2700 yuan/kWh and 4000, respectively [30]. The frequency converter capacity related to the MES is 2 MW. According to the current data of the Chinese labour market, the driver cost is set to 300 yuan/day. Additionally, Station 1, 2 and 3 for the MES are located at nodes 9, 28 and 40, respectively. The distances between Station 1–2, 1–3 and 2–3 are 17 km, 25 km and 16 km, respectively. Additionally, the TES parameters are shown in Tables A2 and A3, where the building thermal resistance and thermal capacity are distributed randomly in their range. The typical outdoor temperature is illustrated in Table A4 [3]. The hourly average value of the net common load (common load minus renewable energy generation) is shown in Table A5. Fig. 4 shows the forecasting curves and uncertainty sets of the RES output, load, electricity price and road congestion delay. Their forecasting errors are [−15%, 15%], [−10%, 10%], [−10%, 10%] and [−30%, 30%] [31–33], respectively. The scheduling time unit is set to 20 min. Based on a 4.0 GHz Windows-based PC with 8 GB of RAM, all numerical simulations are coded in MATLAB and solved using Gurobi called by Yalmip.

4.3. Co-optimization effect analysis of the MES and TES Without losing generality, assume r = 20*15%, l = 25*10%, = 15*10%, d = 20*30% and N = 2 in (10). The MES and TES are cooptimized to improve the operational economy and voltage quality of p

Initialization Set iter=1, itermax, ,

0

Main-problem preprocessing stage Optimize variables x1,titer on a 1-hour scale Obtain slow-action variables x1,titer *,iter iter x1,t = x1,t

Main-problem fine processing stage Optimize variables on a 20-min scale based on x1,t*,iter iter Update LB = Txiter+v iter 1 Feedback cut constraints v iter

x iter

1

Sub-problem processing stage Optimize variables iter and iter on a 20-min scale Update UBiter= Txiter+ iterT(A2 xiter b2) iterb3

4.2. Analysis of basic properties of the MES and TES To exploit the complementary advantages of the MES and TES, the energy storage and power regulation properties are analysed.

No

4.2.1. Energy storage property of the MES/TES If the MES with the initial SOC of 0.50 departs from Station 1 to Station 2 at 0:00, 1:00, 2:00, etc., respectively, the transit time and storage energy variation in the one-hour dispatching cycle are shown in Fig. 5(a). As observed from the figure, MES transit duration varies from 39 to 60 min, and the SOC decreases from 0.4934 to 0.49. The

iter > itermax iter

UB

iter

LB <

Yes Output current UBiter as the final results Fig. 2. Flow chart of algorithm solution. 7

Electrical Power and Energy Systems 118 (2020) 105779

6 2

1

WT3 41

3

40

4 PV1

39

ST3

5 8 37

38

PV2

Load

10 13 17 9 11 12 14 15 16 18 ST1 34 WT1 19 32 23 35 WT2 33 24 25 36 26 28 30 29 ST2 31 MES station PV 7

20

21 22

27

(a) Power property of the MES

Wind turbine

Active power (MW)

interconnection node

Active power (MW)

Q. Sui, et al.

Fig. 3. Structure of the DN integrated the MES and TES.

4.0

(b) Power property of the TES Fig. 6. Power property of the MES and TES.

1.5 1.0

1. 0

0

(3) To make the revenue larger than the life loss cost, the MES is only charged at the electricity-price valley period and discharged at the electricity-price peak period. By contrast, the TES with little life loss frequently exchanges power with the DN in the neighboring period, but fails to store energy in the long period due to its serious heat loss. The results show that the MES and TES are co-optimized to support the DN exploiting their complementary advantages.

0. 5 0

2

4

6

8

10

12 14 Time (h)

16

18

20

0 24

22

RES output uncertainty domain Load uncertainty domain Electricity price uncertainty domain Delay time uncertainty domain Forecasting mean value of load, RES output, electricity price and delay time

SOC

0.500

50

0.496

45

0.492

40 0

2

4

6

8

14 10 12 Time (h)

16

18

20

22

(a)Energy property of the MES SOC

24

0.488

1.000

Outdoor temperature

35

0.500

30

0.000

25

0

2

4

6

8

SOC exceeds lower limit 10 12 14 16 18 Time (h)

20

22

4.4. Comparison with the existing optimization method To verify the superiority of the proposed MES and TES co2000 Price (yuan)

35

40

Temperature ( )

0.504

Transit duration

SOC

60 55

The power distributions of the MES and TES shown in Fig. 8 further verify their complementary properties. As observed from the figure, transiting among stations, the MES has the discrete power regulation property, providing the strong reactive and active power support for nodes 9, 28 and 40 when connected with the power grid. On the other hand, although the power support from the TES based on the inverterAC is limited, the TES continuously regulates the power in numerous nodes. The simulation results indicate the MES and TES have complementary properties in the space.

SOC

Transit duration (min)

Fig. 4. Basic data of the RES output, load, electricity price and traffic delay.

-0.500 24

(b)Energy property of the TES

1000

2

0

1

2 Power (MW)

Fig. 5. Energy property of the MES and TES.

Power (MW)

the DN. The results are analysed as follows. Fig. 7 shows the operational power and storage energy of the MES and TES. From the figure, we have the following observations:

Sotrage energy(MWh)

(1) In the dispatching cycle, the MES spends 4.33 h and 86.67-kWh electricity travelling along the route 3–1–2–1–2–3. The transmit energy consumption is directly reflected in the SOC. On the other hand, the TES generates the reactive power to participate in the DN operation all the time, implying the AC-DC-AC structure of the inverter air-conditioning is well developed and utilized. (2) The MES exchanges the 3.15-MWh active power and 7.14-MWh reactive power with the DN. Meanwhile, the TES provides the 3.71MWh active power and 51.92-MWh reactive power for the DN. We can see that the MES and TES with the same charging, discharging and reactive regulating functions can be uniformly dispatched to carry out peak cutting and valley filling.

3

MES location Electricity price

Active power for MES Reactive power for MES

1 0 -1 -2 1 0

Active power for TES Reactive powerfor TES

-1 -2 -3 2.20

Energy for MES Energy for TES

1.76 1.32 0.88 0.44 0

0

2

4

6

8

10 12 14 Time (h)

16

18

20

22

24

Fig. 7. Operational power and storage energy of the MES and TES. 8

Station

Price (yuan/kW)

0 2.0

delay time (h)

Power (MW)

8.0

Electrical Power and Energy Systems 118 (2020) 105779

Q. Sui, et al.

over the existing method. Fig. 10 and Table 1 show the node voltage of the DN for the four schemes. The maximum voltage deviation is 14% for Scheme IV, implying the DN voltage quality is really poor. In the actual case, the renewable energy generators may even be off grid. Although the maximum voltage deviation is reduced to 8% for Scheme II and 9% for Scheme III, the voltage quality of the DNs fails to meet the operational requirement, indicating that either the MES or TES can improve the DN operation to some extent, but the effect is limited. Unlike Scheme II and Scheme III, the DN voltage deviation for Scheme I is maintained within 5%. The good voltage quality is achieved by co-optimizing the MES and TES. Additionally, Table 1 shows the power loss is 6.11 MWh for Scheme I, 6.67 MWh for Scheme II, 7.54 MWh for Scheme III and 7.85 MWh for Scheme IV, respectively. In other words, the power loss of Scheme I can be reduced by 8%, 19% and 22% compared with that of Scheme II, III and IV, implying the proposed strategy is effective in optimizing power flow distribution and reducing power loss.

Power (kW)

Active power Reactive power

Node (a) MES Active power Reactive power

4.5. Comparison with the existing EMS To analyse the superiority of the proposed EMS that integrates the robust method, three existing EMSs have been evaluated for comparison as follows: Case 1: A deterministic EMS is used. The objective function is shown = 0, N = 0 . in (9) and the parameters in (10) satisfy F , t = F pdc ,t , Case 2: Considering forecasting errors, a classic EMS that integrates robust method is used to evaluate “worst scenario” which occurs at the boundary of the uncertain domain. The objective function is shown in (9) and the parameters in (10) are set as = T, N = 1. Case 3: Considering forecasting errors and the risk preference of the dispatcher, a classic EMS that integrates robust method is used to analyse “worst scenario” which occurs at the boundary of the uncertain domain. The objective function is shown in (9) and the parameters in 15%, l = 25 10%, p = (10) are set as r = 20 15 10%, d = 20 30% and N = 1. Table 2 shows the daily operational costs of the DN as monitored by several EMSs. As observed from the table, the daily operational costs of the DN for Case 1, Case 2, Case 3 and the proposed case in this paper are 70,615 yuan, 76,431 yuan, 73,282 yuan, and 73,355 yuan, respectively. The simulation results show that the operational cost of the DN for Case 1 is lowest among all the cases, implying uncertainty factors will damage the operational economy of the DN. Additionally, compared with that for case 3 and case 4, the operational cost for Case 2 is higher, indicating that the risk preference has a significant impact on the operational economy of the DN, and the more conservative, the worse the economy. Lastly, with the same risk-aversion level, the operational cost using the proposed strategy is higher than that for Case 3. The result show that the proposed EMS in this paper can obtain a more accurate “worst case scenario” to ensure the DN safe operation. The operational result of the DN guided by the proposed EMS is shown in Table 3 to analyse the impact of the risk-aversion level. As depicted in the table, with the uncertain parameters increasing, the purchased electricity from the external grid increases, while the daily

0

-0.1

-0.2 0

4

8

12

16 20 24 0 9

18 27

36 45

(b) TES Fig. 8. Power distribution of the MES and TES.

optimization strategy (defined as Scheme I), three existing schemes are introduced for comparison as follows: (1) Scheme II: The MES is scheduled to optimize the DN operation while the AC acts as an unadjustable load. (2) Scheme III: The DN is optimized using the TES whose reactive power regulation potential is not exploited, and the MES is excluded in the DN operation. (3) Scheme IV: Neither the MES nor the TES is utilized to participate in the DN operation. The DN can only exchange power with the external grid to maintain the power balance. Fig. 9 shows the exchange power between the DN and external grid for Scheme I, Scheme II and Scheme III relative to Scheme IV. As observed from the figure, for Scheme I, both the MES and TES are utilized to carry out peak cutting and valley filling, resulting in an elastic exchange power to the electricity price. For Scheme II, to achieve the energy arbitrage, the MES may charge at the low-price-electricity period and discharge at the high-price-electricity period to compensate for the battery life loss cost. Therefore, the exchange power is relatively concentrated in several periods. For Scheme III, limited by the serious heat leak, the TES cannot store much electricity for a long time. The exchange power floats in the adjacent periods to reduce the electricity consumption cost. In the dispatching cycle, the shifted electricity is 17.79 MWh for Scheme I, 12.27 MWh for Scheme II, 2.67 MWh for Scheme III and 0 MWh for Scheme IV, respectively. Meanwhile, the operational costs are 73,355 yuan for Scheme I, 74,096 yuan for Scheme II, 74,804 yuan for Scheme III and 75,311 yuan for Scheme IV, respectively, that is, compared with Scheme II-IV, the operational economy for Scheme I is improved by 1.0%, 1.9%, 2.6%, respectively, implying the proposed strategy has a significant economic advantage

1

2000 Scheme I Scheme II Scheme III 1500 Price

0

1000

3

Power (MW)

2

-1

500

-2 -3 0

2

4

6

8

10

12 14 Time (h)

16

18

20

22

24

Price (yuan)

Power (kW)

0.1

0

Fig. 9. Exchange power comparison between the DN and external power. 9

Electrical Power and Energy Systems 118 (2020) 105779

Q. Sui, et al.

1.05

Voltage (p.u.)

Voltage (p.u.)

1.05

1.00

1.00 0.95 0.90 0 15

0.95 0 15 30

20

45 24

16 12

4

8

0

45

1.05

1.05

1.00

1.00

0.95

15

30

20

45 24

16

12

24

20

16

12

8

4

0

(c) DN’s node voltage in Scheme III

Voltage (p.u.)

Voltage (p.u.)

(a) DN’s node voltage in Scheme I

0.90 0

30

0.95 0.90 0.85 0

0 8 4

15

(b) DN’s node voltage in Scheme II

30

45 24 20

8

16 12

4

0

(d) DN’s node voltage in Scheme IV

Fig. 10. DN’s node voltage comparison. Table 1 Comparison for voltage deviation and power loss.

Table 3 Risk preference comparison.

Schemes

Maximum voltage deviation

Power loss (MWh)

Scheme Scheme Scheme Scheme

5% 8% 9% 14%

6.11 6.67 7.54 7.85

I II III IV

Table 2 Comparison for daily operational costs. Cases

Daily operational cost (yuan)

Case 1 Case 2 Case 3 The proposed case

70,615 76,431 73,282 73,355

Гr

Гl

Гp

Гd

Purchased electricity (MWh)

Daily operational costs (yuan)

0 15 * 15% 50 * 15% 72 * 15%

0 15 * 10% 50 * 10% 72 * 10%

0 15 * 10% 50 * 10% 72 * 10%

0 15 * 30% 50 * 30% 72 * 30%

133.14 135.90 142.17 145.51

70,615 73,109 75,495 76,595

Table 4 Optimum comparison for iterative problem and original problem. Operational cost (yuan) Iterative optimization problem Original optimization problem

operational cost rises. This is because that to ensure the energy supply for the load, the DN must obtain more electricity support from the external grid to mitigate the impact of the increasing uncertainty, thus the operational economy gradually decreases. Additionally, as the uncertainty decreases, the marginal revenue increases. The result shows the law “the greater the risk, the greater the benefit” can be applied to the DN operation.

73,355 73,008

Table 5 Comparison for existing methods for the TES/MES.

Existing MES method Proposed MES method Existing TES method Proposed TES method

10

Operational cost (yuan)

Computing time (second)

74,108 74,096 74,801 74,804

186 165 127 134

Electrical Power and Energy Systems 118 (2020) 105779

Q. Sui, et al.

4.6. Discussion on the methodology

appropriate place. However, compared with the stationary battery, the MES can serve as an UPS for some intermittent loads in the normal scenario, and provide the power support for important loads in the disaster scenario. Therefore, the MES is valuable.

To verify the feasibility of the proposed three-stage iterative algorithm, the solution to the original non-convex problem maxmin CT y is u U x, y

studied. Because the non-convex optimization model and the strong coupling constraint (15), the strong duality theory cannot be used. We design a two-layer method to solve the original non-convex problem, where the outer layer uses the genetic algorithm to optimize the risk parameter variables u to maximize the function value, and the inner layer uses Gurobi to optimize the variables x and y with the objective min CT y . The optimization results are shown Table 4. As observed from

5. Conclusions In this paper, we have proposed a novel multi-energy storage EMS for DNs to co-optimize the electricity-driven MES and inverter AC-based TES. The MES considering the delay factors and TES with the reactive power regulation have been developed into a unified analytic model capable of charging and discharging. Then a novel robust optimization method is proposed to evaluate the uncertainty by discretizing uncertain domains into multiple virtual bounds, followed by a three-stage iterative solution algorithm. In the end, the performance of the multienergy storage EMS is simulated on a 41-node DN. The results are as follows: (1) The compatible model for the MES and inverter-AC-based TES can facilitate the EMS management; (2) The MES and TES can reduce the power loss by no less than 8%, decrease the operational cost by no less than 1%, and improve the voltage level with a voltage deviation of 5%; (3) The proposed EMS integrated robust method can find a more accurate “worst scenario” to ensure the DN safe operation. Although the good performance of the multi-energy storage EMS has been verified by the simulation, the work ignores the energy storage owners’ service will and conflict interest with the DN operator. Therefore, our future work will focus on the demand response to guide user behaviors. Novel dispatching strategies to guide the MES and TES to participate in power grid frequency regulation will be also our future work.

x,y

the table, the optimality gap between the original non-convex problem and the proposed iterative optimization problem is 347 yuan (0.5%). Additionally, the proposed iterative optimization problem is solved within 241 s. To compared with the proposed methods (please refer to Scheme II and III), the direct optimization method for the MES (the main problem is directly solved) in [8] and the common optimization method for the based-AC TES in [3] are studied. The simulation results are shown in Table 5. As observed from the table, the operational costs of the DN using existing and proposed methods are basically the same. However, for the MES, compared with the direct solution, the computing time using proposed iterative method (the main problem consists of the preprocessing and fine processing stage) can be reduced by 21 s (11%), implying the proposed method provides a fairly accurate result in a relatively short time. Note that a shorter dispatching step such as 10 min will make the solution method advantage more prominent. On the other hand, for the TES, the computing time using the proposed method can be increased by 7 s (6%) compared with that using the existing method. However, the consistent operation models between the TES and MES model may reduce the management difficulty for the multi-energy storage EMS. Meanwhile, the proposed TES model can protect the user privacy to some extent because the users only need to submit the electric power instead of the building parameters when participating in the DN operation. In fact, the MES may not always be better that the stationary energy storage regarding the operational economy. On the one hand, based on the given parameters, a 2.8-MWh/2.5-MW stationary battery located in Station 2 can achieve the voltage deviation within 5%. The operational cost of the DN is 73,786 yuan, which is higher than the operational cost of 73,355 yuan using the MES for Scheme I. On the other hand, when the load is reduced by 10%, the voltage quality of the DN can be satisfied using both the proposed MES and a 2.4-MWh/2.2-MW stationary battery located in Station 2. The 65854-yuan operational cost of the DN for the MES is higher than the 65729-yuan cost for the TES. The above simulation results show that the MES may not be cheaper than an additional battery storage system of an appropriate size and at an

CRediT authorship contribution statement Quan Sui: Conceptualization, Methodology, Formal analysis, Writing - original draft. Fanrong Wei: Validation, Writing - review & editing. Xiangning Lin: Supervision. Chuantao Wu: Resources, Data curation. Zhixun Wang: Software, Visualization. Zhengtian Li: Project administration, Funding acquisition. Declaration of Competing Interest The authors declared that they have no conflicts of interest to this work. Acknowledgements This work was supported by the National Natural Science Foundation of China (51537003).

Appendix A See Tables A1–A5. Table A1 Branch parameters of DN. From

To

Length (km)

R (ohm/km)

X (ohm/km)

l (kA2)

1 2 3 4 5 5 7 9 9

2 3 4 5 6 7 9 10 11

5.7 1.01 0.4 0.38 0.13 0.17 0.26 0.14 0.38

0.169111 0.169111 0.169111 0.169111 0.169111 0.169111 0.169111 0.169111 0.169111

0.418206 0.418206 0.418206 0.418206 0.418206 0.418206 0.418206 0.418206 0.418206

0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5

(continued on next page) 11

Electrical Power and Energy Systems 118 (2020) 105779

Q. Sui, et al.

Table A1 (continued) From

To

Length (km)

R (ohm/km)

X (ohm/km)

l (kA2)

11 12 12 14 16 17 18 19 21 19 23 24 24 26 26 28 29 28 23 32 33 33 35 35 37 38 39

12 13 14 15 17 18 19 20 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

0.56 0.3 3.33 1.03 1.08 1.64 0.47 0.47 0.96 0.19 1.94 2.45 1.63 1.2 2.12 0.73 0.75 2.54 0.36 0.26 3.58 0.77 2.08 4.51 3.24 0.3 0.5

0.169111 0.169111 0.169111 0.169111 0.169111 0.169111 0.169111 0.348124 1.391924 0.348124 0.348124 0.348124 0.348124 0.552276 0.348124 0.552276 0.552276 0.348124 0.276519 0.276519 0.552276 0.276519 0.348124 0.276519 0.169111 0.169111 0.169111

0.418206 0.418206 0.418206 0.418206 0.418206 0.418206 0.418206 0.468482 0.478811 0.468482 0.468482 0.468482 0.468482 0.485241 0.468482 0.485241 0.485241 0.468482 0.45858 0.45858 0.485241 0.45858 0.468482 0.45858 0.418206 0.418206 0.418206

0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5

Table A2 TES parameters. R

2.5–3.5 °C/kW

Ca

2.0–3.0 kWh/°C

k1 l1 Tmax STES

0.03 kW/Hz −0.4 kW 26 °C 3.6 kW

k2 l2 Tmin max PTES

0.78 kW/Hz −0.8 kW 22 °C 3.6 kW

Table A3 TES number at different DN nodes. Node Number Node Number Node Number Node Number

1 0 12 0 23 0 34 50

2 0 13 50 24 0 35 0

3 50 14 50 25 50 36 50

4 0 15 0 26 0 37 50

5 0 16 0 27 50 38 0

6 100 17 0 28 0 39 0

7 0 18 0 29 0 40 0

8 50 19 0 30 50 41 50

9 0 20 0 31 50

10 50 21 0 32 0

11 0 22 50 33 0

Table A4 Outdoor temperature. Time/h T/°C Time/h T/°C

1 28.0 13 37.3

2 27.6 14 37.6

3 27.1 15 38.0

4 26.7 16 37.7

5 26.8 17 37.2

6 27.7 18 36.6

12

7 28.9 19 35.7

8 30.5 20 34.7

9 32.0 21 33.5

10 33.8 22 32.0

11 35.0 23 30.4

12 36.2 24 28.5

Electrical Power and Energy Systems 118 (2020) 105779

Q. Sui, et al.

Table A5 Net common load power (MW) at different DN nodes. Node

1h

2h

3h

4h

5h

6h

7h

8h

9h

10 h

11 h

12 h

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41

0 0 0.22 0 0 0.44 0 0.22 0 0.22 0 0 0.22 0.22 0 0 0 0 0 0 0 0.06 0 0 0.06 0 0.06 0 0 0.06 0.06 0 0 0.69 0 0.69 0.69 0 0 0 0.69

0 0 0.19 0 0 0.38 0 0.19 0 0.19 0 0 0.19 0.19 0 0 0 0 0 0 0 0.06 0 0 0.06 0 0.06 0 0 0.06 0.06 0 0 0.71 0 0.71 0.71 0 0 0 0.71

0 0 0.16 0 0 0.32 0 0.16 0 0.16 0 0 0.16 0.16 0 0 0 0 0 0 0 0.07 0 0 0.07 0 0.07 0 0 0.07 0.07 0 0 0.68 0 0.68 0.68 0 0 0 0.68

0 0 0.13 0 0 0.26 0 0.13 0 0.13 0 0 0.13 0.13 0 0 0 0 0 0 0 0.07 0 0 0.07 0 0.07 0 0 0.07 0.07 0 0 0.62 0 0.62 0.62 0 0 0 0.62

0 0 0.10 0 0 0.20 0 0.10 0 0.10 0 0 0.10 0.10 0 0 0 0 0 0 0 0.08 0 0 0.08 0 0.08 0 0 0.08 0.08 0 0 0.59 0 0.59 0.59 0 0 0 0.59

0 0 0.11 0 0 0.22 0 0.11 0 0.11 0 0 0.11 0.11 0 0 0 0 0 0 0 0.09 0 0 0.09 0 0.09 0 0 0.09 0.09 0 0 0.43 0 0.43 0.43 0 0 0 0.43

0 0 0.12 0 0 0.24 0 0.12 0 0.12 0 0 0.12 0.12 0 0 0 0 0 0 0 0.15 0 0 0.15 0 0.15 0 0 0.15 0.15 0 0 0.41 0 0.41 0.41 0 0 0 0.41

0 0 0.11 0 0 0.22 0 0.11 0 0.11 0 0 0.11 0.11 0 0 0 0 0 0 0 0.27 0 0 0.27 0 0.27 0 0 0.27 0.27 0 0 0.52 0 0.52 0.52 0 0 0 0.52

0 0 0.10 0 0 0.20 0 0.10 0 0.10 0 0 0.10 0.10 0 0 0 0 0 0 0 0.47 0 0 0.47 0 0.47 0 0 0.47 0.47 0 0 0.51 0 0.51 0.51 0 0 0 0.51

0 0 0.12 0 0 0.24 0 0.12 0 0.12 0 0 0.12 0.12 0 0 0 0 0 0 0 0.56 0 0 0.56 0 0.56 0 0 0.56 0.56 0 0 0.36 0 0.36 0.36 0 0 0 0.36

0 0 0.15 0 0 0.30 0 0.15 0 0.15 0 0 0.15 0.15 0 0 0 0 0 0 0 0.56 0 0 0.56 0 0.56 0 0 0.56 0.56 0 0 0.32 0 0.32 0.32 0 0 0 0.32

0 0 0.19 0 0 0.38 0 0.19 0 0.19 0 0 0.19 0.19 0 0 0 0 0 0 0 0.56 0 0 0.56 0 0.56 0 0 0.56 0.56 0 0 0.27 0 0.27 0.27 0 0 0 0.27

Node

13h

14h

15h

16h

17h

18h

19h

20h

21h

22h

23h

24h

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

0 0 0.21 0 0 0.42 0 0.21 0 0.21 0 0 0.21 0.21 0 0 0 0 0 0 0 0.56 0 0 0.56 0 0.06 0 0 0.06

0 0 0.17 0 0 0.34 0 0.17 0 0.17 0 0 0.17 0.17 0 0 0 0 0 0 0 0.54 0 0 0.54 0 0.06 0 0 0.06

0 0 0.14 0 0 0.28 0 0.14 0 0.14 0 0 0.14 0.14 0 0 0 0 0 0 0 0.56 0 0 0.56 0 0.07 0 0 0.07

0 0 0.12 0 0 0.24 0 0.12 0 0.12 0 0 0.12 0.12 0 0 0 0 0 0 0 0.56 0 0 0.56 0 0.07 0 0 0.07

0 0 0.11 0 0 0.22 0 0.11 0 0.11 0 0 0.11 0.11 0 0 0 0 0 0 0 0.54 0 0 0.54 0 0.08 0 0 0.08

0 0 0.15 0 0 0.30 0 0.15 0 0.15 0 0 0.15 0.15 0 0 0 0 0 0 0 0.51 0 0 0.51 0 0.09 0 0 0.09

0 0 0.25 0 0 0.50 0 0.25 0 0.25 0 0 0.25 0.25 0 0 0 0 0 0 0 0.47 0 0 0.47 0 0.15 0 0 0.15

0 0 0.32 0 0 0.64 0 0.32 0 0.32 0 0 0.32 0.32 0 0 0 0 0 0 0 0.40 0 0 0.40 0 0.27 0 0 0.27

0 0 0.37 0 0 0.74 0 0.37 0 0.37 0 0 0.37 0.37 0 0 0 0 0 0 0 0.33 0 0 0.33 0 0.47 0 0 0.47

0 0 0.40 0 0 0.80 0 0.40 0 0.40 0 0 0.40 0.40 0 0 0 0 0 0 0 0.27 0 0 0.27 0 0.56 0 0 0.56

0 0 0.35 0 0 0.70 0 0.35 0 0.35 0 0 0.35 0.35 0 0 0 0 0 0 0 0.17 0 0 0.17 0 0.56 0 0 0.56

0 0 0.29 0 0 0.58 0 0.29 0 0.29 0 0 0.29 0.29 0 0 0 0 0 0 0 0.08 0 0 0.08 0 0.56 0 0 0.56

(continued on next page) 13

Electrical Power and Energy Systems 118 (2020) 105779

Q. Sui, et al.

Table A5 (continued) Node

13h

14h

15h

16h

17h

18h

19h

20h

21h

22h

23h

24h

31 32 33 34 35 36 37 38 39 40 41

0.06 0 0 0.28 0 0.28 0.28 0 0 0 0.69

0.06 0 0 0.41 0 0.41 0.41 0 0 0 0.71

0.07 0 0 0.34 0 0.34 0.34 0 0 0 0.68

0.07 0 0 0.31 0 0.31 0.31 0 0 0 0.62

0.08 0 0 0.26 0 0.26 0.26 0 0 0 0.59

0.09 0 0 0.12 0 0.12 0.12 0 0 0 0.43

0.15 0 0 0.05 0 0.05 0.05 0 0 0 0.41

0.27 0 0 0.02 0 0.02 0.02 0 0 0 0.52

0.47 0 0 0.05 0 0.05 0.05 0 0 0 0.51

0.56 0 0 0.02 0 0.02 0.02 0 0 0 0.36

0.56 0 0 0.21 0 0.21 0.21 0 0 0 0.32

0.56 0 0 0.47 0 0.47 0.47 0 0 0 0.27

[23] Wei Wei, Liu Feng, Mei Shengwei. Distributionally robust co-optimization of energy and reserve dispatch. IEEE Trans Sustain Energy 2016;7(1):289–300. [24] Khodaei Amin, Bahramirad Shay, Shahidehpour Mohammad. Microgrid planning under uncertainty. IEEE Trans Power Syst 2015;30(5):2417–25. [25] Moghaddam EM, Nayeripour M, Aghaei J, Khodaei A, Waffenschmidt E. Interactive robust model for energy service providers integrating demand response programs in wholesale markets. IEEE Trans Smart Grid 2018;9(4):2681–90. [26] Karandeh Roozbeh, Lawanson Tumininu, Cecchi Valentina. Impact of operational decisions and size of battery energy storage systems on demand charge reduction. IEEE PES Power Tech 2019. [27] Naderi E, Pourakbari-Kasmaei M, Abdi H. An efficient particle swarm optimization algorithm to solve optimal power flow problem integrated with FACTS devices. Appl Soft Comput 2019;80:243–62. [28] Floudas Christodoulos A. Nonlinear and mixed-Integer optimization: fundamentals and applications. J Global Optim 1998;12(1):108–10. [29] Atwa YM. Distribution system planning and reliability assessment under high DG penetration. Waterloo, Ontario, Canada: University of Waterloo; 2010. [30] Lawder Matthew T, Suthar Bharatkumar, Northrop Paul WC, Sumitava De C, Hoff Michael. Battery energy storage system (BESS) and battery management system (BMS) for grid-scale applications. Proc IEEE 2014;102(6):1014–30. [31] Wei Fanrong, Sui Quan, Lin Xiangning. Energy control optimization strategy for isolated island microgrid with wind/photovoltaic/diesel/thermal storage under consideration of transferable load efficiency. Proc CSEE 2018;38(4). 10451053+1281. [32] Wang Yiyi. Research on optimal control strategy for multi-point charging / batteryswap station/energy storage system considering the power flow of distribution network. Beijing Jiaotong University; 2016. [33] L.A Traffic.[Online].Available: http://home.znet.com/schester/calculations/ traffic/la/plots/all/index.html.

References [1] Xie Shiwei, Zhijian Hu, Yang Li, Wang Jueying. Expansion planning of active distribution system considering multiple active network managements and the optimal load-shedding direction. Int J Electr Power Energy Syst 2020;115:1–20. [2] Quan S, Fanrong W, Zhang Rui. Xiangning L, Ning T, Zhixun Wang, Zhengtian Li.Optimal use of electric energy system for the building-integrated-photovoltaics community. Appl Energy 2019:1–11. [3] Wei Fanrong, Li Yuanzheng, Sui Quan, Lin Xiangning, Chen Le, Chen Zhe, et al. A novel thermal energy storage system in smart building based on phase change material. IEEE Trans Smart Grid 2019;10(3):2846–57. [4] Sun Yingyun, Li Zuyi, Shahidehpour Mohammad, Ai Bo. Battery-based energy storage transportation for enhancing power system economics and security. IEEE Trans Smart Grid 2015;2395–2402. [5] Sun Yingyun, Li Zuyi, Tian Wei, Shahidehpour Mohammad. A lagrangian decomposition approach to energy storage transportation scheduling in power system. IEEE Trans Power Syst. 2016;31(6):4348–56. [6] Sun Y, Zhong J, Li Z. Stochastic scheduling of battery-based energy storage transportation system with the penetration of wind power. IEEE Trans Sustain Energy 2016:135–44. [7] Chen Y, Zheng Y, Luo F, Wen J, Xu Z. Reliability evaluation of distribution systems with mobile energy storage systems. IET Renew Power Gener 2016;10(10):1562–9. [8] Yao S, Wang P, Zhao T. Transportable energy storage for more resilient distribution systems with multiple microgrids. IEEE Trans Smart Grid 2019;10(3):3331–41. [9] Lei S, Chen C, Li Y, Hou Y. Resilient disaster recovery logistics of distribution systems: co-optimize service restoration with repair crew and mobile power source dispatch. IEEE Trans Smart Grid 2019. [10] Lei S, Chen C, Z H, Hou Y. Routing and scheduling of mobile power sources for distribution system resilience enhancement. IEEE Trans Smart Grid 2018. [11] Abdeltawab HH, Mohamed ARI. Mobile energy storage scheduling and operation in active distribution system. IEEE Trans. Ind. Electron. 2017:1–10. [12] Kim Jip, Dvorkin Yury. Enhancing distribution system resilience with mobile energy storage and microgrids. IEEE Trans Smart Grid 2018. [13] Shah N, Phadke A, Waide P. Cooling the planet: opportunities for deployment of superefficient room air conditioners. In: Tech. Rep., Lawrence Berkeley Nat. Lab., Berkeley, CA, USA ; 2013. [14] Singh Sanjeev, Singh Bhim. A voltage-controlled PFC cuk converter-based PMBLDCM drive for air-conditioners. IEEE Trans Ind Appl 2012;48(2):832–8. [15] Gu B, Park JS, Choi J, Jung I. A design of PWM converter inverter system for building air conditioner compressor. International telecommunications energy conference. 2009. [16] J. K. VSS control of unity power factor. IEEE Trans Ind Electron 1999;46(2):325–32. [17] Kim Y-J, Norford LK, Kirtley JL. Modeling and analysis of a variable speed heat pump for frequency regulation through direct load control. IEEE Trans Power Syst 2015;30(1):397–408. [18] Kim Y-J, Fuentes E, Norford LK. Experimental study of grid frequency regulation ancillary service of a variable speed heat pump. IEEE Trans Power Syst 2016;31(4):3090–9. [19] Song Meng, Gao Ciwei, Yan Huaguang, Yang Jianlin. Thermal battery modeling of inverter air conditioning for demand response. IEEE Trans Smart Grid 2018;9(6):5522–34. [20] Ding X. Regulating strategy and effect evaluation of inverter air-conditioner applied in demand response. Southeast University; 2016. [21] Wang Beibei, Zhang Tianwei, Xiaoqing Hu, Bao Yuqing, Huiling Su. Consensus control strategy of an inverter air conditioning group for renewable energy integration based on the demand response. IET Renew Power Gener 2018;12(14):1633–9. [22] Zhou Anping, Yang Ming, Wang Zhaoyu, Li Peng. A linear solution method of generalized robust chance constrained real-time dispatch. IEEE Trans Power Syst 2018;33(6):7313–6.

Quan Sui received the B.Sc. degree in electrical engineering from Huazhong University of Science and Technology (HUST) in 2017. He is studying for a Ph.D at HUST. His research interests are microgrid planning and scheduling as well as optimal power system.

Fanrong Wei received the Ph.D. degree in electrical engineering, Huazhong University of Science and Technology (HUST). His researches mainly focus on optimal power system/microgrid scheduling and protective relay.

14

Electrical Power and Energy Systems 118 (2020) 105779

Q. Sui, et al. Xiangning Lin is currently a Professor at Huazhong University of Science and Technology (HUST). He received the M.Sc. and Ph.D. degrees in electrical engineering from HUST, Wuhan, China. His research interests are modern signal processing and power system protective relaying.

Zhixun Wang received the B.Sc. degree in electrical engineering from Huazhong University of Science and Technology (HUST) in 2016, where he is currently pursuing the Ph.D. degree in electrical power engineering. His researches mainly focus on micro-gird planning and optimization.

Chuantao Wu received the B.Sc. degree in gineering from Huazhong University of Technology (HUST) in 2018. He is studying HUST. His research interests are microgrid scheduling as well as optimal power system.

Zhengtian Li received the B.Sc. degree from Wuhan University in 2002 and the Ph.D. degree from the Huazhong University of Science and Technology (HUST) in 2011. Currently, he is an Associate Professor at HUST. His research interests are digital protection relaying.

electrical enScience and for a Ph.D at planning and

15