Decentralized saddle-point dynamics solution for optimal power flow of distribution systems with multi-microgrids

Applied Energy 252 (2019) 113361

Contents lists available at ScienceDirect

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

Decentralized saddle-point dynamics solution for optimal power flow of distribution systems with multi-microgrids

T

Chong Tang, Mingbo Liu , Yue Dai, Zhijun Wang, Min Xie ⁎

School of Electric Power Engineering, South China University of Technology, Guangzhou 510640, China

HIGHLIGHTS

optimization model for multi-microgrid systems is established. • AA decentralized linearized Ward equivalent method for optimal power flow is proposed. • Thenovel • decentralized saddle-point dynamics method is applied to solve the problem. ARTICLE INFO

ABSTRACT

Keywords: Distribution system Multi-microgrids Optimal power flow Decentralized saddle-point dynamics method Quadratic programming Ward equivalent

Since various types of distributed renewable energy resources are integrated into distribution systems in the form of microgrids, how to implement the coordinated operation between a distribution system and microgrids is a major concern. This paper focuses on solving the optimal power flow of a distribution system with multiple microgrids based on the decentralized saddle-point dynamics approach. First, the distribution system and the microgrids are regarded as separate entities and their external networks are replaced by Ward equivalent circuits. Hence, the distribution system and microgrids are decoupled so that their linearized power flow models can be built separately. Then, a decentralized quadratic programming model for optimal power flow with the distribution system and microgrids as separate entities is established to achieve low active power loss and high utilization of renewable energy resources. Next, the decentralized saddle-point dynamics approach is applied to solve this model from the viewpoint of the dynamic system control, which transforms the solution of KarshKuhn-Tucker conditions into an asymptotically stable process. The proposed method only needs to exchange the border-bus voltages and equivalent injection powers between the distribution system and each microgrid, which can protect the privacy of different entities and possesses plug-and-play features. Finally, case studies on a real distribution system with two real microgrids are carried out to verify the effectiveness of the proposed method.

1. Introduction Microgrids (MGs) are active clusters of renewable energy resources (RERs), loads, energy storage devices, and other onsite electric components. MGs are typically linked to the downstream section of the medium and low voltage distribution system (DS) [1]. MGs can be considered as intelligent distribution networks with two different modes of operation: islanded mode for the local production of electricity and grid-connected mode. When operating in the grid-connected mode, MGs can sell/buy electricity to/from the DS if necessary [2]. Meanwhile, various types of MGs may exist in the future smart distribution system, where the distribution system owner (DSO) and microgrid owners (MGOs) are usually regarded as independent autonomous entities. Therefore, how to effectively coordinate the optimization ⁎

among networked MGs and a DS is worth further investigating in the future [3]. To date, several works on coordinated optimization between DSs and MGs have been reported. In [4], a two-stage co-optimization model with dynamic electricity pricing was proposed as an incentive scheme for privately-owned MGs to participate in a DS's operation. Fathi et al. considered a distribution network connecting to several MGs in a local area with deterministic demand and other neighbor areas with uncertain demands [5]. An adaptive energy consumption scheduling approach was applied to achieve low-power generation cost and a low peak-to-average ratio. Fang et al. proposed a cooperative operation control strategy in a radial network consisting of multiple autonomous MGs [6]. The cooperation of multiple MGs can be transformed into the problem of optimally allocating the DGs and the critical loads. In [7], a

Corresponding author. E-mail address: [email protected] (M. Liu).

https://doi.org/10.1016/j.apenergy.2019.113361 Received 10 August 2018; Received in revised form 24 January 2019; Accepted 16 May 2019 0306-2619/ © 2019 Elsevier Ltd. All rights reserved.

Applied Energy 252 (2019) 113361

C. Tang, et al.

Nomenclature

e'i, f'i Psell(m)

Abbreviation HV DS MG RER DSO MGO DG PV WT BG ESS SVG OPF ACOPF DSPD ADMM QP KKT

Pbuy(m) Pj/Qj

P lj /Qjl

high voltage distribution system microgrid renewable energy resources distribution system operator MG owner distributed generator photovoltaic wind turbine bio-gas generator energy storage system static Var generator optimal power flow alternating current optimal power flow decentralized saddle-point dynamics alternating direction method of multipliers quadratic programming Karush-Kuhn-Tucker

P dg predicted active power of RER j P jc curtailed active power of RER Qjc reactive power output of inverter of RER m Acod, m , Aco /P dis charging/discharging power of ESS j charging/discharging efficiency of ESS c/ d reactive power output of SVG Qjsvg Wj total energy storage of ESS d d Peq , r / Qeq, r equivalent injection power of DS at the DS side Peqd, r / Qeqd, r equivalent injection power of DS at the MG side m m Peq , s / Qeq, s equivalent injection power of MG at the DS side m Peq, s/ Qeqm, s equivalent injection power of MG at the MG side Sets and indices real and imaginary parts of the impedance matrix for DS with MGs as equivalent networks Rm / X m real and imaginary parts of the impedance matrix for MG with DS and other MGs as equivalent networks DL set of buses in DS BD set of border buses in DS BM set of border buses in MGs PD/PM active power of buses in the original network of DS/MG QD/QM reactive power of buses in the original network of DS/MG M(m) set of buses in the mth MG G(m) set of DG buses in the mth MG S(m) set of SVG buses in the mth MG E(m) set of ESS buses in the mth MG xd/xm primal variables for the subproblem of DS/MG μd, λd/μm, λm dual variables for the subproblem of DS/MG Gd, rd/Gm, rm quadratic and linear term constants of objective function for the subproblem of DS/MG d d m m Aeq , beq /Aeq , beq coefficient matrix and constant term vector of equality constraints for the subproblem of DS/MG Aind , bind/ Ainm , binm coefficient matrix and constant term vector of inequality constraints for the subproblem of DS/MG m Acod, m , Aco coefficient matrix of coupling constraints

R d /X d

Parameters ρloss ρrer ρtrans

weighting constants for active power loss weighting constant for curtailed active power cost of RERs weighting constant for cost/benefit of electricity transaction closs price for active power losses price for selling/buying electricity to/from MGs ctrans ar, br penalty parameters for curtailed active power ijth element of admittance matrix yij Vmin/Vmax lower/upper bounds of voltage magnitude power factor angle of inverter j svg Qmax upper bound of SVG output P jch, rate /P dis j, rate rated charging/discharging power Wj,0 total energy storage of ESS at the initial time Wj,min/Wj,max lower/upper bound of the total energy storage of ESS Variables ei, fi

real and imaginary parts of voltage in MG selling power from DS to the mth MG buying power from the mth MG bus injection active/reactive power active/reactive power of load

real and imaginary parts of voltage in DS

bi-level energy management strategy considering the interaction between a DS and MGs was presented, where an interactive energy game matrix was calculated to guide the interaction among different levels. In [8], a novel bi-level optimization approach for the self-healing resilient operation of distribution networks was presented. Considering a specific fault, the first layer determines the optimal configuration of the distribution network and the second layer determines the optimal operation of the DGs. Du et al. proposed a cooperative game approach to implement multi-MG coordinated operation at the distribution system level [9]. In a typical MG, DGs with inverters, such as photovoltaic (PV) generators and wind turbines (WT), are usually installed. The amount of reactive power injected/consumed by the inverters can be set by drop control [10] or optimal power flow (OPF) [11], which is the focus of this paper. In the medium and low voltage distribution network, line resistance and reactive power flows should not be neglected. Therefore, the traditional direct current optimal power flow model is not suitable in the study of the OPF problem in distribution networks, and the alternating current optimal power flow (ACOPF) model considering constraints on node voltage and branch power should be established,

