Sensitivity analysis of volt-VAR optimization to data changes in distribution networks with distributed energy resources

Applied Energy 261 (2020) 114331

Contents lists available at ScienceDirect

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

Sensitivity analysis of volt-VAR optimization to data changes in distribution networks with distributed energy resources Davye Mak, Daranith Choeum, Dae-Hyun Choi

T

⁎

School of Electrical and Electronics Engineering, Chung-ang University, Dongjak-gu, Seoul 156-756, Republic of Korea

HIGHLIGHTS

framework to quantify the sensitivity of VVO to data change in active distribution networks. • Novel of the closed-form linear sensitivity matrix. • Derivation of the sensitivity matrix using various types of VVO data. • Validation • Fast and accurate analysis of data quality in VVO. ARTICLE INFO

ABSTRACT

Keywords: Volt-VAR optimization Sensitivity analysis Data change Karush-Kuhn-Tucker conditions Active distribution network

The Volt-VAR optimization (VVO) model is combined with distributed energy resources (e.g., solar photovoltaic (PV) systems and energy storage systems (ESSs)), advanced sensors and communication networks, as well as demand response programs for reliable and economical distribution grid operations. Unexpected changes in the data used for VVO can degrade the VVO performance. This study examines the impact of changes in VVO data on the optimal solutions of VVO (e.g., nodal voltage, real/reactive power flow at the substation, and total deviation of voltages from the nominal voltage). Using the perturbed Karush-Kuhn-Tucker conditions of the VVO formulation, we develop a linear matrix that evaluates the sensitivity of VVO with respect to change in various types of data—distribution line parameters, the predicted PV generation outputs, demand reduction, and load exponents in exponential load models. The proposed sensitivity matrix can be used by system operators as an analysis tool to quickly identify the data that influences the optimal VVO solution the most by accurately measuring and prioritizing sensitivities. The results of a simulation study using the developed sensitivity matrix are illustrated and verified in two distribution test systems (IEEE 33-node and 123-node systems) with an onload tap changer (OLTC), capacitor banks (CBs), PV systems, and ESSs.

1. Introduction The Volt-VAR optimization (VVO) is one of the crucial functions in advanced distribution management systems (ADMSs), employed for efficient managing and controlling of active distribution networks with distributed energy resources (DERs) (e.g., solar photovoltaic (PV) system, energy storage system (ESS), and electric vehicle), advanced metering infrastructure (AMI), and demand response programs [1]. In active distribution networks, VVO computes the optimal voltage magnitude for any node in an optimization problem. Thus, the total voltage deviation from the nominal voltage is minimized while considering all operational constraints of distribution systems. Operational constraints in the VVO are formulated using many different types of data—voltage limits, distribution line parameters, operation parameters of PV/ESS, ⁎

the predicted PV generation output and load demand, demand reduction and load exponents in exponential load models. A change in these heterogeneous data may result in a severe malfunction of the VVO. This study attempts to provide an analytical framework for the quantification of the impact of changes in VVO data on the performance of VVO. The primary goal of the VVO is to determine optimal voltage levels along the distribution feeder under all loading conditions through the coordination of conventional voltage regulators such as on-load tap changers (OLTCs) at substations, step-type voltage regulators (SVRs), and capacitor banks (CBs) [2]. Recently, advanced information and communication technology (ICT) such as AMI with smart meters provides voltage regulators with real-time loading and voltage measurements, thereby adjusting the voltage profile more accurately [3]. Furthermore, the smart inverters of DERs can be employed as voltage

Corresponding author. E-mail address: [email protected] (D.-H. Choi).

https://doi.org/10.1016/j.apenergy.2019.114331 Received 22 August 2019; Received in revised form 12 November 2019; Accepted 3 December 2019 0306-2619/ © 2019 Elsevier Ltd. All rights reserved.

Applied Energy 261 (2020) 114331

D. Mak, et al.

• To the best of our knowledge, data sensitivity analysis of mixed

regulators, through the injection or absorption of their real and reactive power, to maintain the desired local voltage profile based on the designed droop control curve when voltage violations occur rapidly owing to intermittent PV generation [4,5]. A mixed-integer second-order cone programming (MISOCP) model was formulated to determine the optimal parameters for the droop control curves of DERs [6]. More recently, with an increasing number of EVs connected to distribution networks, a new multi-agent based VVO method was developed where fast vehicle-to-grid reactive power dispatch, charging coordination and VVO for EVs are carried out in a distributed and fast manner [7]. The impact of different EV penetration levels on an AMI-based quasi realtime VVO was investigated [8]. As the coupling between the VVO and the ICT system with DERs is further strengthened for efficient voltage regulation of an entire active distribution network, a large number of heterogeneous data including sensor measurements and operating parameters for networks and DERs are used for the VVO process as input data. With the aforementioned circumstances, the change in these data due to a natural error or manmade attack may significantly distort the optimal solution of VVO. Recent studies have shown that an adversary can maliciously change the distribution feeder voltage profile by stealthily malfunctioning the voltage regulators through the injection of false data into sensors [9–11]. This study investigates the sensitivity of the VVO with respect to changes in various types of data that are used for the VVO process. We specifically aim to evaluate the impact of these data on the VVO optimal solutions (e.g, nodal voltage and real/reactive power flow at the substation) and the optimal total voltage variation. A large body of literature has been accumulated on the subject of sensitivity analysis regarding data changes for power system applications at the transmission level: (i) locational marginal price (LMP) in the DC and AC optimal power flow models [12,13], (ii) state estimation [14,15], and iii) lookahead security constrained economic dispatch (SCED) [16]. Recently, we conducted data perturbation-based sensitivity analysis of power system applications at the distribution level, including the home energy management system (HEMS) [17] and distribution system state estimation (DSSE) integrated with HEMS [18]. Moreover, voltage sensitivity analysis was performed in distribution systems, including development of the expression for voltage sensitivity that illustrates the relationship between real power and voltage magnitude in a DC distribution system [19], calculation of voltage sensitivity using historical smart meter data without low voltage (LV) grid topology information [20] and of decentralized voltage sensitivity with local and neighborhood measurements only [21], and derivation of the nodal sensitivity of voltage and real/reactive power at nodes with distribution phasor measurement unit (DPMU) [22]. Although much work has been conducted on the subject of sensitivity analysis for power system applications at both the transmission and distribution system levels, to the best of our knowledge, no analytical sensitivity framework has been developed to directly assess the impact of various types of data on the VVO. Previous studies [19–22] have focused on a sensitivity analysis of nodal voltage magnitude considering the change in a few variables such as real and reactive power without explicitly considering the VVO process. Our study is motivated by a desire to investigate the impact of change in all input data applied in VVO on the optimal solution of VVO. In this study, we assume that data change refers to both natural noise and man-made attacks. Our proposed sensitivity framework is based on a perturbation approach using Karush-Kuhn-Tucker (KKT) conditions for general nonlinear optimization problems [23], which have been verified and illustrated in [13,14]. The main contributions of this study are summarized as follows:

•

•

integer nonlinear programming (MINLP)-based VVO via a perturbation approach based on its corresponding KKT equations has not been studied yet; therefore the proposed sensitivity analysis framework will provide system operators with an online analysis tool to quickly assess the impact of data on the performance of the VVO. The key part of the proposed approach is to develop a closed-form linear sensitivity matrix based on the perturbed KKT equations from the VVO formulation in the following three steps: (i) VVO formulation using MINLP, (ii) reformulation of the MINLP-based VVO into nonlinear programming (NLP)-based VVO with an integer optimal solution from the MINLP-based VVO, and (iii) perturbation of KKT equations from NLP-based VVO. Through the performance evaluation of the developed sensitivity matrix in the IEEE 33-node and 123-node systems, we show that the proposed approach can successfully generate the accurate sensitivities within the allowable VVO execution time.

The remainder of this paper is organized as follows. Section 2 introduces a system model and mathematical optimization problem for VVO. In Section 3, we develop a mathematical framework based on the perturbed KKT conditions of the VVO formulation to quantify the sensitivity of VVO to change in data. Section 4 presents numerical examples that illustrate the impact of various types of data on VVO in the IEEE 33-node and 123-node test feeders, along with validation of the accuracy and computation time of the proposed sensitivity framework. We make concluding remarks and suggest future work in Section 5. 2. Preliminary The main notations used throughout this paper are summarized in Table 1. Bold symbols represent vectors or matrices. Hat symbols represent predictions of true parameter values. The other undefined symbols in the nomenclature section are explained in the text. 2.1. System model We consider a radial medium-voltage (MV) distribution system with an OLTC-fitted transformer as shown in Fig. 1 and some nodes in the system are equipped with CBs and PVs/ESSs. Let and be the sets of nodes and lines of the system, respectively. Then, let d be the ,f d be the index of any node that is index of any node and f adjacent to node d. These two nodes form a line df . Using linearized DistFlow equations [24,25], the real/reactive power flow in line df line (Ptline , df , Qt , df ) and voltage magnitude (Vt , d ) for node d at time t are expressed as, respectively,

Ptline , df =

node Ptline , fk + Pt , f , fk k

f

k

f

Qtline , df =

node Qtline , fk + Qt , f , fk

Vt , f = Vt , d

line rdf Ptline , df + x df Qt , df

V0

, fk

, fk

, d, f

df

(1)

df

,d

(2)

f , df

.

(3)

In the right hand side of (1) and (2), the first term represents the line sum of real power flow Ptline , fk and reactive power flow Qt , fk from node f to k where node k belongs to a set of nodes f that includes all nodes connected downstream to node f. Eq. (3) expresses voltage drop between node d and f, and V0 in here is the nominal voltage magnitude (1.0 p.u.) at the substation. In the real and reactive power flow equations above, the nodal balance equations of the real (Ptnode , d ) and reactive power (Qtnode , d ) for node d at time t can be written in terms of the real and/or reactive power of the load, CBs, PVs, and ESSs

• We

present a closed-form analytical framework to evaluate the sensitivity of the optimal VVO solution with respect to changes in various types of data used for VVO.

Ptnode = Ptload ,d ,d 2

PV

Pt , d

PtESS,dch + PtESS,ch ,d ,d

(4)

Applied Energy 261 (2020) 114331

D. Mak, et al.

Table 1 Notation. Resistance of the line df

rdf

Reactance of the line df

x df

Real power flow from node d to f at period t

Ptline , df

Qtline , df

Reactive power flow from node d to f at period t

Ptnode ,d

Net real power consumption for node d at period t Net reactive power consumption for node d at period t

Qtnode ,d

Voltage magnitude for node d at period t Minimum(Maximum) limit of the allowed voltage range for node d Nominal voltage magnitude Predicted PV real power output for node d at period t

Vt, d V min(max)

V nom PV

Pt , d

QtPV ,d

PV reactive power output for node d at period t

QdPV,min(max) PtESS,ch(dch) ,d QtESS ,d

Charging(Discharging) real power of ESS for node d at period t

Minimum(Maximum) limit of reactive power of PV for node d

Fig. 1. Typical radial distribution system.

Charging or discharging reactive power of ESS for node d at period t Minimum(Maximum) limit of reactive power of ESS for node d

QdESS,min(max)

subsection. 2.2. Volt-VAR optimization model

Maximum charging real power of ESS for node d

PdESS,ch,max

Minimum charging real power of ESS for node d

PdESS,ch,min

In this study, we consider the scenario where the VVO schedules the optimal voltage profile along the distribution feeder through the coordination of OLTC, CBs, and the smart inverters of PV systems/ESSs as shown in Fig. 1. We assume that VVO obtains all input data prior to the optimal voltage scheduling, including: (i) line resistance and reactance, (ii) predicted load demand and PV generation output, (iii) operational parameters of DERs and voltage regulators, and (iv) parameters for load exponential models. with a scheduling period t and prediction For each node d horizon Nt , the VVO problem is generally formulated as the MINLP model:

Maximum discharging real power of ESS for node d

PdESS,dch,max

Minimum discharging real power of ESS for node d

PdESS,dch,min

Maximum energy capacity of ESS for node d

EdESS,max

State of the charge of ESS for node d at period t Maximum(Minimum) state of the charge of ESS for node d

SOCt, d SOCdmax(min)

Charging(Discharging) efficiency of ESS for node d

ch(dch) d

Reactive power output of CB for node d at period t

QtCAP ,d

Real load consumption for node d at period t

Ptload ,d

Reactive load consumption for node d at period t

Qtload ,d

Predicted real load consumption at nominal voltage for node d at period t Predicted reactive load consumption at nominal voltage for node d at period t Load reduction ratio in [0, 1] by demand response for node d at period t Exponent of real load consumption for node d in exponential load model Exponent of reactive load consumption for node d in exponential load model Integer tap position of OLTC at period t

load,nom Pt , d load,nom

Qt, d

XtDR ,d d d

TaptOLTC

Nt

NCapmax

aOLTC btCAP ,d

Nt Nd

Qtnode = Qtload ,d ,d

QtPV ,d

QtESS ,d

QtCAP ,d ,

: Ptline , df =

load,nom

Ptload , d = Pt , d

(Vt , d )

load,nom

Qtload , d = Qt , d

d

(Vt , d) d.

(Qtload ,d )

node Ptline , fk + Pt , f , fk k

f

k

f

: Qtline , df =

, fk

node Qtline , fk + Qt , f , fk

df

, fk

line rdf Ptline , df + x df Qt , df

: Vt , f = Vt , d

, d, f

V0

(9)

df

(10)

,d

f , df

: Vt ,1 = V nom + aOLTCTaptOLTC

: Ptnode = Ptload ,d ,d

PV

QtPV ,d

load,nom

: Ptload , d = Pt , d

are ex-

XtDR , d )(Vt , d )

(1

(13)

QtCAP ,d

XtDR , d )(Vt , d )

(1

load,nom

: Qtload , d = Qt , d

QtESS ,d

(14) (15)

d

(16)

d

CAP CAP,nom : QtCAP , d = bt , d Qd

(6) (7)

: SOCt , d = SOCt

Here, the real and reactive load consumptions at each node d depends on the corresponding voltage and certain specified exponents, d and d , respectively. Based on the system model illustrated above, the detailed VVO formulation along with the operational constraints of PV/ESS, OLTC, and CB along with the demand reduction is introduced in the following

: SOCdmin

SOCt , d

: PdESS,ch,min btESS ,d : PdESS,dch,min (1 3

1, d

(17) ch ESS,ch d Pt , d EdESS,max

+

PtESS,dch ,d dch ESS,max d Ed

(18) (19)

SOCdmax

PtESS,ch ,d btESS ,d )

(11) (12)

PtESS,dch + PtESS,ch ,d ,d

Pt , d

µ : Qtnode = Qtload ,d ,d

(5)

where the real and reactive load consumption pressed using the following load exponential model [26,27].

(Ptload ,d )

(8)

s.t.

Binary switch status of the capacitor for node d at period t; “1” for ON and “0” OFF Binary charging state of ESS for node d at period t.; “1” for charging and “0” otherwise. Total number of prediction horizon steps Total number of nodes

btESS ,d

V nom )2

(Vt , d t=1 d=1

Maximum switching operations of OLTC during the prediction horizon Nt Maximum switching operations of CB during the prediction horizon Nt Step size of change in OLTC tap position

NTapOLTC,max

Nd

minJ =

PdESS,ch,max btESS ,d PtESS,dch ,d

PdESS,dch,max (1

(20)

btESS ,d )

(21)

Applied Energy 261 (2020) 114331

D. Mak, et al. Nt

TaptOLTC

TaptOLTC 1

OLTC NTapmax

t=1

min J x, a

(22)

(27)

x

s.t.

Nt

btCAP ,d

btCAP 1, d

NCapmax

t=1

: QdPV,min

QtPV ,d

QdPV,max

(24)

: QdESS,min

QtESS ,d

QdESS,max

(25)

: V min

Vt , d

V max.

: f (x, a) = 0,

(23)

: g (x, a)

(28)

0,

where J (x, a) is a scalar objective function, and f (x, a) and g (x, a) depict the equality and inequality constraints with the corresponding Lagrangian multipliers and , respectively. x is the decision variable vector and a is the data vector for the optimization problem. 3.1.1. KKT conditions Given the Lagrangian function ( (x, a) = J (x, a) + + g (x, a) ), the KKT conditions are written as follows:

(26)

In this formulation, the objective function is to minimize the total deviation of voltages from the nominal voltage for all nodes during the prediction horizon in (8). Eqs. (9)–(11) represent the linearized distribution real power flow, reactive power flow, and voltage drop for line df at the scheduling period t, respectively. Eq. (12) represents substation voltage, which can be determined by the OLTC tap position TaptOLTC , along with the step size for changing tap positions aOLTC . Eqs. (13) and (14) express nodal real and reactive power balance. In (15) and (16) which define the real and reactive powers of load consumption based on an exponential load model, XtDR , d indicates the amount of load reduction within [0, 1] obtained via the demand response program. Eq. (17) represents the reactive output that is supported by the CB for node d at scheduling period t. Here, QdCAP,nom is the capacitors’ size, and btCAP ,d is a binary decision variable that determines the switch status of these capacitors. For the ESS at node d, (18) defines the operational dynamics of its SOC at time t in terms of its SOC at a previous time t 1, its battery capacity EdESS,max , its charging and discharging efficiency, dch and ddch , and its charging and discharging power, PtESS,ch and PtESS,dch . ,d ,d The capacity constraint of the SOC is presented in (19). Eqs. (20) and (21) depict constraints on the charging and discharging powers of the ESS, where btESS , d represents a binary decision variable that determines the on/off status of the ESS. During the prediction horizon Nt , the total number of switching operations for the OLTC and CBs is limited by their OLTC corresponding switching thresholds, NTapmax and NCapmax in (22) and (23). Eqs. (24) and (25) represents the reactive power capability of the PV system and the ESS at node d. The range of allowable voltages for all nodes can be expressed in (26) where V min and V max are selected to be 0.95 p.u. and 1.05 p.u., respectively. The decision variable vectors from the MINLP-based VVO problem V , Pline , Qline , Pnode , Qnode , P load , Qload , QCAP , are denoted as ESS,ch ESS,dch PV SOC, P ,P , Q , QESS , TapOLTC , and bCAP . The other variable , , ) corresponding vectors ( , , , , , µ , , , , , , , , to the equality/inequality constraints are the Lagrangian multipliers with non-negative values. Indeed, no Lagrangian multipliers exist in the MINLP optimization problem. In this paper, the Lagrangian multiplier vectors above are used to derive the perturbed KKT conditions of the NLP-based VVO problem that is converted from the MINLP-based VVO problem, which is illustrated in the next section.

K1 x (x, a) = 0 K2 f (x, a) = 0 , K3 g (x, a) 0 , K4 g (x, a) =0,

xJ

(x, a) +

x f (x ,

a) +

x g (x ,

f (x , a)

a) =0 ,

0,

where (K1) depicts the first-order optimality conditions, (K2) and (K3) denote the primal feasibility conditions, and (K4) depicts either the dual feasibility or the complimentary slackness conditions. In this study, only the binding inequality constraints ( > 0 ) are considered for our sensitivity analysis. 3.1.2. Perturbed KKT conditions The perturbed equations of the previously described KKT conditions subject to x , a, , and , are written as follows: P1 [ x J (x, a)]T d x + [ a J (x, a)]T da -dJ = 0 P2 [ xx J (x, a) + xx f (x , a) + xx g (x , a) ]dx +[ xa J (x, a) + xa f (x , a) + xa g (x , a) ]da + x (f (x, a)) d + x (g (x , a)) d = 0 , P3 [ x f (x, a)]T d x + [ a f (x, a)]T d a = 0 , P4 [ x g (x , a)]T d x + [ a g (x, a)]T d a = 0 . Employing (P1) (P4), the following linear matrix equation is constructed, (29)

S1 [d xd d dJ ]T = S 2 d a

where S1 and S2 are matrices whose elements are the coefficients of d x, d a, d , d , and dJ in the perturbed KKT conditions (detailed expressions of S1 and S2 are provided in [23]). Finally, the following matrix S = S1 1 S 2 offers the sensitivity of the decision variables, the Lagrangian multipliers, and the objective function with respect to the data changes.

S = S1 1 S 2 =

d x d d dJ da da da da

T

.

(30)

3.2. Proposed sensitivity matrix

3. Proposed sensitivity framework for VVO

The main goal of this study is to derive the following two linear sensitivity matrices of the optimal solutions x and objective function J for VVO with respect to the change in different types of data a :

This section introduces the proposed analysis framework used to quantify the sensitivity of the VVO with respect to changes in the data. In Section 3.1, the sensitivity framework for a general NLP optimization problem is presented, which can be obtained by perturbing the decision variables and data in the KKT equations from the NLP optimization problem. In Section 3.2, the proposed sensitivity matrix for the NLPbased VVO is proposed along with the classification of both the optimal solution and various types of data for VVO.

S(x) =

dx dJ , S (J ) = . da da

(31)

For t = 1, …, Nt and d = 1, …, Nd , the optimal solution vector x , the data vector a , and the Lagrangian multiplier vector lp are written as, respectively, line node node load load PV ESS CAP ESS,ch x = [Vt , d, Ptline , , df , Qt , df , Pt , d , Qt , d , Pt , d , Qt , d , Qt , d , Qt , d , Qt , d , Pt , d

3.1. Sensitivity analysis in an optimization problem

PtESS,dch , SOCt , d ]T ,d (32)

A general NLP optimization problem is formulated as follows: 4

Applied Energy 261 (2020) 114331

D. Mak, et al.

Step 1 The optimal continuous and discrete decision variable vector xMINLP from the original MINLP-based VVO problem is calculated. Step 2 The MINLP-based VVO problem is transformed into the NLPbased VVO problem by assuming the discrete decision variables obtained from Step 1) as the fixed data. However, all the constraints that include only discrete decision variables are removed (e.g, (22), (23)). Subsequently, the optimal continuous decision variable vector xNLP for the NLP optimization problem is calculated. Step 3 The KKT equations of the NLP optimization problem in Step 2) are derived with xNLP . All the KKT equations together with the objective function J are perturbed with respect to the decision variables (32), Lagrange multipliers (34), and data (33) while the KKT conditions still hold. Finally, the desired sensitivity matrix is constructed with the perturbed KKT equations.