which is a nonlinear and non-convex optimization problem including nonlinear equality and inequality constraints [12]. To solve the ACOPF model of distribution networks, many related works have been reported. Instead of the traditional alternating current power flow model, a balanced radial network can be represented using the DistFlow model [13]. Since the original DistFlow model contains non-convex constraints, a second-order cone programming model was proposed [14]. In [15], another common relaxation technique, i.e., semidefinite programming relaxation, was proposed. The second-order cone programming relaxation is much simpler computationally than the semidefinite programming relaxation, while the semidefinite programming relaxation is tighter for general networks [16]. Unlike the linearized model in direct current OPF, an alternative linearization approach was applied in the ACOPF of distribution networks under the assumptions that shunt impedances are negligible and voltage magnitudes are represented near nominal values [17]. However, the voltage magnitude and phase angle of each bus correspond to the injection powers at all buses. Thus, it is difficult to apply directly in a decentralized coordinated optimization between DSs and MGs. Considering different benefits and optimization targets, a 2

Applied Energy 252 (2019) 113361

C. Tang, et al.

decentralized scheduling model tends to be more reasonable in the practical application [18]. A decentralized solution to the coordinated dispatch of integrated electrical and heating systems was proposed with guaranteed convergence for convex programs [19]. In [20], a multiagent system of a stand-alone MG was built, and it was capable of learning to control different components of the MG in a distributed model. Liu et al. proposed a decentralized optimization model for the coordination of transmission and distribution networks [21]. In [22], a fully distributed voltage control method in DS was proposed, and each component was regarded as an agent. Furthermore, since the privacy of data will become increasingly important in the future, decentralized algorithms have potential advantages over centralized algorithms in a coordinated scenario of a DS with multiple MGs. In [23], decentralized optimal voltage regulation for DSs with multiple MGs was established, which can execute real-time control and reduce the calculation and communication burden. In [24], an autonomous optimized dynamic economic dispatch model was established considering the DS and MGs as different interest entities. In [25], a decentralized energy management framework was proposed to coordinate the power exchange between a DS and MGs in a fully decentralized mode. In decentralized convex optimization research, the alternating direction method of multipliers (ADMM) algorithm and its variants have been extensively applied. In [26], a fully distributed algorithm for reactive power optimization based on the ADMM algorithm was proposed, which transformed the non-convex power flow equations into convex equations utilizing the second-order cone programming relaxation. The penalty factor regulation technique, which can accelerate convergence, was also presented. In [27], a distributed OPF method for radial networks based on an improved ADMM was proposed, which can greatly speed up each iteration. Communication in this distributed algorithm is only required between adjacent buses. In addition to ADMM, other decentralized optimization techniques have been investigated for OPF. A bi-level analytical target cascading algorithm was applied in [28] to coordinate the operation of a DS with multiple MGs. A parallel OPF approach for large multi-regional interconnected power systems was presented in [29], which provided guidance regarding the choice of penalty parameters in the auxiliary problem principle formulation. However, the aforementioned decentralized optimization techniques usually require tolerances as the stop criteria of the algorithm, which may greatly influence their calculation precision and time [30]. The saddle-point dynamics method transforms the traditional optimization problem into a gradual stabilization process based on Gradient Dynamics programming [31]. In [32], the saddle-point dynamics method was applied to solve the nonconvex OPF problem without a guarantee that the dynamic system can converge to a global optimal solution. Besides, the information exchanged among units is high, which may greatly slow its convergence speed. In [33], a discontinuous decentralized saddle-point dynamics (DSPD) model was established to solve a linear programming problem distributedly. It has been reported that this dynamic model can converge to the global optimal solution in the case of a sudden disturbance. This paper studies a novel decentralized OPF optimization method in a DS with multiple MGs. In our study, the DSO and MGOs are regarded as separate entities with different interest and privacy. Power exchange occurs between the DS and each MG to maximize their own profits. We used a linearized AC power flow model and Ward equivalent model to build a decentralized QP model for the OPF of a DS with multiple MGs. The proposed decentralized optimization was implemented through the DSPD approach. The main contributions of this paper are as follows:

Fig. 1. Framework of DS with two networked MGs.

global optimal solution and also solves the problem of the inability to decouple the DS and MGs in the current linearized power flow model. (2) For the first time, DSPD is applied to solve a decentralized OPF of a DS with multiple MGs, which can obtain the global optimum solution and possesses plug-and-play functions considering network reconfiguration, bus power variation, and MGs plugging out. The rest of this paper is organized as follows. Section 2 depicts the linearized AC power flow models with Ward equivalent for the DS and MGs. A decentralized optimization model of the DS with multiple MGs is established in Section 3, and DSPD for solving this proposed model is introduced in Section 4. Case studies, results, and discussion are presented in Section 5. Section 6 concludes this paper. 2. Linearized AC power flow model with Ward equivalent for DS and MGs As a typical example, two networked MGs are integrated into a DS in Fig. 1. In our model, the DS is directly connected to the HV system, and MGs are linked to the DS. The main objective of the DS and each MG is to minimize their own operation costs. Dispatchable DGs, such as bio-gas generators (BGs), PVs, WTs, and energy storage systems (ESSs), can also be deployed in MGs based on their own planning schemes. In the linearized AC power flow model [17], the voltage magnitude and phase angle can be approximated as expressions of bus injection powers and impedance matrix. Since the impedance matrix is a full matrix, the voltage magnitude and phase angle of each bus correspond to the injection powers at all buses. However, for an independent entity, the information about local bus injection powers and impedance should not be obtained by other entities to protect privacy. Therefore, it is not feasible to apply the linearized AC power flow model directly for the decentralized optimization of multi-MG DSs. To solve this problem, the Ward equivalent method [34] is introduced here to decouple the DS and MGs. 2.1. DS model with MGs as equivalent networks The DS with two MGs as equivalent networks is illustrated in Fig. 2. By retaining the root bus, the network of the DS remains invariant while the networks of MG1 and MG2 are reduced to two border buses (i.e., BM1 and BM2) with the corresponding equivalent injection powers. It is noted that the Ward equivalent depicted in Fig. 2 is true under assumptions that shunt impedances can be negligible in the DS and MGs. On the basis of Ward equivalent, the original networks of the DS and MGs are transformed to novel networks with the root bus remaining and shunt impedances negligible. Then, the voltage magnitudes and phase angles of buses in the DS together with BM1 and BM2 can be approximated as functions of the bus injection active and reactive powers as follows [17]:

(1) A novel linearized Ward equivalent method for ACOPF problems is proposed. It can transform a conventional nonlinear ACOPF model into a distributed QP model. The proposed model avoids the problem that the conventional ACOPF model cannot guarantee the 3

Applied Energy 252 (2019) 113361

C. Tang, et al.

3. Decentralized QP model of DS with multi-MGs In this section, a decentralized OPF model for the DS with multiple MGs is introduced from three respects: QP subproblem of the DS, QP subproblem of MGs, and coupling constraints between QP subproblems of the DS and MGs. 3.1. QP subproblem of DS A. Objective The objective of the DS is to minimize the active power loss and maximize the revenue of supplying power to MGs, that is, DS loss · Closs

min OFDS =

Fig. 2. DS with MGs as equivalent networks.

DS Closs = closs·

|V |i

(Rijd Pj

1+

+

Xijd Qj ),

i

DL

i

Rijd Qj ),

i

DL

BM .

(1a)

BM .

(1b)

i

M (m )

j Mm

(Xijm Pj

i j Mm

Rijm Qj ),

i

M (m )

BD .

(3b)

ctrans · Psell (m)

BD .

(3c)

er2

Psell (m) = (er es + fr fs

r

fr2 )·Re{yrs } + (er fs BD , s

fr es )·Im{yrs }

BM .

(3d)

B. Constraints (1) Bus voltage constraints

The MG with the DS and other MGs as equivalent networks can be illustrated in Fig. 3. By retaining the root bus in the DS, the network of MG2 remains invariant while the networks of the DS and MG2 are reduced to a border bus (i.e., BD2) with the corresponding equivalent injection power. To protect the privacy of different entities, the detailed network and bus injection power information of MG1 should not be known by the owner of the DS and MG2. Therefore, the network of MG1 is reduced to a border bus (i.e., BM1) with an equivalent injection power firstly. Then, the network of the DS and BM1 with the corresponding injection power is reduced to an equivalent impedance with a novel equivalent injection power at BD2. It is noted that the Ward equivalent depicted in Fig. 3 is true under the assumptions that shunt impedances can be negligible in the DS and MGs. On the basis of the Ward equivalent, the original networks of the DS and MGs are converted to novel networks with the root bus remaining and shunt impedances negligible. The voltage magnitudes and phase angles of buses in MG2 together with BD2 can be approximated as functions of the bus injection active and reactive powers as follows:

(Rijm Pj + Xijm Qj ),

f j )2)

m=1

2.2. MG model with DS and other MGs as equivalent networks

1+