Table 2 Classification of data in VVO. Group (G1)

Description

Data

Line parameter

(G2)

Predicted PV generation output/load

(G3)

Load model

(G4)

Operational parameter of PV/ESS

load,nom PV load,nom a = [rdf , x df , P^t , d , P^t , d , Q^t , d , XtDR ,d ,

rdf , x df PV

load,nom

Pd , Pd

load,nom

, Qd

DR d , d , Xd ch(dch) ESS,ch(dch),max , Pd d QdPV,max , QdESS,max , EdESS,max

ch dch PV,min , d , d , d , Qd ESS,max min max PV , max ESS, min min max Qd , Qd , Qd , SOCd , SOCd , V , V , PdESS, ch, min, PdESS, ch, max , PdESS,dch,min, PdESS,dch,max , EdESS,min, EdESS,max ]T d,

(33)

lp = [ , ,

, , , µ, , , , ,

,

, ,

,

,

]T .

Due to space limitations, all KKT equations and their corresponding perturbed equations are not shown in this paper. Instead, we provide a simple example, which aids our proposed approach to be understood clearly. For example, we consider the situation where the MINLP-based VVO determines the charging status of the ESS with btESS , d = 1 (Step 1). Using btESS , d = 1, the capacity constraint (20) in the MINLP problem is PdESS,ch,max in the NLP problem. With modified to PdESS,ch,min PtESS,ch ,d this modified constraint, the NLP problem recalculates the optimal that is binding at PdESS,ch,max (Step 2). The binding charging power PtESS,ch ,d ESS,ch ESS,ch,max = Pd constraint Pt , d is one of the KKT equations in the NLP problem, and the perturbed equation of this constraint in regards to the decision variable PtESS,ch and data PdESS,ch,max is expressed as ,d ESS,ch ESS,ch,max dPt , d = dPd (Step 3). Using the aforementioned approach, all perturbed KKT equations are finally written in a linear matrix equation form (35).

(34)

Here, the data vector a in (33) is categorized into four groups ((G1) (G4)) as shown in Table 2. The first group includes the line resistance and reactance parameters. In the second group, the data represent the predicted values of the PV generation real output and the real/reactive load consumption. The third group includes the parameters of demand reduction and exponent values in exponential load model. The fourth group corresponds to the operation parameters of PV and ESS, including the capacity limit of PV reactive power, the charging/discharging efficiency and capacity limit of real and/or reactive power along with the battery capacity limit for ESS. To develop the desired linear sensitivity matrix, the MINLP-based VVO problem in Section 2.2 is first converted into the NLP optimization problem with the fixed discrete decision variables from the MINLP optimization problem. The KKT equations of the NLP optimization problem are subsequently derived assuming that a small change at the right sides of the constraints in the MINLP optimization problem results in no change of a discrete variable [28]. Finally, employing the perturbation approach illustrated in Section 3.1, the desired sensitivity matrix is obtained by simultaneously perturbing both the continuous decision variables and data in the KKT equations. Fig. 2 illustrates the procedure of the proposed sensitivity approach with the classification of VVO input data. As shown in Fig. 2, we note that the proposed sensitivity analysis approach works when an optimal VVO solution is known. The detailed procedure for the derivation of the sensitivity matrix is summarized as follows:

M 1 2 3

0 T 2

0 0

0 T 3

0 0

1 0 0 0

dx d lp = dJ

1 2 3

[d a ].

4

(35)

Taking the inverse of the left-hand side matrix in (35) on both sides of (35), we obtain the sensitivity matrices that evaluates the impact of data on the optimal VVO solutions: S(x) = d x/ d a and S(J ) = dJ / d a . The derivations of the submatrices in the linear matrix Eq. (35) are provided in Appendix A.

Fig. 2. Proposed sensitivity analysis framework for VVO along with data classification. 5

Applied Energy 261 (2020) 114331

D. Mak, et al.

Fig. 3. Modified IEEE distribution test feeders with OLTC, CBs, PV systems and ESSs: (a) 33-node system, and (b) 123-node system.

4. Numerical results

Qtline ,1 XdDR

In this section, the developed sensitivity matrices and in Section 3 are tested and validated to assess the impact of changes in various types of data on VVO. The impact assessment is illustrated in the modified IEEE 33-node [25] and 123-node test system [29] as shown in Fig. 3. A simulation environment for the sensitivity analysis is described in Section 4.1, which is followed by the subsequent five subsections that include the sensitivity results for different data groups in Table 2:

S(x)

S (J )

Vt , d XdDR Ptline ,1 rdf Ptline ,1 XdDR

Qtline ,1 rdf

,

Vt , d

Vt , d

,

d

,

,

(x ) j

(J ) j

,

Qtline ,1 xdf

Ptline ,1

PV

Pt , d

d

,

J

(puV) .

d

x / ai

j

ai

J / ai

j

(36)

(37)

(J ) j

As shown in Fig. 3(a), the modified IEEE 33-node test feeder is equipped with an OLTC, nine CBs, and four smart inverters for both PV systems and ESSs. An integer tap position of OLTC at the substation ranges from −10 to 10 with its step change aOLTC = 0.005 p.u. The maximum number of tap changes for OLTC during the predicted horOLTC = 3. The maximum output of each CB is izon Nt = 1 is set to NTapmax CAP,nom Qd = 30 kVAr . Each ESS and PV system have a size of 30 kVA and 100 kVA, respectively, and a battery capacity of ESS is EdESS,max = 100 kWh . For the ESSs, the maximum charging and discharging power is 27 kW, and the minimum charging and discharging powers is 0 kW. For the simulation study in the IEEE 123-node test feeder as shown Fig. 3(b), the three-phase unbalanced feeder is assumed to be the three-phase balanced feeder that is equipped with an OLTC at the node 0 and three voltage regulators at node 17, 25, and 73, respectively. Each of nodes 7, 33, and 79 is equipped with ESS with a size of 50 kVA with 200 kWh battery capacity. For the ESSs, the maximum charging and discharging power is 45 kW, and the minimum charging and discharging power is 0 kW. Nodes 7, 33, and 79 are also equipped with PV system with a size of 345 kVA, 345 kVA, and 690 kVA, respectively [30]. Four CBs with a size of 200 kVAr, 50 kVAr, 50

(puMW)

PV P^t , d

J

4.1. Simulation setup

(nounit)

Qtline ,1

J , XdDR

where and represent the average of the absolute normalized sensitivities of the decision variable and the objective function with respect to any data ai in the set j , respectively. Here, j is the set for type j data. For example, rd depicts the set for resistance at any line d. The cardinality of the set j implies the number of elements in j . Finally, the accuracy and computation time of the proposed sensitivity framework is validated in Section 4.6.

d

(puMVAr/puOhm),

ai j

= (x ) j

puV/puMW

Ptline ,1

(puV) 2/puOhm ,

j

=

i=1

(puV)

puMW/puOhm ,

xdf

d

,

PV

(puMVAr).

d

i=1

d

Ptline ,1

Ptline ,1

Vt , d

Qtline ,1

To fairly compare the sensitivities with different units, the following normalized sensitivity metrics are proposed, and the corresponding sensitivity results are provided in Section 4.5:

In this simulation study, for better illustrating the results of the proposed approach, Sections 4.2, 4.3 and 4.4 provide the sensitivity results in the IEEE 33-node system, whereas Sections 4.5 and 4.6 provide the sensitivity results and validate the accuracy and computation time of the proposed framework in both the IEEE 33-node and 123-node systems. In Section 4.2, the following sensitivities for voltage magnitude (Vt , d ) for any node d at any period t and the injected real/reactive power line flow (Ptline ,1 , Qt ,1 ) for node 1 at any period t are computed:

Pt , d

d

J J , rdf xdf

proposed sensitivity framework

puV/puOhm ,

,

Similarly, the following sensitivities for the total voltage deviation and real power loss J are presented in Sections 4.3 and 4.4:

• Section 4.2: Sensitivity analysis of the decision variable values to (G1) (G3) data change • Section 4.3: Sensitivity analysis of the objective function value (total voltage deviation) to (G1) and (G3) data change • Section 4.4: Sensitivity analysis of the objective function value (total real power loss) to (G1) and (G3) data change • Section 4.5: Normalized sensitivity analysis of the decision variable and objective function values to (G1) (G4) data change • Section 4.6: Validation of the accuracy and computation time of the

Vt , d Vt , d , rdf xdf

Qtline ,1

,

(puMVAr/puMW)

6

Applied Energy 261 (2020) 114331

D. Mak, et al.

Fig. 4. Schedule of the optimal voltage magnitude for VVO at 9 A.M. in a modified IEEE 33-node test feeder.

kVAr, and 50 kVAr are connected to nodes 111, 116, 118, and 120, respectively. For both the IEEE 33-node and 123-node test feeders, the initial, minimum, and maximum SOCs are 0.5, 0.3, and 1.0, respectively, and the charging and discharging efficiencies dch and dch are both 95%. The two exponents of the exponential load model at any node are identically set to d = 1.5 and d = 3.15. The allowable range of voltage magnitude for all nodes in the system is [0.95 p.u., 1.05 p.u.]. We assume that there is an identical demand reduction XdDR = 0.1 at any node d. All sensitivity results are based on the optimal solution that is calculated at 9 A.M. by VVO, and the optimal voltage magnitudes in the IEEE 33-node test system are shown in Fig. 4. 4.2. Sensitivity of the optimal voltage magnitude and real/reactive power flow at the substation to data changes In this subsection, the sensitivities of the feeder voltage magnitudes and the injected real and reacitve power at the substation are evaluated with respect to the changes in the (G1) (G3) data. To better quantify the voltage magnitude sensitivities, 3 out of 33 nodes are selected: node 11 with the PV system/ESS and two end nodes 18 and 33 without the PV system/ESS for different feeder sections. Fig. 5 show and compare the sensitivities of V11, V18, and V33 to data for line resistance/reactance, demand reduction, and load exponent for voltage-dependent load model. We first observe from Fig. 5(a) that the change in line parameter results in the sensitivities with negative value. This result is due to the increase of the line resistance and reactance values leading to a higher voltage drop according to the voltage Eq. (11). We also observe that the change of line parameters yields the voltage sensitivities with an opposite sign at both end nodes of the line. Specifically, when the parameter at some line increases, the voltage magnitude decreases significantly at the downward node connected to the line whereas the voltage magnitude increases marginally at the upward node connected to the line. For example, V11 significantly decreases with an increasing parameter value of line 10 whereas V11 marginally increases with an increasing parameter value of line 11. In addition, we can list from Fig. 5(a) the top five most influential line data to three voltage magnitudes in a decreasing order of sensitivity as follows: r1 (x1) > r2 (x2 ) > r3 (x3) > r4 (x 4 ) > r5 (x5) . In this order, it is verified that data change in line 1 among all lines has the most significant impact on voltage magnitude. This order of sensitivity can be used to prioritize a distribution system upgrade for the design of reliable network planning. For example, the upgrade of line 1 can be considered first to maintain a robust VVO operation to the potential line data corruption under a certain operating condition. Another observation is that in general the change of line resistance has a larger impact on voltage magnitude than that of line reactance. This is because the amount of the scheduled real power flow is greater than the scheduled reactive power flow by the VVO and hence

Fig. 5. Sensitivities of Vt, d at 9 A.M. with respect to the change in: (a) rdf and x df , (b) XdDR , and (c) d and d .

the voltage magnitude becomes more distorted with the corruption of line resistance data than with the corruption of line reactance data in (11). Fig. 5(b) illustrates the sensitivities of V11, V18, and V33 with respect to the change of demand reduction XdDR at node d. We observe from this figure that any sensitivity value to demand reduction at any node is positive. This observation is consistent with our expectation that as the load consumption at node d decreases with an increasing XdDR , the voltage magnitude at node d increases. Moreover, the increase in 7

Applied Energy 261 (2020) 114331

D. Mak, et al.

demand reduction at some node can increase the voltage magnitude at the other nodes even without demand reduction. On the other hand, the change of demand reduction at some node does not always affect the voltage magnitude at the same node the most. For example, in Fig. 5(b), we can identify the most influential demand reduction-enabled node to the voltage magnitude change: node 30 for V11 and V33 and node 14 for V18 . In particular, it is observed that demand reduction at node 30 has a quite large impact on V11, V18, and V33. Indeed, the sensitivity associated with demand reduction relies on the size of load at each node. In this simulation, a larger load is connected to node 30 rather than the other nodes and hence the load change at node 30 leads to significant voltage distortions at nodes 11, 18, and 33. Fig. 5(c) shows the sensitivity of V33 to the change of real and reactive load exponent at any node. Compared to the sensitivity results regarding demand reduction change, the changes of load exponent values have a relatively smaller impact on the voltage magnitude. In this figure, we can classify all nodes into two sensitivity groups with a different sign of sensitivity: Group I with a positive sensitivity (nodes 7–18, nodes 27–33); and Group II with a negative sensitivity (nodes 2–6, nodes 19–26). The sign of the sensitivity is determined by the value of the optimal voltage magnitude from VVO, as shown in Fig. 4. For example, if some node has the optimal voltage magnitude less than 1 p.u., the increase of its load exponent value results in the decrease of its load consumption, consequently leading to the increase of the voltage magnitude (Group I). Otherwise, the increase of load exponent value at node leads to the decrease of the voltage magnitude at node (Group II). That is, the number of nodes for each sensitivity group depends on the optimal voltage magnitude that is calculated by VVO at every dispatch interval. In addition, similar to the result of Fig. 5(b), the load exponent data at node 30 has the most significant impact on the voltage magnitude sensitivity. Figs. 6 show and compare the sensitivities of real power flow P1line and reactive power flow Q1line at the substation to data for line resistance/reactance, demand reduction and the exponent for the exponential load model. In these figures, the left and right y-axis correspond to the sensitivity of the P1line and Q1line , respectively. From Fig. 6(a), we can verify the negative sensitivities at most nodes. This implies that the increase of the line resistance or reactance causes P1line and Q1line to decrease. This is because the increase of the line parameter results in a more voltage drop as shown in Fig. 5(a), which in turn reduces real and reactive power consumption at each load node and hence requires less real and reactive power injected from the substation. Two observations from Fig. 5(a) still hold true for the sensitivity analysis of P1line and Q1line : (i) the parameter of line 1 is the most influential to the sensitivity of the real and reactive power flows at the substation, and (ii) the sensitivity to the line resistance data is larger than the sensitivity to the line reactance data. Fig. 6(b) shows the impact of demand reduction XdDR on the P1line and Q1line . This figure shows that two sensitivity plots show similar behavior to each other where the sensitivity plot of the P1line is lower than that of the Q1line . This observation derives from the fact that more demand reduction requires less real power from the grid and allows consumers to reduce their electricity bills successfully. Fig. 6(c) shows the sensitivity of P1line and Q1line to the change of load exponents in real and reactive load models. As expected, it is observed from this figure that the sensitivities of P1line and Q1line to their corresponding d and d , respectively, are more fluctuating than the other two types of sensitivities that show almost identical values at any node. Similar to the observation from Fig. 5(c), the sensitivity of P1line also has positive or negative value at each node. This is because the sign of this sensitivity relies on the optimal voltage magnitude from VVO. When the optimal voltage magnitude at a certain node is smaller than 1 p.u., the increase of load exponent at the node results in the decrease of the load consumption, consequently leading to the decrease of real and reactive power flow at the substation. With this reason, the sensitivity of the P1line has opposite sign compared to the sensitivity of voltage in Groups I

Fig. 6. Sensitivities of P1line and Q1line at 9 A.M. with respect to the change in: (a) rdf and x df , (b) XdDR , and (c) d and d .