ej )2 + (fi

Re{yij }·((ei

M DS Bsell =

j DL

|V |i

(3a)

(i, j ) D

j DL

(Xijd Pj

DS trans · Bsell

(Rijd Pj + Xijd Qj ),

ei = 1 +

i, j

DL

BM ,

(4a)

j

(Xijd Pj

fi =

Rijd Qj ),

i, j

DL

BM ,

(4b)

j

Vmin

(Rijd Pj + Xijd Qj )

1+

Vmax

i

j

DL

BM ,

(4c)

where (4a) and (4b) capture the voltage across the network [23], and (4c) represents the constraint of voltage magnitudes according to (1a). (2) Bus injection power balance constraints m Ps = Peq ,s,

s

BM

(5a)

Qs =

m Qeq ,s,

Pj =

P lj,

j

DL

(5c)

(2a)

Qj =

Qlj ,

j

DL .

(5d)

(2b)

(3) Equivalent injection power constraints

s

BM

Fig. 3. MG with DS and other MGs as equivalent networks. 4

(5b)

Applied Energy 252 (2019) 113361

C. Tang, et al. ds = ((R d )2 + (X d )2) 1·(R d R d + X d X d )· P + ((R d ) 2 + (X d ) 2) 1·(R d X d Peq D ,r rr rr rr rD rr rD rr rr rr rD d )· Q , Xrrd RrD D

r

BD

ds d 2 d 2 1 d d Qeq , r = ((Rrr ) + (Xrr ) ) ·(Xrr RrD

(6a)

d Rrrd XrD )·PD

Pj = P jdg

P jc

Qj = Qjc

Qlj ,

P lj,

j

j

(9a)

G (m )

(9b)

G (m )

Pj = P jdis

P jch

where (6a) and (6b) denote the detailed expression of equivalent injection powers of the DS side at the border buses. Their derivation is detailed in the Appendix A.

Qj = Qjsvg

Qlj ,

3.2. QP subproblem of MGs

Qr = Qeqd, r ,

A. Objective

(3) DG and SVG power output constraints

The objective of the MG is to minimize the active power loss, the cost of charging and discharging [35], and the cost of purchasing electricity from the DS and to maximize RER utilization, that is,

0

d d + ((Rrrd )2 + (Xrrd ) 2) 1·(Rrrd RrD + Xrrd XrD )·QD ,

min OF MG (m) =

MG loss · Closs (m)

MG Closs (m) = closs·

MG rer · Crer (m )

+

+

ej )2 + (fi

Re{yij }·((ei (i, j) M (m)

MG Crer (m ) =

MG ess · Cess (m)

r

+

BD ,

MG trans · Cbuy (m)

f j ) 2)

(7a)

MG Cbuy (m ) = ctrans· Pbuy (m),

r

(fr ) 2)·Re{yrs } + (fr es BD , s

er f s )·Im{yrs }

BM

i, j

M (m )

Vmin

(Rijm Pj + Xijm Qj )

j

(10c)

S (m ).

min P jch, rate,

1 t

P jdis

min P jdis , rate ,

d

t

(Wj,max

(Wj,0

Wj,0) ,

Wj,min ) ,

j

j

E (m )

(11a)

E (m )

(11b)

i, j

M (m )

Vmax ,

j

BD

i

M (m )

m m = ((Rssm ) 2 + (Xssm ) 2) 1·(Rssm RsM + Xssm XsM )· PM + ((Rssm ) 2 + (Xssm )2)

(8a)

m ·(Rssm XsM

(8b)

j

1+

svg Qmax ,

(10b)

G (m )

Peqm, s

BD

j

Rijm Qj ),

j

(5) Equivalent injection power constraints

(Rijm Pj + Xijm Qj ), (Xijm Pj

P jc ),

(10a)

G (m )

As depicted in (11a)–(11b), a common battery model is applied to describe the charging/discharging power limit and energy balance of ESS [36]. Meanwhile, it can be proved that the optimal pair = 0, (Pich , Pidis ) must satisfy the complementary constraint P jch ·P dis j given an increasing positive operational cost function for charging and discharging in (7d). The detailed proof can be seen in Appendix B.

B. Constraints (1) Bus voltage constraints

fi =

(9f)

BD .

tan j ·(P jdg

P jch

0

(7f)

ei = 1 +

(9e)

BD

r

Qjs

(9d)

S (m )

c

(7e)

(er ) 2

r

(9c)

E (m )

(4) ESS constraints

(7d)

i E (m )

j

j

svg Qmin

(7c)

(ach (Pich ) 2 + adis (Pidis )2).

Pbuy (m ) = (er es + fr f s

|Q jc|

j

P jdg ,

P jc

0

i G (m )

P lj,

Pr = Peqd, r ,

(7b)

(ar (Pic )2 + br Pic )

MG Cess (m ) =

(6b)

BD .

m Xssm RsM )·QM ,

s

m Qeqm, s = ((Rssm )2 + (Xssm ) 2) 1·(Xssm RsM

(12a)

BM m Rssm XsM )·PM

m m + ((Rssm )2 + (Xssm )2) 1·(Rssm RsM + Xssm XsM )·QM ,

(8c)

1

s

BM . (12b)

(2) Bus injection power balance constraints

Fig. 4. Exchange process of coupling variables for DS and MGs. 5

Applied Energy 252 (2019) 113361

C. Tang, et al.

3.3. Coupling constraints between QP subproblems of DS and MGs

4. Implement of the DSPD approach

er = er ,

r

BD

(13a)

fr = fr ,

r

BD

(13b)

es = es ,

s

BM

(13c)

fs = f s ,

s

BM

(13d)

m Peq ,s m Qeq ,s

d Peq ,r

=

Peqm, s,

s

BM

(13e)

=

Qeqm, s ,

s

BM

=

Peqd, r ,

r

BD

d d Qeq , r = Qeq, r ,

r

BD .

To simplify the discussion on the decentralized solution, we rewrite the decentralized QP model in Section 3 as follows: M OF MG (x m) m=1

min f = OF DS (xd) +

x d , xm 1

= 2 xdT Gd xd + rdT xd + d Aeq xd

s. t

Aind xd

M m=1

(

1 T x G x 2 m m m

d beq =0

bind

+ rmT xm

)

(14a) (14b)

0

(14c)

(13f)

m Aeq xm

m beq

=0

(14d)

(13g)

Ainm xm

binm

0

(14e)

(13h)

Acod, m xd

m Aco xm = 0,

(14f)

where (14a) denotes the operation cost of the DS and MGs, (14b)–(14c) denote the equality and inequality constraints in the DS, (14d)–(14e) denote the equality and inequality constraints in the MGs, and (14f) denotes the coupling constraints between the DS and MGs. In this paper, we adopt the DSPD approach to solve the decentralized QP model in (14). First, we build the augmented Lagrange function L. Then, we transform the solution of the decentralized QP model in (14) into an asymptotically stable process of the KKT conditions, which can be obtained according to the augmented Lagrange function. Thus, we construct the corresponding differential equations based on the SPD method as follows:

The coupling constraints between the DS and MGs are expressed in (13a)–(13h), where the entries on the right of the equations represent the information of voltages and equivalent injection powers at the border buses (i.e., BD and BM) transmitted by MGs, and the entries on the left represent the information of voltages and equivalent injection powers at the border buses transmitted by the DS. According to Figs. 2 and 3, the exchange process of coupling variables for QP subproblems of the DS and MGs are illustrated in Fig. 4. The border variables consist of the voltages and the equivalent injection powers at the border buses (i.e., BD2 and BM2). In the subproblem of the DS, the equivalent injection power at BD2 is calculated based on (6a)–(6b), which is only transmitted to MG2 and not involved in the optimization during each iteration. The equivalent injection power at BM2 is optimized for the subproblem of the DS with the coupling constraints in (13a)–(13d) and (13e)–(13f). In the subproblem of MG2, the equivalent injection power at BM2 is calculated based on (12a)–(12b), which is only transmitted to the DS and not involved in the optimization during each iteration. The equivalent injection power at BD2 is optimized for the subproblem of MG2 with the coupling constraints in (13a)–(13d) and (13g)–(13h).

L xd

xd = µd =

L µd L

d

=

d

(16)

,

0

Fig. 5. Flowchart of DSPD for OPF. 6

(15)

d

> 0 or Aind xd otherwise .

bind > 0 (17)

Applied Energy 252 (2019) 113361

C. Tang, et al.

L xm

xm = µm =

L µm L

m

=

m

generation power of renewable energy at 11:00 a.m. in a day in a city of China. To verify the convergence of the proposed model, Fig. 9 shows parts of dynamic trajectories of the bus voltage magnitudes in the DS. As depicted in Fig. 9, the trajectories of voltage magnitudes start out getting gradually close to 1p.u., then oscillate near 1 p.u., and finally become stable in certain values. In other words, after a certain number of iterations, the dynamic trajectories of the voltage magnitudes become smooth lines. Similar to the trajectories of voltage magnitudes, the trajectories of power outputs tend to be smooth lines after a certain amount of oscillation. Figs. 10 and 11 respectively show the trajectories of active and reactive power outputs in MG1 and MG2. Fig. 12 shows the dynamic trajectories of powers exchanged between the DS side and the MG side. The interactive power of the DS side is not equal to that of the MG side at the beginning. As the information exchange proceeds in the dynamic system, the trajectories will converge to the agreement on both the DS side and the MG side.

(18) (19)

,

m

0

µco, m = Acod

> 0 or Ainm xm

binm > 0

otherwise . m

xd

Acom xm

(20) (21)

where Eqs. (15)–(17) denote the dynamic process of primal and dual variables in the DS, while Eqs. (18)–(20) denote the dynamic process of primal and dual variables in the MGs. Eq. (21) shows the dynamic process of dual variables in coupling constraints. The detailed calculation equations can be seen in Appendix C. Eqs. (17) and (20) ensure that the complementary slackness constraints of the KKT conditions can be satisfied during the whole iterative process if the initial values of λd and λm are non-negative. According to [32], a saddle point must satisfy the KKT conditions. Since the decentralized QP model in (14) is convex, a saddle point is exactly a global optimal solution. As can be seen from (15)–(21), the update of variables of the DS and MGs is only related to their own primal variables and dual variables, and MGs are completely independent of each other and completely confidential to each other. In each iteration, the DS receives the information of voltages and equivalent injection powers at the border buses from the MGs, updates its local bus voltages and bus injection powers, and then sends the information of voltages and equivalent injection powers at the border buses to each MG. Each MG receives and sends the information in a similar way. The detailed calculation procedure of DSPD is illustrated in Fig. 5.

B. Accuracy analysis of DSPD for decentralized OPF To verify the accuracy of the proposed model, a comparison of voltage magnitudes and angles is presented in Fig. 13. It can be found that the results obtained by the DSPD algorithm are almost the same as those obtained by the centralized algorithm applying CPLEX solver. Furthermore, as depicted in Fig. 13, the voltage magnitudes of buses in MG1 (buses 36 to 43) achieve approximate values mainly because of the single bus radiation structure of MG1. The voltage magnitudes of buses in MG2 (buses 44 to 56) obtain different values because most buses in MG2 are not connected to a public bus. To evaluate the performance of the proposed method from another perspective, we introduce one measure that determines the relative error of the objective function compared with the optimal value calculated in the centralized algorithm over the iterations. The detailed expression can be given as follows:

5. Numerical results and discussion 5.1. 35-bus system with 2 MGs

rel =

The simulation was carried out on a real distribution system connecting two MGs. This distribution system consisted of 35 buses and 37 lines, three of which were contact lines (the green dotted lines in Fig. 6). The reference voltage in the headend of network was 10.5 kV, and the basis power was 100 MW. As depicted in Fig. 6, MGs were located at buses 12 and 33. The structures of MG1 and MG2 can be seen in Figs. 7 and 8, respectively. MG1 consisted of three BGs, one SVG, one PV, one ESS, and one load, while MG2 consisted of one BG, one WT, two PVs, two ESSs, and four loads. Relevant data of partial electric components is depicted in Table 1 and Table 2. Other parameters for the proposed optimization problem are as follows: ρloss = 1, ρrer = 1, ρess = 1, ρtrans = 1, ach = 0.1, adis = 0.15, ar = 0.1, br = 1 [17], closs = 2 $/kWh [37], ctrans = 0.3 $/kW [3], ηc = 0.95, and ηd = 0.95 [36]. The dynamic optimizaion system illustrated in Fig. 5 was established in a MATLAB/Simulink platform utilizing ode15s solver to calculate the simulation results. The computer was equipped with an Intel (R) Xeon(R) E3-1225 v5 CPU with 16G of RAM. The numerical results can be analyzed from the following aspects: convergence analysis of DSDP for decentralized OPF, accuracy analysis of DSPD for decentralized OPF, and plug-and-play features analysis of DSDP for decentralized OPF.

|f

f | (22)

f

where f* is the optimal objective function value calculated by solving the centralized OPF problem. Fig. 14 depicts the trajectories of the relative error to the centralized solution. It is obvious that the relative error will converge to a stable value less than 10-6, which validates the accuracy of the algorithm effectively. Additionally, the calculation time of the proposed algorithm is only about 9.5 s, which can fully meet the actual engineering requirements. More details of calculation results can be found in Table 4. C. Plug-and-play features analysis of DSPD for decentralized OPF

A. Convergence analysis of DSPD for decentralized OPF In this paper, the system load and the predicted power generated from renewable energy were selected as input parameters. Table 3 lists the distribution network load, microgrid load, and the predicted

Fig. 6. Initial structure of a real DS with two real MGs. 7

Applied Energy 252 (2019) 113361

C. Tang, et al.

Table 3 Parameters of load demand and renewable generation power at 11:00 a.m.

Load in DS Load in MG1 Load in MG2 PV in MG1 PV1 in MG2 PV2 in MG2 WT in MG2

Active power/kW

Power factor

3112 1525 503 143 541 28 89

0.97 0.80 0.80 0.85 0.85 0.85 0.85

Fig. 7. Structure of MG1.

Fig. 9. Trajectories of voltage magnitudes.

18 and 20, and buses 20 and 21 are disconnected. Fig. 16 shows the trajectories of voltage magnitude before and after network reconfiguration. When the network information is updated, the trajectories of voltage magnitudes change from smooth straight lines to oscillating curves. Then, after a small number of iterations, the trajectories of voltage magnitudes stabilize, forming new straight lines. Furthermore, Fig. 16 shows that voltage magnitudes at buses 19, 20, and 21 tend to achieve lower values than before, which is mainly because those buses are closer to the end of the network after network reconfiguration. From the preceding discussion, it can be concluded that when the network structure in the DS changes during the dynamic optimization process, the proposed dynamic optimization problem will still converge to a novel optimal solution. Besides, the new convergence process takes fewer iteration steps than the initial convergence process, which indicates that the number of iterations can be small when the initial value is close to the optimal value. Scenario 2: DSPD for OPF considering bus power variation Table 5 lists the distribution network load value, microgrid load value, and the predicted generation power of renewable energy at 11:15 a.m. As depicted in Figs. 17 and 18, after bus power information is updated, the trajectories of power outputs in MGs change from smooth straight lines to oscillating curves. Then, after a small number of iterations, the trajectories stabilize as new straight lines. Since the parameters of bus power have changed, optimal power outputs of each electric components may also change. In Fig. 17, the active power outputs of PV and ESS and the reactive power outputs of PV achieve new values and other electric components remain invariant when reaching the stable state. In Fig. 18, the active power and the reactive power outputs of PV1, PV2, ESS and WT achieve new values and other electric components remain invariant when reaching the stable state. Similar to Scenario 1, it can be concluded that when bus power changes during the dynamic optimization process, the proposed dynamic optimization problem will still converge to a novel optimal solution. Scenario 3: DSPD for OPF considering one MG plugs out

Fig. 8. Structure of MG2. Table 1 Parameters of partial electric components for MG1. BG1

Rated active power/kW 495

Power factor 0.9

BG3

BG2

Rated active power/kW 495

Power factor 0.9

ESS

W0/kWh

Wmin/ kWh 200

500

Rated active power/kW 495

Power factor

SVG

Qmin/kVar

Qmax/kVar

−420

420

Wmax/kWh

ch Pmax /kW

dis Pmax /kW

1000

250

250

0.9

Table 2 Parameters of partial electric components for MG2. BG ESS1 ESS2

Rated active power/kW 600

Power factor 0.9

W0/kWh

Wmin/kWh

Wmax/kWh

650

260

1300

ch Pmax /kW 300

dis Pmax /kW 300

W0/kWh

Wmin/kWh

Wmax/kWh

250

100

500

ch Pmax /kW 150

dis Pmax /kW 150

Scenario 1: DSPD for OPF considering network structure reconfiguration As depicted in Fig. 6, there are several contact lines in DS for network reconfiguration. Fig. 15 shows the reconfiguration structure of the DS. The contact lines in a closed state are represented by a solid green line. Bus 19 is connected to buses 11 and 21. Bus 21 is connected to buses 19 and 26. Meanwhile, the lines between buses 18 and 19, buses 8

Applied Energy 252 (2019) 113361

C. Tang, et al.

Fig. 10. Trajectories of power outputs in MG1.

Fig. 11. Trajectories of power outputs in MG2.

Fig. 12. Trajectories of powers exchanged between DS and MGs.

Fig. 13. Comparison of voltage magnitudes and angles in centralized and decentralized methods for 35-bus system with two MGs.

When one or more MGs plug out from DS, the topology of communication will undoubtedly be changed. Fig. 19 illustrates the selfhealing feature of DSPD. When MG2 plugs out because of failure or communication interruption during the dynamic optimization process, power exchange and communication between the DS and MG2 are no longer carried out. DS will no longer receive the information of voltage

and equivalent injection power at the border buses from the MG2. However, the communication between the DS and MG1 is still linked. Therefore, the DS and MG1 continue to update their local bus voltage and bus injection power until convergence. As shown in Fig. 20, the optimization algorithm will converge to a 9

Applied Energy 252 (2019) 113361

C. Tang, et al.

Fig. 16. Trajectories of voltage magnitudes before and after network reconfiguration.

Fig. 14. Trajectories of relative error to centralized solution for 35-bus system with two MGs.

Table 5 Parameters of load demand and renewable generation power at 11:15 a.m.

Table 4 Comparison of decentralized and centralized methods for 35-bus system with two MGs. Method

Entity

Operation Cost ($)

Exchange power cost ($)

Decentralized (DSPD)

DS MG1 MG2 Total

60.3452 8.0839 11.6925 80.1217

From DS to MG1

−93.6464

From DS to MG2

−318.4142

Centralized (CPLEX)

DS MG1 MG2 Total

60.3452 8.0839 11.6925 80.1217

From DS to MG1

−93.6464

From DS to MG2

−318.4142

Load in DS Load in MG1 Load in MG2 PV in MG1 PV1 in MG2 PV2 in MG2 WT in MG2

Active power/kW

Power factor

3504 1688 545 305 428 22 200

0.95 0.80 0.80 0.85 0.85 0.85 0.85

Fig. 17. Trajectories of active power outputs in MG1 considering bus power variation.

Fig. 15. Reconfiguration structure of a real DS with two real MGs.

new global optimal solution after MG2 plugs out. When the signal of plugging out starts, the trajectory of power that the DS supplies to MG2 approaches zero while the trajectory of power that the DS supplies to MG1 tends to be stable again after a slight oscillation. In general, it can be concluded that when one MG plugs out, the dynamic optimization problem will still converge to a novel optimal solution. Therefore, the proposed method in this article possesses selfhealing features considering failure or communication interruption. 5.2. 110-bus system with four MGs To verify the scalability of the proposed method, we applied the proposed method to a large-scale practical 110-bus active DS in China, and its topology is presented in Fig. 21. This DS contained 109 branches and 110 buses. The total load power is 10.35 + 2.27i MVA. The

Fig. 18. Trajectories of active power outputs in MG2 considering bus power variation.

10

Applied Energy 252 (2019) 113361

C. Tang, et al.

and relevant data of MG4 was the same as that of MG2. The nodes in MG1 are numbered from 111 to 118, the nodes in MG2 are numbered from 119 to 131, the nodes in MG2 are numbered from 132 to 139 and the nodes in MG2 are numbered from 140 to 152. A. Results analysis for 110-bus system with four MGs Fig. 22 shows the dynamic trajectories of powers exchanged between the DS side and the MG side. The interactive power of the DS side is not equal to that of the MG side at the beginning. As the information exchange proceeds in the dynamic system, the trajectories will converge to the agreement on both the DS side and the MG side. A comparison of voltage magnitudes and angles is presented in Fig. 23. It can be found that the results obtained by the DSPD algorithm are almost the same as those obtained by the centralized algorithm applying CPLEX solver. Fig. 24 depicts the trajectories of the relative error to the centralized solution. It is obvious that the relative error will converge to a stable value less than 10-7, which validates the accuracy of the algorithm. Additionally, the calculation time of the proposed algorithm is only about 152 s, which can fully meet the actual engineering requirements. More details of the calculation results can be found in Table 6.

Fig. 19. Illustration of self-healing feature in a real DS with two real MGs.

B. Comparison between DSPD and ADMM In order to illustrate the characteristics of the proposed method and its difference from common decentralized algorithms, the ADMM algorithm is applied in this case for comparison with the proposed algorithm. As described in [30], the convergence process of the ADMM algorithm is a process in which the primal residual and dual residual gradually decrease to zero. The primal residual reflects the gap between coupling variables, while the dual residual reflects the gap of optimization variables between two adjacent iterations. Fig. 25 shows the trajectories of the primal residual and dual residual in the iteration process of ADMM. As depicted in Fig. 25, it can be further inferred that the primal residual and dual residual will approach zero, and the ADMM calculation results will be completely consistent with the centralized calculation results when the iteration number approaches infinity. However, excessive iterations of the ADMM algorithm can consume a lot of calculation time, which is not conducive to practical engineering applications. At the same time, when the primal residual and dual residual are small enough, the calculation results of ADMM can usually reach a

Fig. 20. Trajectories of powers exchanged between DS and MGs considering MG2 plugs out.

reference voltage in the headend of the network was 10.5 kV, and the basis power was 100 MW. As depicted in Fig. 21, there were four MGs located at buses 18, 47, 78, and 103. Since there were just two sets of real MG data in our study, we copied MG1 and MG2 as MG3 and MG4. Therefore, the structure and relevant data of MG3 was the same as that of MG1, and the structure

Fig. 21. Structure of 110-bus DS with four MGs. 11

Applied Energy 252 (2019) 113361

C. Tang, et al.

Fig. 22. Trajectories of powers exchanged among DS and MGs.

Fig. 23. Comparison of voltage magnitudes and angles in centralized and decentralized methods for 110-bus system with four MGs. Table 6 Comparison of decentralized and centralized methods for 110-bus system with four MGs.

Fig. 24. Trajectories of relative error to centralized solution for 110-bus system with four MGs.

Method

Entity

Operation Cost ($)

Exchange power cost ($)

Decentralized (DSPD)

DS MG1 MG2 MG3 MG4 Total

169.8352 5.0114 2.2216 1.2846 5.7967 184.7153

From From From From

DS DS DS DS

to to to to

MG1 MG2 MG3 MG4

−78.5849 −206.5228 −52.0755 −135.2870

Centralized (CPLEX)

DS MG1 MG2 MG3 MG4 Total

169.8352 5.0114 2.2216 1.8503 5.7967 184.7153

From From From From

DS DS DS DS

to to to to

MG1 MG2 MG3 MG4

−78.5849 −206.5228 −52.0755 −135.2870

represent the voltage amplitudes of node 18 and node 111 on the DS side, and e18 and e111 represent the voltage amplitudes of node 18 and node 111 on the MG1 side. It can be found from Fig. 26 that after a period of disturbance, the trajectories of the residuals between coupling variables will stabilize at the horizontal line equal to zero. The values of the coupling variables between DS and MG1 are equal and no longer

high precision. Therefore, in practical application, a proper tolerance needs to be set as the stop criteria of ADMM iterations. Fig. 26 shows the trajectories of the residuals between some coupling variables in the process of DSPD iterations, where e18 and e111 12

Applied Energy 252 (2019) 113361

C. Tang, et al.

Table 7 Comparison of different decentralized methods. Method

Tolerance

Total Cost ($)

Relative Error to Centralized Solution

Calculation Time (s)

ADMM

1e−3 1e−4 1e−5

185.4950 184.7168 184.7151

0.0042 8.12e−6 1.08e−6

111 247 532

DSPD

/

184.7153

5.76e−8

152

termination tolerance before iterations begin. Therefore, its calculation results have better performance in terms of stability. In addition, Table 7 also shows that the DSPD algorithm can obtain nearly the same result as the centralized algorithm without consuming too much calculation time. 6. Conclusion

Fig. 25. Trajectories of primal and dual residuals in ADMM iterations.

This paper proposes a decentralized saddle-point dynamics method for a decentralized optimal power flow optimization model for a distribution system with multiple microgrids. On the premise of network structure, the distribution system and each microgrid are treated as different entities. The advantages of the proposed method can be concluded as follows: (1) Instead of the original nonlinear and non-convex power flow equations, we proposed a linearized alternating current power flow model based on the Ward equivalent method. This power flow model can implement decoupling between the distribution system and microgrids. (2) Compared with the traditional centralized models of coordination between the distribution system and microgrids, the proposed decentralized optimization model does not require global information. Each entity only exchanges the information of voltages and powers at the border buses. Then, the distribution system and each microgrid can be optimized independently. Different entities are completely confidential to each other for the information of local bus voltages and injection powers. (3) Unlike adopting common decentralized optimization methods, the proposed optimization model based on the decentralized saddlepoint dynamics approach is solved from the perspective of dynamic system control without global coordination. The convergence, accuracy, and plug-and-play features of the proposed method are well verified by the numerical results.

Fig. 26. Trajectories of residual of coupling variables in DSPD iterations.

change. At this moment, we can say the decentralized optimization converges to the saddle point. Unlike the ADMM algorithm, the termination criterion is unnecessary for the DSPD algorithm. Table 7 shows the calculation results of the ADMM algorithm and the DSPD algorithm with different tolerances. It can be found that the ADMM calculation results approach the centralized calculation results as the tolerance decreases. However, the calculation time significantly increases at the same time. In practical applications, the tolerance must be set based on experience before ADMM iterations begin, but the value of the tolerance can have a major impact on the performance of the algorithm. If the tolerance is too large, accurate calculation results are not guaranteed. If the tolerance is too small, it may consume a lot of calculation time. The DSPD algorithm does not need to set a

Acknowledgement This work was supported by the National Key Basic Research Program of China (973 Program) under Grant 2013CB228205.

Appendix A Here, we take Fig. 2 for example to obtain an expression of equivalent injection powers. MG2 is regarded as the equivalent networks first. To apply the linearized AC power flow model in the original network, the following admittance matrix should be achieved:

YII YIB 0 YLL = YBI YBB YBE , 0 YEB YEE

(A1)

where I represents the set of buses for DS and MG1, B represents BM2, and E represents the set of buses for MG2 without BM2. YII, YIB, YBI, YBB, YBE, YEB, and YEE are the corresponding entries of the admittance matrix. The inverse of YLL can be calculated based on the Gaussian elimination method, which is detailed as follows:

(YLL )

1

ZII ZIB ZIE = ZBI ZBB ZBE ZEI ZEB ZEE

(A2a) 13

Applied Energy 252 (2019) 113361

C. Tang, et al.

YBI (YII ) 1YIB

YBB = YBB

ZBI =

1Y BI

(YBB )

ZBB = (YBB )

(YII )

(A2d)

1

1

ZII = (YII )

(A2c)

1

(YBB ) 1YBE (YEE )

ZBE =

(A2b)

YBE (YEE ) 1YEB

(A2e)

1

(A2f)

1Y Z IB BI

(YII )

ZIB =

(YII ) 1YIB ZBB

ZIE =

(YII ) 1YIB ZBE

(A2g) (A2h) (A2i)

1

(YEE )

ZEB =

(YEE )

1Y Z EB BB

(A2j)

ZEI =

(YEE ) 1YEB ZBI .

(A2k)

ZEE = (YEE )

1Y Z EB BE

Next, the following formulas can be derived when shunt impedances are negligible:

YBE =

YBE (YEE ) ZBE =

(A3a)

eET YEE 1

(A3b)

eET

=

(YBB )

1Y BE

(YEE )

1

= (YBB )

1·e T E

=

(A3c)

ZBB ·eET ,

where eE denotes a column vector of E rows in which all the elements are equal to one. Since ZBB ∈ R1×1, Eq. (A3c) denotes that each element in ZBE possesses the same numerical value of ZBB. Then, we can obtain the following equations according to (A2c)–(A2k): (A4a)

RIE = RIB · eET

(A4b)

XIE = XIB · eET

(A4c)

1Y R IB BB

RIB =

(YII )

XIB =

(YII ) 1YIB XBB .

(A4d)

Considering the consistency of power flow state, the voltage should be the same before and after the Ward equivalent. Then, we can obtain the following formulas:

VLre V Lim

VLre V Lim

=

=

RLE XLE

RLB XLB

XLE RLB RLE XLB

XLB RLI RLB XLI

XLB RLI RLB XLI

XLI RLI

XLI RLI

PE QE PB . QB PI QI

(A5a)

PB QB . PI QI

(A5b)

RLE ·PE + XLE · QE + RLB ·PB + XLB ·QB = RLB · PB + XLB ·QB

(A6a)

XLE · PE

(A6b)

RLE ·QE + XLB ·PB

RLB · QB = XLB · PB

RLB ·QB .

(A5a) and (A5b) represent the deviation value exerted by the external network before and after the Ward equivalent, where L ∈ I ∪ B. (A6a) and (A6b) are the results of the simultaneous solution of (A5a) and (A5b). PB and QB denote the original active and reactive power on the border bus, PE and QE denote the original active and reactive injection power in MG1 without the border bus, and PI and QI denote the original active and reactive power in DS and MG2. P'B and Q'B represent the equivalent injection power at the border bus. According to (A4a) and (A4b), (A6a) and (A6b) can be transformed into the following equations:

RLB ·eET · PE + XLB · eET · QE + RLB ·PB + XLB ·QB = RLB · PB + XLB ·QB .

(A7a)

XLB ·eET ·PE

(A7b)

RLB · eET ·QE + XLB ·PB

RLB · QB = XLB · PB

RLB ·QB .

According to (A4c) and (A4d), we can obtain (A8a) and (A8b) as follows: T (RLB RLB

T T T T T + XLB XLB ) 1RLB RLB · eET = (RBB RBB + XBB XBB ) 1RBB RBB ·eET T T T = (RBB RBB + XBB XBB ) 1RBB RBE

(A8a)

T T T T T T (RLB RLB + XLB XLB ) 1XLB XLB ·eET = (RBB RBB + XBB XBB ) 1XBB XBB ·eET T T T = (RBB RBB + XBB XBB ) 1XBB XBE .

(A8b) 14

Applied Energy 252 (2019) 113361

C. Tang, et al.

Finally, based on (A7a)–(A7b) and (A8a)–(A8b), the detailed calculation of P'B and Q'B are represented: 2 2 PB = PB + (RBB + XBB ) 1·(RBB RBE + XBB XBE )·PE 2 2 + (RBB + XBB ) 1·(RBB XBE

(A9a)

XBB RBE )· QE

2 2 QB = QB + (RBB + XBB ) 1·(XBB RBE

RBB XBE )·PE

2 2 + (RBB + XBB ) 1·(RBB RBE + XBB XBE )·QE .

(A9b)

Furthermore, if we define F ∈ E ∪ B, (A9a) and (A9b) can be simplified as follows: 2 2 PB = (RBB + XBB ) 1·(RBB RBF + XBB XBF )·PF 2 2 + (RBB + XBB ) 1·(RBB XBF

(A10a)

XBB RBF )· QF

2 2 QB = (RBB + XBB ) 1·(XBB RBF

RBB XBF )· PF

2 2 + (RBB + XBB ) 1·(RBB RBF + XBB XBF )·QF .

(A10b)

Appendix B A common battery model usually contains the charging/discharging power limit constraint and energy limit constraint:

0

Pich

Pich , rate .

(B1a)

0

Pidis

Pidis , rate .

(B1b)

Wi,min

1

ch c Pi

Wi,0 +

Pidis

t

Wi,max ,

i

E (m).

(B1c)

d

For a single-period OPF problem, the initial energy of ESS Wi,0 in the current time can be achieved as a given parameter. Therefore, the constraint descripted in (B1c) can be rewritten as follows: ch c Pi

Wi,0 +

1

Wi,0

t

Pidis t

Wi,max ,

i

E (m).

Wi,min,

i

E (m).

(B2a) (B2b)

d

For the simultaneous Eqs. (B1a), (B2b), (B2a), and (B2b), the upper bound of charging and discharging power should be expressed as follows:

0

Pich

1

ch Pich ,max = min Pi, rate ,

t

c

0

Pidis

d

dis Pidis ,max = min Pi, rate ,

t

(Wi,max

(Wi,0

Wi,0) ,

Wi,min ) ,

i

i

E (m ),

(B3a)

E (m),

(B3b)

where and are defined as the upper bound for the charge and discharge, respectively. ch Since Wi,0 is a given parameter, Pich ,max and Pi,max are given parameters as well. Since the optimization problem without the complementary constraint is a QP problem, the optimal solution is unique. Suppose there exists an optimal pair of charge and discharge satisfying bothPich > 0 andPidis > 0. Then we can find another feasible pair satisfying P¯ich = Pich* min(Pich*. Pidis*) and P¯idis = Pidis* min(Pich*. Pidis*) , which follows the equations:

Pidis ,max

Pich ,max

P¯idis

P¯ich = Pidis*

P¯ich· P¯idis = 0,

0

P¯ich

0

P¯idis

i

<

Pich*

<

Pidis*

Since MG ¯ ch Cess (Pi ,

P¯ich

Pich*,

<

(B4b)

i

Pidis ,max ,

Pich* ,

(B4a)

E (m).

E (m).

Pich ,max ,

<

P¯idis )

i

P¯idis

i

<

MG Cess (Pich*,

(B4c)

E (m).

Pidis* ,

(B4d)

E (m). and the cost function is an increasing function for

(Pich.

Pidis) ,

we then have (B5)

Pidis*).

A lower objective value of MG can be achieved, which contradicts that the optimal pair (Pich*. Pidis*) must satisfy the complementary constraint Pich*. Pidis* = 0 .

(Pich*.

Pidis*)

is unique. Therefore, the optimal pair

Appendix C In this paper, we adopt the DSPD approach to solve the decentralized QP model in (14). The augmented Lagrange function is detailed as follows:

15

Applied Energy 252 (2019) 113361

C. Tang, et al.

L (xd , µd ,

d ,xm,

= 2 xdT Gd xd + rdT xd +

M m=1

1

+

µdT

M m=1

+

d (Aeq xd

(µ

T m

d beq )

m (Aeq xm

T d

+

m beq )+

(Aind xd

bind)

T m m (Ain xm

T d, m + µco , m (Aco xd

+

µm ,

(

1 2

d Aeq xd

binm) +

Acom xm ) +

m)

1 T x G x 2 m m m

1 2

1 2

+ rmT x m +

m beq

2 2

+

Acom xm

2 2

m Aeq xm

Acod, m xd

)

1 2

d 2 beq 2

[Aind xd 1 2

bind]+

[Ainm xm

),

2 2

binm]+

2 2

(C1)

+

where the symbol [·] represents the nonnegative projection. Then, we transform the solution of the decentralized QP model in (14) into an asymptotically stable process of the KKT conditions, which can be obtained according to the augmented Lagrange function. Thus, we construct the corresponding differential equations based on the SPD method as follows: L xd

xd =

=

d T (Gd xd + rd + (Aeq ) µd + (Aind )T d T d + (Aeq ) (Aeq xd

µd =

L d = Aeq xd µd L

d

=

xm =

d

= Aind xd

bind ,

d

> 0 or Aind xd

=

bind > 0

otherwise . m T (Gm x m + rm + (Aeq ) µm + (Ainm )T

=

L m = Aeq xm µm L

m

d beq ) + (Aind )T [Aind xd

m

binm , 0

m

xd

m T ) µco, m

bind]+ ).

(C2)

.

m

m beq ) + (Ainm )T [Ainm x m

(C4) m T (Aco ) µco

binm]+ )

(C5)

m beq

= Ainm xm

µco, m = Acod

(Acod

(C3)

0 L xm

M m=1

+

d beq .

m T m + (Aeq ) (Aeq xm

µm =

d

(C6) m

> 0 or Ainm xm

binm > 0 (C7)

otherwise .

Acom xm

(C8)

References [17]

[1] Dall'Anese E, Hao Z, Giannakis GB. Distributed optimal power flow for smart microgrids. IEEE Trans Smart Grid 2013;4:1464–75. [2] Wang Z, Chen B, Wang J, Kim J, Begovic MM. Robust optimization based optimal DG placement in microgrids. IEEE Trans Smart Grid 2014;5:2173–82. [3] Wang Z, Chen B, Wang J, Begovic MM, Chen C. Coordinated energy management of networked microgrids in distribution systems. IEEE Trans Smart Grid 2015;6:45–53. [4] Hu X, Liu T. Co-optimisation for distribution networks with multi-microgrids based on a two-stage optimisation model with dynamic electricity pricing. IET Gener Transm Distrib 2017;11:2251–9. [5] Fathi M, Bevrani H. Adaptive energy consumption scheduling for connected microgrids under demand uncertainty. IEEE Trans Power Delivery 2013;28:1576–83. [6] Fang X, Yang Q, Wang J, Yan W. Coordinated dispatch in multiple cooperative autonomous islanded microgrids. Appl Energy 2016;162:40–8. [7] Lv T, Ai Q. Interactive energy management of networked microgrids-based active distribution system considering large-scale integration of renewable energy resources. Appl Energy 2016;163:408–22. [8] Zadsar M, Haghifam MR, Miri Larimi SM. Approach for self-healing resilient operation of active distribution network with microgrid. IET Gener Transm Distrib 2017;11:4633–43. [9] Du Y, Wang Z, Liu G, Chen X, Yuan H, Wei Y, et al. A cooperative game approach for coordinating multi-microgrid operation within distribution systems. Appl Energy 2018;222:383–95. [10] Jahangiri P, Aliprantis DC. Distributed Volt/VAr control by PV inverters. IEEE Trans Power Syst 2013;28:3429–39. [11] Farivar M, Neal R, Clarke C, Low S. Optimal inverter VAR control in distribution systems with high PV penetration. In: 2012 IEEE power and energy society general meeting; 2012. p. 1–7. [12] Ghaddar B, Marecek J, Mevissen M. Optimal power flow as a polynomial optimization problem. IEEE Trans Power Syst 2016;31:539–46. [13] Baran ME, Wu FF. Optimal capacitor placement on radial distribution systems. IEEE Trans Power Delivery 1989;4:725–34. [14] Farivar M, Low SH. Branch flow model: relaxations and convexification—Part I. IEEE Trans Power Syst 2013;28:2554–64. [15] Lavaei J, Low SH. Zero duality gap in optimal power flow problem. IEEE Trans Power Syst 2012;27:92–107. [16] Molzahn DK, Dorfler F, Sandberg H, Low SH, Chakrabarti S, Baldick R, et al. A

[18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31]

16

survey of distributed optimization and control algorithms for electric power systems. IEEE Trans Smart Grid 2017;8:2941–62. Guggilam SS, Dall'Anese E, Chen YC, Dhople SV, Giannakis GB. Scalable optimization methods for distribution networks with high PV integration. IEEE Trans Smart Grid 2016;7:2061–70. Xydas E, Marmaras C, Cipcigan LM. A multi-agent based scheduling algorithm for adaptive electric vehicles charging. Appl Energy 2016;177:354–65. Huang J, Li Z, Wu QH. Coordinated dispatch of electric power and district heating networks: a decentralized solution using optimality condition decomposition. Appl Energy 2017;206:1508–22. Kofinas P, Dounis AI, Vouros GA. Fuzzy Q-learning for multi-agent decentralized energy management in microgrids. Appl Energy 2018;219:53–67. Liu J, Cheng H, Zeng P, Yao L, Shang C, Tian Y. Decentralized stochastic optimization based planning of integrated transmission and distribution networks with distributed generation penetration. Appl Energy 2018;220:800–13. Wang X, Wang C, Xu T, Meng H, Li P, Yu L. Distributed voltage control for active distribution networks based on distribution phasor measurement units. Appl Energy 2018;229:804–13. Wang X, Wang C, Xu T, Guo L, Li P, Yu L, et al. Optimal voltage regulation for distribution networks with multi-microgrids. Appl Energy 2018;210:1027–36. Xie M, Ji X, Hu X, Cheng P, Du Y, Liu M. Autonomous optimized economic dispatch of active distribution system with multi-microgrids. Energy. 2018;153:479–89. Gao H, Liu J, Wang L, Wei Z. Decentralized energy management for networked microgrids in future distribution systems. IEEE Trans Power Syst 2018;33:3599–610. Zheng W, Wu W, Zhang B, Sun H, Liu Y. A fully distributed reactive power optimization and control method for active distribution networks. IEEE Trans Smart Grid 2015:1-. Peng Q, Low SH. Distributed optimal power flow algorithm for radial networks, I: Balanced single phase case. IEEE Trans Smart Grid 2018;9:111–21. Kargarian Marvasti A, Yong F, DorMohammadi S, Rais-Rohani M. Optimal operation of active distribution grids: a system of systems framework. IEEE Trans Smart Grid 2014;5:1228–37. Hur D, Park JK, Kim BH. Evaluation of convergence rate in the auxiliary problem principle for distributed optimal power flow. IEE Proc – Gen Transm Distrib 2002;149:525–32. Distributed Boyd S. Optimization and statistical learning via the alternating direction method of multipliers. Found Trends® Mach Learn 2010;3:1–122. Feijer D, Paganini F. Stability of primal–dual gradient dynamics and applications to

Applied Energy 252 (2019) 113361

C. Tang, et al. network optimization. Automatica 2010;46:1974–81. [32] Ma X, Elia N. Convergence analysis for the primal-dual gradient dynamics associated with optimal power flow problems. 2015 European control conference (ECC). 2015. p. 1261–6. [33] Richert D, Cortes J. Robust distributed linear programming. IEEE Trans Autom Control 2015;60:2567–82. [34] Lo KL, Peng LJ, Macqueen JF, Ekwue AO, Cheng DTY. An extended ward equivalent approach for power system security assessment. Electr Power Syst Res 1997;42:181–8.

[35] Wang Z, Wu W, Zhang B. A fully distributed power dispatch method for fast frequency recovery and minimal generation cost in autonomous microgrids. IEEE Trans Smart Grid 2016;7:19–31. [36] Hu B, Wang H, Yao S. Optimal economic operation of isolated community microgrid incorporating temperature controlling devices. Prot Contr Mod Power Syst 2017;2:7. [37] Gao H, Liu J, Wang L. Robust coordinated optimization of active and reactive power in active distribution systems. IEEE Trans Smart Grid 2018;9:4436–47.

17