and II. Table 3 shows the sensitivity of V11, V18, V33 and P1line, Q1line to the change in the predicted PV real power output. From the second, third and fourth column of this table, we observe that all sensitivity values of three voltage magnitudes are positive. This observation is natural because the PV real power injection into the grid increase the voltage profile along the distribution feeder. In addition, we also verify the localized effects where the sensitivity of the voltage magnitude at a certain node increases with the increase in data for PV adjacent to the node. For example, in the third column of Table 3, the change of PV data at node 16 has the most influential impact on the sensitivity of V18 (bolded number) where node 16 is the nearest to node 18 among the four nodes with the PV system. This localized effect is verified for the other sensitivities of V11 and V33 as well. In the fifth and sixth column of Table 3, we identify that the sensitivities of P1line and Q1line are negative 8

Applied Energy 261 (2020) 114331

D. Mak, et al.

to the demand reduction at nodes 20 22 are positive, and the total voltage magnitude deviation becomes greater with an increasing demand reduction at these nodes. That is, demand reduction at these nodes may reduce consumers’ electricity cost at the expense of degrading the voltage quality along the feeder. In addition, it is verified from Fig. 7(b) that most of nodes have negative sensitivities for J in regards to load exponents. In summary, the aforementioned sensitivity results can be used by system operators to more flatten the feeder voltage profile by tuning these data for line resistance/reactance, demand reduction, and load exponents.

Table 3 Sensitivity of voltage magnitude and real/reactive power flow at the substation at 9 A.M. to predicted PV real power output. Data

V11/ a

V18/ a

V33/ a

P1line/ a

Q1line/ a 0.07419

PV

0.04091

0.03991

0.01510

−0.94155

PV P16 PV P27 PV P31

0.04013

0.06927

0.01482

−0.93099

0.08311

0.01584

0.01546

0.01480

−0.95588

0.06502

0.01517

0.01480

0.04491

−0.93987

0.10058

P11

4.4. Sensitivity of the total real power loss to data changes

and positive, respectively. Indeed, more PV real power injection into the grid decrease the injected real power at the substation, consequently yielding a negative sensitivity for P1line . In contrast, the sensitivity of the Q1line is positive, and the increase in PV real power output increases the voltage magnitude, which causes a more reactive power consumption of loads. Regarding the PV data change, the most influential nodes to the sensitivity of real and reactive power flow at the substation are identified at nodes 27 and 31, respectively.

In this subsection, we conduct another case study where the objective function of the VVO in Section 2.2 is modified to the total real power loss in the network. This implies that the proposed sensitivity framework can be applicable to a general optimization problem where its KKT conditions exist. The VVO problem of this case study has the same constraints as the one in Section 2.2. The only difference is a modified objective function, which is expressed by

4.3. Sensitivity of the optimal total voltage deviation to data changes

minJ Loss =

rdf t = 1 df

In this subsection, the sensitivities of the optimal total voltage deviation J to the changes in (G1) and (G3) data are illustrated in Figs. 7. Fig. 7(a) shows the sensitivity of J with respect to the change of the resistance and reactance for all lines. We classify the sensitivities into two line groups: i) Group I with a negative sensitivity (lines 1, 18 24); and ii) Group II with a positive sensitivity (lines 2–17, 25–32). Using the sensitivity results, system operators can, from a perspective of the minimization of the total voltage deviation, upgrade the distribution system by increasing and decreasing the value of line data in Groups I and II, respectively. In this figure, we also identify the top three most influential lines to the value of J and list them in decreasing order of sensitivity: r5 > r2 > r4 for line resistance and x5 > x 4 > x3 for line reactance. This order information can also be used to prioritize a distribution system upgrade for maintaining good quality of the voltage profile. Fig. 7(b) shows the sensitivity of J with respect to the data change in the load model such as demand reduction and load exponents. In this figure, the left and right y-axis correspond to the sensitivity regarding demand reduction and load exponents, respectively. We can verify from Fig. 7(b) that demand reduction at most of nodes reduces the total deviation of voltage magnitudes with negative sensitivities. It is noted that demand reduction at node 30 yields the lowest sensitivity, which implies that among all nodes, demand reduction at node 30 can minimize the total voltage magnitude deviation and flatten feeder voltage profiles the most. However, the sensitivities with regard

2

Nt

line (Ptline , df ) + (Qt , df )

(Vt ,1)2

2

.

(38)

It is noted that the modified VVO formulation neglects power loss terms in its power flow and voltage drop equations (Eqs. (9)–(11)). However, based on the results in [31,32], it is reasonable to formulate the total real power loss as the objective function since it is written in term of real and reactive power flow of each line obtained from Eqs. (9)–(11). This case study is validated in IEEE 33-node test system under the same simulation environment as described in Section 4.1. Fig. 8(a) illustrates the sensitivity of the total real power loss to the changes of resistance and reactance of all lines in IEEE 33-node test system. Unlike the sensitivity results of voltage deviation to these two parameters, all the sensitivities of the total real power loss are positive. This is natural because an increase in line resistance and reactance results in an increase in the total real power loss in the network. We can also observe from Fig. 8(a) that these parameters of line 1 impact the most on the real power loss. In addition, it is observed that a change in resistance has a larger impact than in reactance. This is due to the fact that real power loss is directly related to line resistance. Similar to Fig. 7(b), Fig. 8(b) shows the sensitivity of J Loss with respect to change of data in load models including demand reduction and load exponents. Regarding a change in demand reduction, the sensitivity values shown in the left y-axis are all negative, which mean that an increase in demand reduction (a decrease in loads) yields a lower real power loss. This observation is natural and easily expected, however, this

Fig. 7. Sensitivities of J at 9 A.M. with respect to the change in: (a) rdf and x df , and (b) XdDR , 9

d

and

d.

Applied Energy 261 (2020) 114331

D. Mak, et al.

Fig. 8. Sensitivities of J Loss at 9 A.M. with respect to the change in: (a) rdf and x df , and (b) XdDR ,

sensitivity result can provide an accurate assessment on how much demand reduction influences the power loss in the network. Similarly, we can also see from Fig. 8(b) that most of the sensitivities values of J Loss with respect to load exponents d and d are negative in the right yaxis except the sensitivities at nodes 2 and 19. This is because the optimal voltage magnitudes at these two nodes are above 1 p.u., which increases power consumption of loads when load exponents increase, subsequently leading to an increase of power losses. On the other hand, the optimal voltage magnitudes at the other nodes are below 1 p.u., which yields a contrary effect when load exponent increases and hence reduces the power losses in the network. These sensitivity results can provide some guideline to system operators for assessing the impacts of these data changes on the total real power loss in the network as well as for upgrading or tuning the data in order to improve the reliable distribution network operation in terms of power loss.

d

and

Table 4 Average of the absolute normalized sensitivities

d.

x j

and

33-node system (all values are multiplied by 10 4 ). Data

Variables (

rdf x df

d d dch d load,nom Pd load,nom Qd ESS,dch,max Pd QdPV,max QdESS,max EdESS,max

In this subsection, we quantify the average of the absolute normalized sensitivity of the optimal solution and objective function in VVO to changes in all types of VVO data in both the IEEE 33-node and 123-node test feeders. For the IEEE 33-node test feeder, Table 4 provides two types of the average of the absolute normalized sensitivity, (x ) and (jJ ) , regarding: i) the optimal solution V33, P1line, Q1line , and j SOC11, and ii) the total voltage deviation J subject to all types of VVO data that are defined in Table 2. In this table, the bolded number at each column represent the largest sensitivity value regarding its corresponding data among the selected thirteen types of data. We first observe from the second and sixth columns of Table 4 that the change in rdf and x df among the selected data has the most and second most influential impact on both V33 and J, respectively. Under this VVO scheduling condition, system operators may consider the line upgrade with the highest priority to ensure a good quality of voltage profiles throughout the VVO process. In the third and fourth columns of Table 4, it is observed that the P1line and Q1line are influenced more by the change of XdDR than by that of the other data. This is natural because the change of load demand directly affects the real and reactive power injected from the grid, leading to consumers’ cost savings and voltage regulation, respectively. This result may offer guidelines for a system operator to prioritize the tuning of data for guaranteeing an efficient and robust demand response program while the VVO is carried out. Lastly, regarding the sensitivity of SOC11, the fifth column of Table 4 shows that SOC11 is dominantly affected by only three data with decreasing order of sensitivity as follows: ddch > PdESS,dch,max = EdESS,max . Furthermore, these data have a more significant impact on the SOC11 sensitivity among the selected sensitivity target variables as shown in

P1line

10.95

27.84

0.07

2.97

4.76

PV Pd XdDR

4.5. Normalized sensitivity of the optimal voltage magnitude, real/reactive power, and total voltage deviation to changes in all types of data

V33

x j)

Q1line 38.77

J j

at 9 A.M. in the IEEE

Objective Function (

SOC11

J

0.00

4.93

0.00

0.04

8.92

11.85

0.00

2.15

85.43

48.24

0.00

0.90

0.28 0.35

17.83 0.66

0.93 24.03 0.00

0.00 0.00

710.53

0.12 0.13

0.12

7.83

0.45

0.00

0.06

0.07

0.14

4.79

0.00

0.03

0.06

2.54

0.22

7.11

0.03

0.07

0.17

4.12

0.00

0.04

0.02

0.05

1.24

0.00

0.01

0.00

0.00

0.00

7.11

0.00

0.00

0.00

0.25

J j)

1.70

0.00

the ninth, twelfth, and fifteenth rows of Table 4. These sensitivity results inform that, through tuning only ESS related data in the VVO process, a system operator can maintain an adequate SOC level of ESS without being concerned about their impacts on overloads of the distribution line and nodal voltage violation. For the IEEE 123-node test feeder, Table 5 provides two types of the average of the absolute normalized sensitivity, (jx ) and (jJ ) , regarding: (i) the optimal solution V114, P1line, Q1line , and SOC97 , and (ii) the total voltage deviation J subject to all types of VVO data. It is verified that all the observations and their interpretations from Table 5 are consistent with those from Table 4. Finally, the normalized sensitivity indices (36) and (37) can be modified to provide system operators with system-wide operational or planning metric to assess the impact of data on VVO, on average, with respect to a total of j data and T execution intervals for VVO. (x ) j (T )

T

j

=

ai

x (t ) / ai

j

T

ai

J (t ) / ai

j

T

t=1 i =1

(J ) j (T )

T

j

= t=1 i =1

(39)

(40)

Here, the indices (39) and (40) with T = 1 and T 1 correspond to system-wide operational or planning index, respectively. 10

Applied Energy 261 (2020) 114331

D. Mak, et al.

Table 5 Average of the absolute normalized sensitivities

x j and 4

123-node system (all values are multiplied by 10 ). Data

rdf x df PV

Pd

XdDR d d dch d load,nom Pd load,nom Qd ESS,dch,max Pd QdPV,max QdESS,max EdESS,max

Variables (

V114

P1line

0.657

1.842

0.538

x j)

J j

Objective Function (

Q1line 1.973

SOC97

J

0.00

0.183

0.627

0.636

0.00

0.084

5.419

3.378

0.00

0.019

0.009 0.012

1.141 0.042

0.066 1.655 0.00

0.00 0.00

22.453

0.001 0.003

0.004

0.525

0.030

0.00

0.029

0.002

0.009

0.446

0.00

0.016

0.00

0.00

0.00

0.00

0.00

0.002

0.009

0.389

0.00

0.001

0.001

0.003

0.256

0.00

0.002

0.00

0.00

0.00

2.245

0.00

0.00

0.00

0.00

0.00

0.00

0.00

Table 6 Mean absolute percentage error of the proposed method at 9 A.M. in the IEEE 33-node system with respect to various types of data (all values are in %).

at 9 A.M. in the IEEE

j

j

× i=1

Si f ,per ( ) Si

f ,per

0.563

( = 0.0001)

XdDR ( = 0.1) d ( = 0.1) d ( = 0.1)

0.00

dch ( d load,nom Pd ( load,nom Qd ( PdESS,dch,max QdPV,max ( QdESS,max ( EdESS,max (

0.00

= 0.01) = 0.0001) = 0.0001) ( = 0.001)

= 0.001) = 0.001) = 0.001)

0.426

P1line

Q1line

0.820

1.517

0.928

1.522

SOC11 0

J 13.625

0

12.813

0.085

0.026

0.812

0

0.089

0.039

0.165

0

0.644

0.231 0.232

0.171 0.015

0.377 0.173 0

0 0

4.302

12.405 12.414

0.060

0.00006

0.548

0

11.765

0

0

2.776

0

0.061

0.187

0.0008

0

11.977

9.932

9.981

10.065

5.086

12.379

0.091

0.424

0.043

0

2.328

0.094

0.421

0.046

0

2.325

0

0

0

4.722

0

type j data, Si f ,pro is the sensitivity obtained from the proposed method, and j is a total number of type j data. Tables 6 and 7 show the values of MAPE for the proposed method in the IEEE 33-node and 123-node systems with different amount of data change ( = 0.0001, 0.001, 0.01, 0.1), respectively. We observe from these tables that the MAPEs for both test systems are quite small. Another observation is that in general the MAPE for the objective function J is higher than for the decision variables. This observation is natural because the MAPE for the objective function J involves the sum of the MAPEs for voltage magnitude sensitivities for all nodes, which consequently leads to a higher MAPE. This phenomenon would become more prominent with larger scale distribution systems. Next, Table 8 shows the total computation time (i.e., the sum of the computation times for MINLP and NLP) in the IEEE 33-node and 123node test systems. Since the computation time for calculating the sensitivity using (35) is very small, it is excluded in the total computation time. The results in Table 8 are calculated in computer (Intel Core i5 CPU 3.0 GHz with 8 GB of RAM) using BONMIN solver via OPTI toolbox version 2.28. From Table 8, we can see that the total computation times are around three minutes (190 s) and ten minutes (615 s) for the IEEE 33-node and 123-node test systems, respectively. This result justifies that in view of the total computation time, the proposed sensitivity framework can be well applied to the hourly-based VVO algorithm formulated in this paper. Considering a high performance computer in distribution management system, the proposed sensitivity framework can also be applicable to VVO with a faster execution time (e.g., 15 minbased VVO algorithm) in large-scale distribution networks.

Si f ,pro ( )

rdf ( = 0.1) x df ( = 0.1) PV

In this subsection, the accuracy and computation time of the proposed sensitivity framework are quantified and validated. To assess the accuracy of the proposed approach, the perturbation method is selected as a benchmarking method, which is expressed as f (x ) ([f (x + ) f (x )]/ ) . Fig. 9(a) and (b) compare voltage magDR DR = 0.1) between nitude sensitivities regarding the change in X11 ( X11 the perturbation method and the proposed method in the IEEE 33-node and 123-node test systems, respectively. To better clarify the comparison of the sensitivities between both methods in the IEEE 123-node system, the sensitivities from node 1 to node 20 are illustrated as shown in Fig. 9(b). We can observe from Fig. 9(a) and (b) that the sensitivities of the perturbation method and the proposed method in both test systems are almost identical. To show the overall accuracy of the proposed approach compared to the conventional perturbation approach, the following mean absolute percentage error (MAPE) is selected as a performance index

100%

V33

Pd

0.178

4.6. Validation of the accuracy and computation time of the proposed sensitivity framework

MAPEfj ( ) =

Data

J j)

(41)

where Si f ,per ( ) is the sensitivity of the VVO optimal solution f ( f = x orJ ) from the perturbed method with respect to change in

Fig. 9. Comparison of voltage magnitudes sensitivities at 9 A.M. between the proposed framework and the perturbation method with respect to change in DR DR X11 ( X11 = 0.1) in: (a) the IEEE 33-node test system, and (b) the IEEE 123-node test system. 11

Applied Energy 261 (2020) 114331

D. Mak, et al.

• The total computational times for executing the VVO and calcu-

Table 7 Mean absolute percentage error of the proposed method at 9 A.M. in the IEEE 123-node system with respect to various types of data (all values are in %). Data

V114

rdf ( = 0.1) x df ( = 0.1)

1.045

1.113

1.282

1.918

1.582

0.254 0.431

PV

Pd

1.124

dch ( d load,nom Pd ( load,nom Qd ( PdESS,dch,max QdPV,max ( QdESS,max ( EdESS,max (

= 0.0001) = 0.0001) ( = 0.001)

= 0.001) = 0.001) = 0.001)

SOC9

1.223

1.784

0

14.112

5.356

0

1.153

0

1.278

0.197 0.112

0.691 0.287 0

0 0

5.531

12.931 13.193

2.861

0.998

2.952

0

13.923

1.051

1.973

0.849

0

12.877

0

0

0

0

0

0.902

1.003

0.126

0

5.432

1.007

1.404

0.224

0

6.221

0

0

0

5.204

0

0

1.989

J

13.251

0

= 0.01)

Q1line

0

1.342

( = 0.0001)

XdDR ( = 0.1) d ( = 0.1 ) d ( = 0.1)

P1line

lating the sensitivity are around 3 min and 10 min for IEEE 33-node and 123-node test systems, respectively. This result justifies that in view of the total computation time, the proposed sensitivity framework can be well applied to the hourly-based VVO algorithm formulated in this paper.

Finally, in this study, the proposed approach is limited to calculate a local sensitivity of VVO to the change in a small single data corruption without considering multiple data corruptions. However, we emphasize that our study is the first step toward understanding the impact of data corruption on VVO. An important extension of our study here would be to conduct a global sensitivity analysis of VVO considering multiple data corruptions. To this end, a possible direction would use a multiparametric programming [33], and it is referred to as a future work.

9.732

0

5. Conclusions The Volt-VAR optimization is formulated with many different types of input data that involve distribution line parameters, PV and load forecasts with demand reduction, voltage-dependent load models, and operating conditions for PV system/ESS and voltage regulators. Since data change may take place due to a natural error or a man-made attack, maintaining data quality in Volt-VAR optimization is of vital importance for reliable distribution system operations. In this context, sensitivity analysis can help to rapidly assess the impact of data change on the performance of Volt-VAR optimization. In this study, we propose an analytical framework to measure the sensitivity of Volt-VAR optimization with regards to changes in various types of data. We build a linear sensitivity matrix based on the perturbed KKT equations from the nonlinear programming-based VVO formulation that is transformed from the original mixed integer nonlinear programming-based VVO formulation. Using the built sensitivity matrix, we then investigate how much the optimal solution of Volt-VAR optimization (i.e., voltage magnitude and total voltage deviation) varies with changing data. In light of maintaining voltage quality during the Volt-VAR optimization process, the analytical sensitivity results can potentially be used for identification of the most influential data to Volt-VAR optimization through listing different sensitivity values and robust distribution system operation to cyber data attack into Volt-VAR optimization. In the future, we plan to develop the sensitivity analysis framework for Volt-VAR optimization in more realistic three-phase unbalanced distribution networks. The developed sensitivity framework should be validated and tested in larger test systems under various operating conditions. Another interesting future direction is to develop a general sensitivity framework with which system operators can perform a global sensitivity analysis of Volt-VAR optimization considering multiple data corruptions.

Table 8 Total computation time in the IEEE 33-node and 123-node test systems. System

MINLP Computation time (s)

NLP Computation time (s)

Total time (s)

IEEE-33 IEEE-123

180 600

10 15

190 615

This study is the first assessment of the sensitivity of VVO with respect to change in all the VVO data. The advantages and main results of the proposed sensitivity analysis framework can be summarized as follows:

• The developed sensitivity matrix is an online analysis tool to quickly

•

•

quantify the impact of data on the performance of the VVO. This sensitivity approach requires no repetitive execution of VVO with varying data. Furthermore, from an implementation perspective, the developed sensitivity framework can be readily integrated as an analysis tool into VVO without the need to change the VVO formulation. The sensitivity matrix rapidly informs system operators of the most influential data to the VVO solution. Analytical results based on the sensitivity matrix will provide system operators some guideline to, through ordering the VVO sensitivity results, enhance the performance of VVO from a planning (e.g., line upgrade) and operational (e.g., efficient demand response program) perspective. Furthermore, the sensitivity results can also be used to mitigate the impact of cyber attack on VVO by prioritizing the upgrade of sensors that monitor power consumption and PV generation output, such as smart meters and smart inverters. The accuracy test results for the proposed approach show that the mean absolute percentage errors of the sensitivity of voltage magnitude and the total real/reactive power flow are mostly below 1% for IEEE 33-node test system and around 1 2% for IEEE 123-node test system, respectively. On the other hand, the maximum mean absolute percentage errors of sensitivity of the objective function to all data are 12.4% for IEEE 33-node test system and 14.1% for IEEE 123-node test system, respectively. These errors derive from the fact that the objective function J involves the sum of the mean absolute percentage errors for voltage magnitude sensitivities for all nodes.

Declaration of Competing Interest None. Acknowledgments This work was supported in part by the Korea Electric Power Corporation under Grant R17XA05-75, and in part by the Korea Government (MSIP) through the National Research Foundation of Korea (NRF) under Grant 2018R1C1B6000965.

12

Applied Energy 261 (2020) 114331

D. Mak, et al.

Appendix A. Perturbation submatrices The submatrices in (35) are expressed in the following two subsections. In these subsections, Diag(x )a represents the a × a diagonal matrix with an element x. I(a × b) and 0(a × b) are the a × b identity and zero matrices, respectively. Np and Na are a total number of decision variables and data, respectively. Nb and N are a total number of equality and binding inequality constraints, respectively. Here, Nb and N include the number of equality and binding inequality constraints corresponding to their Lagrangian multiplier, respectively. For example, Nb is a total number of equality constraints associated with Lagrangian multiplier . N is a total number of binding inequality constraints associated with Lagrangian multiplier . A.1. Submatrices on the Left-hand Side of (35) The following matrix corresponds to the perturbed objective functions with respect to the decision variables:

V nom )(1 × N N ) 0 (1× [Np

M(1 × Np) = [2(Vt , d

t d

The other three matrices ( 1, Lagrangian multipliers: 1(Np × Np)

Diag(z ) Nt Nd

=

0 ([Np

3)

2,

0 (Nt Nd × [Np 0 ([Np

Nt Nd] × Nt Nd )

Nt Nd ]) ].

(42)

are associated with the perturbed KKT condition equations with respect to the decision variables and Nt Nd])

Nt Nd ] × [Np Nt Nd])

where load,nom (1 t , d Pt , d

z=2+[

+[

load,nom (1 t , d Qt , d

XtDR ,d )

XtDR ,d )

d (Vt , d )

d (Vt , d )

d 2(

d 2(

(43)

1)]

d

(44)

1)]

d

and 2(Nb × Np) =

[

T 20

T 21

T 22

T 23

T 24

I(Nb

× Nt Nd )

T 25

T 26

T 27

T 28

T 29 ]

(45)

where

0

20

= [ 0 (Nb

× Nt Nd )

21

= [ 0 (Nb

× 2Nt Nd )

22

= I

24

× Nt Nd )

0 (Nb

(

)

rdf V nom

(46)

× 10Nt Nd ) ]

(47)

× 9Nt Nd) ]

Nb × Nt Nd

Diag

(

xdf V nom

)

Nb × Nt Nd

]

1

=

Diag

Nb × Nt Nd

Nb × 9Nt Nd

23

I(Nb

0 (Nb

(48) Nb × Nt Nd

= [0(Nb × 2Nt Nd)

0

Nb × 11Nt Nd

I(Nb × Nt Nd)

0(Nb × Nt Nd)

(49)

0(Nb × 2Nt Nd) I(Nb × Nt Nd)

I(Nb × Nt Nd)

0(Nb × 4Nt Nd) ]

(50)

25=

[0(Nbµ × 4Nt Nd) I(Nbµ × Nt Nd)

I(Nbµ × Nt Nd) 0(Nbµ× Nt Nd) I(Nbµ × Nt Nd) 0(Nbµ× Nt Nd) I(Nbµ× Nt Nd) 0(Nbµ × Nt Nd)

I(Nbµ × Nt Nd) ]

(51)

26= load, nom [Diag (P^t , d (1

I(Nb

× Nt Nd )

0(Nb

XtDR ,d )

d (Vt , d )

d 1)

(Nb × Nt Nd )

0(Nb

× 4Nt Nd)

× 6Nt Nd ) ]

(52)

27= load, nom [Diag (Q^t , d (1

I(Nb 28

0(Nb

× 5Nt Nd ) ]

× 11Nt Nd)

I(Nb

× Nt Nd )

= [0(Nb

XtDR ,d )

d (Vt , d )

d 1) (N × N N ) b t d

0(Nb

× 5Nt Nd )

(53) × Nt Nd )]

(54)

13

Applied Energy 261 (2020) 114331

D. Mak, et al.

= [0

29

I

Nb × 7Nt Nd

1

Diag

0

dch ESS,max d Ed

0

Nb × Nt Nd

Nb × 2Nt Nd

Nb × Nt Nd

Nb × Nt Nd

(55)

and 3(Nk × Np)

=[

T

T ]T

T

(56)

where

= [ 0 (Nk

× 9Nt Nd )

I(Nk

× Nt Nd )

0 (Nk

× 2Nt Nd )) ]

(57)

= [ 0 (Nk

× 8Nt Nd )

I(Nk

× Nt Nd)

0 (N k

× 3Nt Nd ) ]

(58)

= [ 0 (N k

× 7Nt Nd )

I( N k

× 4Nt Nd ) ].

(59)

0 (Nk

× Nt Nd)

A.2. Submatrices on the right-hand side of (35) The following matrix correspond to the perturbed objective functions with respect to data: (1 × Na)

(60)

= 0.

The former is the row vector with zeros because the total voltage deviation function has no data associated with the sensitivity analysis. The other three matrices ( 2, 3, 4 ) are associated with the perturbed KKT condition equations with respect to data 2(Np× Na)

=[

T 21

T 22

0T

T 23

T 24

0T

0T ]T

(61)

where 21

= [0 (Nt Nd × [2Nd + Nt Nd]) Diag(k1) (Nt Nd × Nt Nd) Diag(k2)(Nt Nd × Nt Nd) (62)

Diag(k3)(Nt Nd × Nt Nd) Diag(k4 )(Nt Nd × Nd) Diag(k5) (Nt Nd × Nd) 0 (Nt Nd × 5Nd])] where

k1 =

t , d (1

XtDR ,d )

d (Vt , d )

k2 =

t , d (1

XtDR ,d )

d (Vt , d )

d 1 d 1

load,nom d 1 t , d Qt , d d (Vt , d )

load,nom d 1 t , d Pt , d d (Vt , d )

k3 = k4 =

load,nom (1 t , d Pt , d

XtDR , d )[(Vt , d )

d 1

+

d (Vt , d )

k5 =

load,nom (1 t , d Qt , d

d 1 XtDR , d )[(Vt , d )

+

d (Vt , d )

Diag

( ) t, d

22

=

23

= 0 (Nt Nd × Nd)

24

= 0 (Nt Nd × [4Nd + 4Nt Nd]) Diag

Diag

( ) t ,d

V nom (N N × N ) t d d

t, d dch ESS,max )2 d (Ed (Nt Nd × Nd )

t , d )]

d 1log(Vt , d )]

(63)

0 (Nt Nd × [8Nd + 4Nt Nd])

V nom (N N × N ) t d d

Diag

d 1log(V

(64)

0 (Nt Nd × [4Nt Nd + 7Nd])

(65)

t ,d 2 ( ddch ) EdESS,max

(Nt Nd × Nd)

0 (Nt Nd × 3Nd)

(66)

and 3(Nb × Na)

=[

T 30

T 31

T 32

T 33

T 34

T 35

T 36

T 37

T 38

T T 39 ]

(67)

where 32

34

=

Diag

= [ 0 (Nb × 2Nd)

Ptline , df V nom

Diag Nb × Nd

Qtline , df V nom

0 Nb × Nd

Nb × 4Nt Nd + 7Nd

(68)

I(Nb × Nt Nd) 0 (Nb × [3Nt Nd + 7Nd]) ]

(69)

14

Applied Energy 261 (2020) 114331

D. Mak, et al. 36=

[0(Nb

× [2Nd + Nt Nd])

Diag (

load, nom Diag (P^t , d (1

37

= [0

load,nom

30

= 0

=

× Nt Nd )

0(Nb

× Nt Nd )

b × Nt Nd )

PtESS,dch ,d dch ESS,max )2 d (Ed

=

33

=[

T 41

0 (Nb

0 Nb × Nd

=

35

=

38

× 6Nd) ]

(70)

× Nd) b × Nd )

Diag

0(Nb

b × Nt Nd )

XtDR , d )(Vt , d ) d log(Vt , d ))(N

Nb × 4Nd + 4Nt Nd

31

× Nd )

XtDR , d )(Vt , d ) d )(N

Diag((1

(Vt , d ) d)(N

load,nom Diag(Qt , d (1

Diag

d XtDR , d )(Vt , d ) log (Vt , d )) (Nb

Nb × 2Nd + Nt Nd

Diag( Qt , d

39

d XtDR , d )(Vt , d ) )(Nb

Diag ((1

load, nom P^t , d (Vt , d ) d ) (Nb × Nt Nd)

PtESS,dch ,d 2 dch ( d ) EdESS,max

0 (Nb

× 5Nd)]

(71)

Nb × Nd

Nb × 3Nd

(72) (73)

=0

and 4(Nk × Na)

T 42

T T 43 ]

(74)

where 41

= [ 0 (N k

× [7Nd + 4Nt Nd])

I(Nk

× Nd )

42

= [ 0 (N k

× [8Nd + 4Nt Nd])

I( N k

× Nd ) ]

43

= [ 0 (N k

× [6Nd + 4Nt Nd])

I(Nk

× Nd)

0 (Nk

(75)

× Nd) ]

(76)

0 (Nk

(77)

× 2Nd) ].

References [13]

[1] Veda S, Wu H, Martin M, Baggu M. Developing use cases for the evaluation of adms applications to accelerate technology adoption. In: 2017 Ninth Annual IEEE Green Technologies Conference (GreenTech); 2017. p. 132–8. doi:https://doi.org/10. 1109/GreenTech.2017.25. [2] Ahmadi H, Martí JR, Dommel HW. A framework for volt-var optimization in distribution systems. IEEE Trans Smart Grid 2015;6(3):1473–83. https://doi.org/10. 1109/TSG.2014.2374613. [3] Dao TV, Chaitusaney S, Nguyen HTN. Linear least-squares method for conservation voltage reduction in distribution systems with photovoltaic inverters. IEEE Trans Smart Grid 2017;8(3):1252–63. https://doi.org/10.1109/TSG.2016.2536782. [4] Weckx S, Gonzalez C, Driesen J. Combined central and local active and reactive power control of PV inverters. IEEE Trans Sustain Energy 2014;5(3):776–84. https://doi.org/10.1109/TSTE.2014.2300934. [5] Wang L, Bai F, Yan R, Saha TK. Real-time coordinated voltage control of PV inverters and energy storage for weak networks with high PV penetration. IEEE Trans Power Syst 2018;33(3):3383–95. https://doi.org/10.1109/TPWRS.2018.2789897. [6] Ji H, Wang C, Li P, Zhao J, Song G, Ding F, et al. A centralized-based method to determine the local voltage control strategies of distributed generator operation in active distribution networks. Appl Energy 2018;228:2024–36. https://doi.org/10. 1016/j.apenergy.2018.07.065. [7] Zhang W, Gandhi O, Quan H, Rodríguez-Gallegos CD, Srinivasana D. A multi-agent based integrated volt-var optimization engine for fast vehicle-to-grid reactive power dispatch and electric vehicle coordination. Appl Energy 2018;229:96–110. https:// doi.org/10.1016/j.apenergy.2018.07.092. [8] Manbachi M, et al. Impact of EV penetration on volt–var optimization of distribution networks using real-time co-simulation monitoring platform. Appl Energy 2016;169:28–39. https://doi.org/10.1016/j.apenergy.2016.01.084. [9] Isozaki Y, et al. Detection of cyber attacks against voltage control in distribution power grids with PVs. IEEE Trans Smart Grid 2016;7(4):1824–35. https://doi.org/ 10.1109/TSG.2015.2427380. [10] Teixeira A et al. Security of smart distribution grids: Data integrity attacks on integrated volt/var control and countermeasures. In: 2014 American Control Conference (ACC); 2014. p. 4372–8. https://doi.org/10.1109/ACC.2014.6859265. [11] Choeum D, Choi D-H. OLTC-induced false data injection attack on Volt/VAR optimization in distribution systems. IEEE Access 2019;7:34508–20. https://doi.org/ 10.1109/ACCESS.2019.2904959. [12] Zhou Q, Tesfatsion L, Liu CC. Global sensitivity analysis for the short-term

[14] [15]

[16] [17] [18] [19] [20] [21] [22] [23] [24] [25]

15

prediction of system variables. In: 2010 IEEE PES general meeting; 2010. p. 1–8. https://doi.org/10.1109/PES.2010.5588174. Conejo AJ, Castillo E, Minguez R, Milano F. Locational marginal price sensitivities. IEEE Trans Power Syst 2005;20(4):2026–33. https://doi.org/10.1109/TPWRS. 2005.857918. Mínguez R, Conejo AJ. State estimation sensitivity analysis. IEEE Trans Power Syst 2007;22(3):1080–91. https://doi.org/10.1109/TPWRS.2007.901488. Molin A, Sandberg H, Johansson M. Study on the sensitivity matrix in power system state estimation by using sparse principal component analysis. 2016 IEEE 55th conference on decision and control (CDC) 2016. p. 1529–35. https://doi.org/10. 1109/CDC.2016.7798483. Choi D-H, Xie L. Data perturbation-based sensitivity analysis of real-time lookahead economic dispatch. IEEE Trans Power Syst 2017;32(3):2072–82. https://doi. org/10.1109/TPWRS.2016.2598874. Choi D-H, Xie L. A framework for sensitivity analysis of data errors on home energy management system. Energy 2016;Part1(117):166–75. https://doi.org/10.1016/j. energy.2016.10.062. Kang J-W, Choi D-H, Xie L. Impact of data quality in home energy management system on distribution system state estimation. IEEE Access 2018;6:11024–37. https://doi.org/10.1109/ACCESS.2018.2804380. Jeong H-Y, Choi J-C, Won D-J, Ahn S-J, Moon S-I. Formulation and analysis of an approximate expression for voltage sensitivity in radial DC distribution systems. Energies 2015;8:9296–319. https://doi.org/10.3390/en8099296. Weckx S, D’Hulst R, Driesen J. Voltage sensitivity analysis of a laboratory distribution grid with incomplete data. IEEE Trans Smart Grid 2015;6(3):1271–80. https://doi.org/10.1109/TSG.2014.2380642. Wang X, et al. Optimal voltage regulation for distribution networks with multimicrogrids. Appl Energy 2018;210:1027–36. https://doi.org/10.1016/j.apenergy. 2017.08.113. 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. https://doi.org/10.1016/j.apenergy.2018.08.042. Castillo E, Conejo AJ, Castillo C, Mínguez R, Ortigosa D. Perturbation approach to sensitivity analysis in mathematical programming. J Optim Theory Appl 2006;128(1):49–74. https://doi.org/10.1007/s10957-005-7557-y. Wang Z, Chen H, Wang J, Begovic M. Inverter-less hybrid voltage/var control for distribution circuits with photovoltaic generators. IEEE Trans Smart Grid 2014;5(6):2718–28. https://doi.org/10.1109/TSG.2014.2324569. Baran ME, Wu FF. Network reconfiguration in distribution systems for loss reduction and load balancing. IEEE Trans Power Deliv 1989;4(2):1401–7. https://doi.

Applied Energy 261 (2020) 114331

D. Mak, et al. org/10.1109/61.25627. [26] Price W, Casper SG, Nwankpa C, Bradish R, Chiang H-D, Concordia C, et al. Bibliography on load models for power flow and dynamic performance simulation. IEEE Trans Power Syst 1995;10(1):523–38. https://doi.org/10.1109/59.373979. [27] Tseng K-H, Kao W-S, Lin J-R. Load model effects on ditance relay setting. IEEE Trans Power Deliv 2003;18(4):1140–6. https://doi.org/10.1109/TPWRD.2003. 817507. [28] Huang S, Wu Q, Cheng L, Liu Z. Optimal reconfiguration-based dynamic tariff for congestion management and line loss reduction in distribution networks. IEEE Trans Smart Grid 2016;7(3):1295–303. https://doi.org/10.1109/TSG.2015. 2419080. [29] Kersting WH. Radial distribution test feeders. IEEE Trans Power Syst 1991;6(3):975–85. https://doi.org/10.1109/59.119237.

[30] Jha RR, Dubey A, Liu C-C, Schneider KP. Bi-level Volt-VAR optimization to coordinate smart inverters with voltage control devices. IEEE Trans Power Syst 2019;34(3):1801–13. https://doi.org/10.1109/59.119237. [31] Wang Z, Wang J, Chen B, Begvic MM, He Y. Mpc-based voltage/var optimization for distribution circuit wth distributed generators and exponential load models. IEEE Trans Smart Grid 2014;5(5):2412–20. https://doi.org/10.1109/TSG.2014. 2329842. [32] Hassan A, Mieth R, Chertkov M, Deka D, Dvorkin Y. Optimal load ensemble control in chance-constrained optimal power flow. IEEE Trans Smart Grid 2019;10(5):5186–95. https://doi.org/10.1109/TSG.2018.2878757. [33] Ji Y, Thomas RJ, Tong L. Probabilistic forecasting of real-time LMP and network congestion. IEEE Trans Power Syst 2017;32(2):831–41. https://doi.org/10.1109/ TPWRS.2016.2592380.

16