Nội dung

Journal of Advanced Research 20 (2019) 51–60 Contents lists available at ScienceDirect Journal of Advanced Research journal homepage: www.elsevier.com/locate/jare Original article Optimization unit for real-time applications in unbalanced smart distribution networks M. F. Shaaban a,b,⇑, M. H. Ahmed b, M. M.A. Salama b, Ashkan Rahimi-Kian c a Department of Electrical Engineering, American University of Sharjah, Sharjah 26666, UAE Electrical and Computer Engineering Department, University of Waterloo, Waterloo, ON N2M2C7, Canada c Opusone Energy Corp., Toronto, ON M5C 1P, Canada b h i g h l i g h t s g r a p h i c a l a b s t r a c t
This work proposes a new generalized formulation for the optimal power flow.
The proposed formulation is suitable for real-time energy management systems.
Unlike previous research, the proposed formulation tackles practical unbalanced systems.
Detailed analysis with simulation results of the proposed formulation are provided.
The proposed approach is characterized by scalability and low computational time. a r t i c l e i n f o Article history: Received 12 December 2018 Revised 6 April 2019 Accepted 6 April 2019 Available online 9 April 2019 Keywords: Distribution systems Energy management systems Jacobian matrix Optimal power flow Smart Grids Unbalanced systems Inputs Present and forecasted parameters Real-me unbalanced opmal power flow Objecve funcon Constraints Gradient of the objecve Jacobian of the constraints Hessian of the constraints Outputs Opmal decisions or variables opmal values a b s t r a c t This paper presents a new generic approach for developing a Jacobian matrix for use with the optimization unit in real-time energy management systems (EMS) for unbalanced smart distribution systems. The proposed formulation can replace approximated calculations for real-time optimal power flow in an optimization unit while providing greater accuracy and requiring less computational time, which is critical for real-time EMS. The effectiveness and robustness of the proposed approach have been tested through simulations with different distribution networks. The simulation results demonstrate a significant reduction in the computational time with the new proposed formulation. Moreover, the results demonstrate the scalability of the proposed approach as the reduction in the computational time is more significant for large practical systems. The proposed approach is characterized by evaluating the scalability and low computational time; thus, it can be used by grid operators in real-time energy management applications for large-scale practical distribution systems. Ó 2019 THE AUTHORS. Published by Elsevier BV on behalf of Cairo University. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Introduction For a number of environmental and techno-economic reasons, smart grids have become the subject of significant interest over Peer review under responsibility of Cairo University. ⇑ Corresponding author. E-mail address: mshaaban@aus.edu (M. F. Shaaban). the last few years [1,2]. A particular focus is related to the application of smart grids in distribution systems, where most power system losses and failures occur. For this reason, distribution systems are currently undergoing a significant transition to a new structure with respect to information, control, and power flow. Power flow represents an essential engine in application software for distribution systems [3]. Power flow methods are divided into two categories based on the state variable that is employed for https://doi.org/10.1016/j.jare.2019.04.001 2090-1232/Ó 2019 THE AUTHORS. Published by Elsevier BV on behalf of Cairo University. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). 52 M.F. Shaaban et al. / Journal of Advanced Research 20 (2019) 51–60 solving the power flow problem: either the node voltage or the branch current [3]. In the first category, numerous power flow methods have been developed since the 1960s, including network equivalence methods [4], Newton-Raphson methods [5], and fast decoupled methods [6]. The second category is associated with the development of techniques such as backward/forward sweep methods [7], which are widely used in commercial software [8,9]. Optimizing distribution system assets is a fundamental functionality of smart distribution networks (SDNs). However, even if the convergence of power flow methods is guaranteed, they are not suitable for real-time applications that require an optimal power flow (OPF). These applications require not only a robust OPF but also a fast OPF algorithm, which is the main focus of this work. The OPF problem involves solving for an optimal operating point for a distribution system that minimizes a predefined cost function, such as system losses or generation scenarios that are subject to specific techno-economic and/or environmental constraints on system variables [10]. The OPF usually involves two convergence problems: (1) the convergence of the power flow to a feasible solution and (2) the convergence of the optimization to a global optimum. Many optimization techniques guarantee the feasibility of the solution, which in turn guarantees the convergence of the power flow problem [11]. However, with regard to the optimization and because the unbalance power flow problem is a nonlinear highly nonconvex problem [12], the convergence to a global optimum is not guaranteed. To address this limitation, several methods have been proposed for the global optimization of nonlinear problems, such as convexification [13] and branch-and-reduce [14] methods. OPF convergence problems have thus been adequately addressed in the research. To perform a real-time process, such as energy storage system scheduling, the whole process should be completed in the range of half of a minute up to a few minutes. This process involves not only OPF but also other algorithms, such as topology processing, load allocation, and forecasting. These algorithms are components of the future smart grid energy management systems, which should be designed to address various grid and customer technologies. Most of the recent work in the area of smart grids utilizes optimization software for solving OPF problems, but a major aspect of real-time applications, the computational time required by the proposed algorithms, has not yet been examined. Previous works [15–27] proposed a variety of approaches for coordinating different SDN components/assets in real time: energy storage systems (charging/discharging), electric vehicle (EV) charging, distributed generation (DG), volt/var control equipment, and/or residential load consumption. Some of the abovementioned research has failed to consider power flow [15,16]. Other studies considered the power flow constraints for balanced distribution systems [17–22]. Linearization of the power flow constraints was proposed elsewhere [23,24], where the power flow formulas were linearized around an estimated operation point. This method can lead to inaccurate results, which in turn can have severe consequences in real time. For example, the voltage magnitude or the limits for voltage unbalance could be violated in the actual distribution system on unmonitored nodes but not in the linearized approximated model. The work presented by Maffei et al. [25] proposed a real-time OPF for energy management in smart grids utilizing the preceding horizon technique. The approaches reported previously [26,27] are based on the real-time energy management of controllable loads and EVs in unbalanced distribution systems. However, the approaches described earlier [25–27] utilized commercial software for solving the required optimization problem, which cannot be relied upon in real-time applications, as explained in this section. Moreover, no mention is made of the computational time, which is a highly critical aspect and can be a barrier to the practical implementation of any energy management system. Some work utilized different techniques to solve large scale OPF problems. In Ref. [28], the authors utilized Benders decomposition to reduce the problem size. The work in [29] proposed using a genetic algorithm (GA) and fuzzy clustering to increase the computational speed. Mostafa et al. [30] proposed a multi-objective technique for optimizing the unbalanced system operation under a high penetration of renewable DG by minimizing the system losses and improving the voltage profile. The work presented by others [29,30] relied on the utilization of a GA to solve the mixedinteger nonlinear programming (MINLP) problem. However, meta-heuristic techniques are unsuitable for real-time applications in smart grids for two reasons: (1) the computational time required for these techniques is unpredictable, and (2) for large systems in which the decision variables might be in the range of hundreds or thousands, these techniques are likely to be slow compared to the gradient-based techniques. To solve the unbalanced OPF problem, researchers usually employ optimization software tools, which normally find a solution for the optimization problem using a procedure that includes an evaluation of the approximated values for the constraint derivatives, which is known as the Jacobian matrix. For large networks, the use of an approximated Jacobian matrix, based, for example, on finite differences, significantly increases the computational time. It is well known that providing the derivatives of the constraints to the solver produces a more reliable solution and decreases the computational time. The solver can also find a feasible point for the problem; however, the finite differences around that feasible point may result in an infeasible solution that causes the solver to terminate prematurely without reaching an optimal feasible solution [31]. The work presented in this paper is focused on the optimization unit, which is a core element of any energy management system (EMS). The optimization unit is responsible for solving the assigned OPF problem by the EMS. In this work, we propose a new generalized Jacobian matrix formulation for improving accuracy and reducing the computational time in an unbalanced OPF for SDNs. Energy management system for SDN To optimize the SDN operation in real time, local distribution companies (LDCs) usually adopt a general procedure, whereby a central controller receives data from the users and the metering system via a communication infrastructure. This step is defined as stage 1 in Fig. 1. These data is preprocessed according to the measured values, type of measurement, and time stamp of the measurement; the data is then stored in the database (DB). Further, the data are processed in stage 2 by the forecasting unit to predict future consumptions or generations. In stage 3, all the parameters are sent to the optimization engine to process the real-time OPF to develop the optimal operation decisions of the SDN, which are returned back to the DB. These results are then used in stage 4 for (1) updating the situational awareness about the SDN conditions and (2) updating the control actions of the local controllers (LCs) within the SDN. For real-time SDN system/market operations, this process should be completed as quickly as possible; a few seconds for small systems and up to a few minutes for large systems is considered an appropriate range. The computational time needed for solving the OPF problem represents a major challenge. This is the time it takes for the parameters to be sent to the optimization unit until when the optimal decisions are received, as shown in Fig. 1. The OPF run-time increases exponentially as the size of the system increases. The research presented in this paper addresses Jacobian 53 M.F. Shaaban et al. / Journal of Advanced Research 20 (2019) 51–60 Fig. 1. Main components of a real-time energy management system that utilizes the OPF in an SDN. matrix evaluation, which is a core element for solving the OPF problem. The derivation of the proposed form of the Jacobian matrix for unbalanced distribution systems is explained in the next section.  bði;p Þ ðp Þ ðp Þ ðp Þ ðp Þ ðp Þ 1 Q IN1ðiÞ ¼ Q Gð1iÞ  Q L01ðiÞ V ðiÞ1  Q L11ðiÞ   X X ðp Þ ðp Þ ðp ;p Þ ðp ;p Þ ðp Þ ðp Þ ¼ V ðiÞ1 V ðjÞ2 Y ði;j1Þ 2 sin hði;j1Þ 2 þ dðjÞ2  dðiÞ1 j2I  Proposed generalized OPF ð1Þ ð2Þ hðwÞ ðxðv Þ Þ  0 ð3Þ where v , u, and w are the indices for the variables, equality con
 straints, and inequality constraints, respectively; f xðv Þ is the objec

 tive function; g xðv Þ and hðwÞ xðv Þ are the equality and inequality constraints, respectively. Any OPF problem includes the power flow mismatch constraints expressed in Eqs. (4) and (5) as the equality constraints. The remaining equality and inequality constraints define other technical, environmental, and economical constraints, which are selected based on the system requirements. Eqs. (4) and (5) are as follows:  aði;p Þ ðp Þ ðp Þ ðp Þ ðp Þ ðp Þ 1 P IN1ðiÞ ¼ P Gð1iÞ  P L01ðiÞ V ðiÞ1  P L11ðiÞ   X X ðp Þ ðp Þ ðp ;p Þ ðp ;p Þ ðp Þ ðp Þ ¼ V ðiÞ1 V ðiÞ2 Y ði;j1Þ 2 cos hði;j1Þ 2 þ dðiÞ2  dðiÞ1 u2 ðp Þ ðp Þ ðp Þ ðp Þ ðp Þ 1 1 and Q GðiÞ are the measured at bus i for phase p1 , respectively; P GðiÞ or known generated active and reactive power at bus i for phase ðp Þ g ðuÞ ðxðv Þ Þ ¼ 0 ðp Þ where i and j are the buses indices; p1 and p2 2 U are the buses 1 1 and Q INðiÞ are the injected active and reactive power indices; PINðiÞ ðp Þ 1 1 and Q L0ðiÞ are the nominal measured or known p1 , respectively; PL0ðiÞ and is subject to j2I   ðp ;p Þ ðp Þ ðp Þ sin hði;j1Þ 2 þ dðiÞ2  dðiÞ1 8i 2 I ; p1 2 HðiÞ ð5Þ This section presents the derivation of a Jacobian matrix for the proposed generalized OPF. In general, any OPF problem can be formulated as follows: Minf ðxðv Þ Þ u2 ðp Þ ðu Þ ðp ;p Þ V ðiÞ1 V ðiÞ 2 Y ði;j1Þ 2 ðp ;p2 Þ  V ðiÞ1 V ðjÞ2 Y ði;j1Þ   ðp ;p Þ ðp Þ ðp Þ cos hði;j1Þ 2 þ dðjÞ2  dðiÞ1 8i 2 I; p1 2 HðiÞ ð4Þ active and reactive power of the load at bus i for phase p1 , respecðp Þ ðp Þ 1 1 and Q L1ðiÞ are the nominal unknown active and reactive tively; PL1ðiÞ power of the load at bus i for phase p1 , respectively; aði;p1 Þ and bði;p1 Þ are the active and reactive power exponents of the load at bus i, ðp Þ ðp Þ respectively; V ðiÞ1 and dðiÞ1 are the magnitude and angle of the voltðp ;p2 Þ age at bus i for phase p1 , respectively; Y ði;jÞ1 ðp ;p2 Þ and hði;jÞ1 are the mag- nitude and angle of the admittance element in the branch admittance matrix, respectively; I is the set of all system buses; U ¼ f1; 2; 3g is the set of phases; HðiÞ  U is the subset of the existing phases at bus i. The active and reactive power mismatch equations (4) and (5), are formulated based on a branch admittance matrix, which can be calculated using Kron’s reduction on Carson’s equations for the self and mutual impedances [1]. Since some of the load bus quantities might be measured or known, while others are unknown, it is assumed that the measured or known quantities are voltage dependent. In contrast, the unknown quantities are assumed to be voltage independent because their values are unknown or variable in any case. Solving this nonlinear programming problem requires three sets of derivatives: (1) the gradient, (2) the Jacobian matrix, and 54 M.F. Shaaban et al. / Journal of Advanced Research 20 (2019) 51–60 (3) the Hessian matrices. The focus of the research presented in this paper is the Jacobian matrix parametric definition and the fast numerical computation for the SDN OPF runs. In this study, all the SDN active and reactive powers, voltages, and voltage angles are assumed as variables, and all variables and parameters are considered in per-unit values. It is assumed that the Jacobian matrix size is m  n, where m is the number of constraints and n is the number of variables. Only the power flow constraints and variables are considered in this paper. The Jacobian matrix J can be given as in (6). Thus, each entry J ðu; v Þ in the Jacobian matrix is defined as in (7), which corresponds to the differentiation of constraint u with respect to (w.r.t.) variable v . These considerations mean that each constraint must be differentiated w.r.t. all variables. 2 dg ð1Þ dxð1Þ 6 6 . JðmnÞ ¼ 6 6 .. 4 dg ðmÞ dxð1Þ  .. .  dg ð1Þ dxðnÞ 3 7 7 .. 7 . 7 5 dg ð6Þ ðmÞ dxðnÞ   dg ðuÞ J ðu; v Þ ¼ J g ðuÞ ; xðv Þ ¼ dxðv Þ 8u  m; v n ð7Þ Assuming the total number of buses in the system is N bus , the total number of variables can be defined as 3  6  N bus , where 3 ðp Þ ðp Þ ðp Þ 1 1 , Q L1ðiÞ , P Gð1iÞ , refers to the number of phases, and 6 refers to P L1ðiÞ ðp1 Þ , Q GðiÞ ðp Þ V ðiÞ1 ðp Þ anddðiÞ1 . Moreover, the number of constraints for any OPF can be identified as 3  2  N bus , where 3 refers to the number of phases, and 2 refers to the two type of constraints: the active P ðp Q ðp Þ Þ power mismatch g ðiÞ 1 and the reactive power mismatch g ðiÞ 1 . Therefore, the Jacobian matrix would have 6N bus rows, corresponding to the constraints, and 18N bus columns, corresponding to the variables. The Jacobian matrix thus can be defined as in (8), where   the entry J g ðuÞ ; xðv Þ refers to dg ðuÞ =dxðv Þ . The 6N bus constraints and their derivatives, which represent the Jacobian matrix entries, are defined in (9)–(27), as explained below. As a first step in forming the Jacobian matrix, the Jacobian   entries are initialized to zero, i.e.; J g ðuÞ ; xðv Þ ¼ 08u  m; v  n. This step is included because most Jacobian entries are zeros due to the radial structure of distribution systems; each bus is usually connected to two or three buses. For example, assume that bus i is connected only to buses i þ 1 and i  1. This arrangement means that all the entries for the power mismatch constraint for bus i in the Jacobian matrix will be zeros except for those corresponding to buses i  1; i;and i þ 1. After all the Jacobian entries have been initialized to zeros, as a second step, each entry is updated. The active power mismatch entries are derived first and followed by the reactive power mismatch entries. Jacobian matrix entries for active power constraints To update the Jacobian matrix entries for the active power conP ðp Þ straints, i.e., dg ðiÞ 1 =dxðv Þ , the procedure shown in Fig. 2 can be used, which is explained as follows. Based on (4), each bus has up to three active power mismatch constraints, as in (9). The Jacobian ðp Þ ðp Þ matrix entries corresponding to the variables PL11ðiÞ and PGð1iÞ can thus be derived as in (10) and (11). Pðp Fig. 2. The procedure to update the Jacobian matrix entries for the active power constraints. Þ As indicated in (9), g ðiÞ 1 is composed of two major terms: a positive double summation and a negative double summation. Each summation can be broken down into three components so that P ðp Þ P ðp Þ g ðiÞ 1 can be rewritten as in (12), where g ðiÞ 1 is composed of six terms. The third and sixth terms cancel each other out. 55 M.F. Shaaban et al. / Journal of Advanced Research 20 (2019) 51–60 P ðp Þ For the derivative of g ðiÞ 1 w.r.t. the magnitude of the voltage, three different cases can be defined as follows: Pðp X X þ j2I p2 –p1 Þ
The derivative of g ðiÞ 1 w.r.t. the same bus i and the same phase p1 voltage magnitude, as in (13); P ðp Pðp Þ g ðiÞ 1  aði;p Þ ðp Þ ðp Þ ðp Þ ðp Þ 1 ¼ P L01ðiÞ V ðiÞ1 þ P L11ðiÞ  PGð1iÞ Þ
The derivative of g ðiÞ 1 w.r.t. the same bus i and a different phase p2 –p1 voltage magnitude, as in (14); Pðp Þ
The derivative of g ðiÞ 1 w.r.t. a different bus voltage magnitude j–i, as in (15). X þ  ðp Þ ðp Þ ðp ;p2 Þ V ðiÞ1 V ðiÞ2 Y ði;j1Þ ðp Þ ðp Þ ðp ;p1 Þ V ðiÞ1 V ðiÞ1 Y ði;j1Þ   ðp ;p Þ ðp Þ ðp Þ cos hði;j1Þ 2 þ dðiÞ2  dðiÞ1   ðp ;p Þ cos hði;j1Þ 1 ! j–i   ðp Þ ðp Þ ðp ;p Þ ðp ;p Þ þ V ðiÞ1 V ðiÞ1 Y ði;i1Þ 1 cos hði;i1Þ 1 0 1   X X ðp Þ ðp Þ ðp ;p Þ ðp1 ;p2 Þ ðp2 Þ ðp1 Þ A 1 2 1 2 @ V ðiÞ V ðjÞ Y ði;jÞ cos hði;jÞ þ dðjÞ  dðiÞ  j–i p2 2HðjÞ Pðp Þ Similarly, for the derivative of g ðiÞ 1 w.r.t. the voltage angle, three cases can be identified, as in (16)–(18): Pðp Þ
The derivative of g ðiÞ 1 w.r.t. the same bus i and the same phase p1 voltage angle, as in (16); P ðp  X p2 –p1   ðp ;p Þ ðp Þ ðp Þ cos hði;iÞ1 2 þ dðiÞ2  dðiÞ1 ðp Þ ðp Þ ðp ;p Þ V ðiÞ1 V ðiÞ2 Y ði;iÞ1 2 !    ðp Þ ðp Þ ðp ;p Þ ðp ;p Þ  V ðiÞ1 V ðiÞ1 Y ði;iÞ1 1 cos hði;iÞ1 1 8i 2 I ; p1 2 HðiÞ  P   aði;p Þ 1 ðp Þ ðp Þ ðp Þ ðp Þ 1 J g ðiÞ 1 ; V ðiÞ1 ¼ aði;p1 Þ P L01ðiÞ V ðiÞ1 Þ
The derivative of g ðiÞ 1 w.r.t. the same bus i and a different phase þ p2 –p1 voltage angle, as in (17); Pðp ! X X j2I p2 –p1 Þ
The derivative of g ðiÞ 1 w.r.t. a different bus j–i voltage angle, as þ in (18). X ðp Þ ðp ;p Þ V ðiÞ2 Y ði;j1Þ 2 ðp Þ ðp ;p Þ 2V ðiÞ1 Y ði;j1Þ 1 ð12Þ !   ðp1 ;p2 Þ ðp2 Þ ðp1 Þ cos hði;jÞ þ dðiÞ  dðiÞ !   ðp1 ;p1 Þ cos hði;jÞ j–i 0 1   X X ðp Þ ðp ;p Þ ð Þ ð Þ ð Þ p ;p p p V ðjÞ2 Y ði;j1Þ 2 cos hði;j1Þ 2 þ dðjÞ2  dðiÞ1 A @ Jacobian matrix entries for reactive power constraints j–i p2 2HðjÞ The Jacobian matrix entries that correspond to the reactive X ðp Þ ðp ;p2 Þ power mismatch constraints g ðiÞ 1 in (19) can be derived in a sim-  ilar way since they are derived for the active power. These entries are defined by (20)–(27). 8i 2 I; p1 2 HðiÞ Q ðp Þ p2 –p1 V ðiÞ2 Y ði;i1Þ !   ðp ;p Þ ðp Þ ðp Þ cos hði;i1Þ 2 þ dðiÞ2  dðiÞ1 ð13Þ ð8Þ  aði;p Þ P ðp Þ ðp Þ ðp Þ ðp Þ ðp Þ 1 g ðiÞ 1 ¼ P L01ðiÞ V ðiÞ1 þ P L11ðiÞ  P Gð1iÞ 0 1   X X ðp Þ ðp Þ ðp ;p Þ ðp1 ;p2 Þ ðp2 Þ ð p1 Þ A 1 2 1 2 @ þ V ðiÞ V ðiÞ Y ði;jÞ cos hði;jÞ þ dðiÞ  dðiÞ 0 @ j2I p2 2HðjÞ X X j2I p2 2HðjÞ  P  ðp Þ ðp Þ J g ðiÞ 1 ; PL11ðiÞ ¼ 1 8i 2 I ; p1 2 HðiÞ  P  ðp Þ ðp Þ J g ðiÞ 1 ; PGðiÞ1 ¼ 1 8i 2 I ; p1 2 HðiÞ j2I p2 –p1  1   ðp1 Þ ðp2 Þ ðp1 ;p2 Þ ðp1 ;p2 Þ ðp2 Þ ðp1 Þ A V ðiÞ V ðjÞ Y ði;jÞ cos hði;jÞ þ dðjÞ  dðiÞ 8i 2 I ; p1 2 HðiÞ !   X X  ðp Þ ðp ;p Þ ðp1 ;p1 Þ ðp2 Þ ðp1 Þ 1 1 1 V ðiÞ Y ði;jÞ cos hði;jÞ þ dðiÞ  dðiÞ  P  ðp Þ ðp Þ J g ðiÞ 1 ; V ðiÞ2 ¼ !   X  ðp Þ ðp ;p Þ ðp1 ;p1 Þ ðp2 Þ ðp1 Þ 1 1 1 V ðiÞ Y ði;jÞ cos hði;jÞ þ dðiÞ  dðiÞ p2 –p1 8i 2 I; p1 2 HðiÞ ; p2 –p1 ð9Þ ð10Þ  P    X X  ðp Þ ðp ;p Þ ðp Þ ðp Þ ðp ;p Þ ðp Þ ðp Þ J g ðiÞ 1 ; V ðjÞ2 ¼  V ðiÞ1 Y ði;j1Þ 1 cos hði;j1Þ 1 þ dðjÞ2  dðiÞ1 j–i p2 2HðjÞ 8i 2 I; p1 2 HðiÞ ; j–i ð11Þ ð14Þ ð15Þ 56 M.F. Shaaban et al. / Journal of Advanced Research 20 (2019) 51–60  P  ðp Þ ðp Þ J g ðiÞ 1 ; dðiÞ1 ¼ X X j2I p –p ðp Þ ðp Þ ðp ;p2 Þ V ðiÞ1 V ðiÞ2 Y ði;j1Þ   ðp ;p Þ ðp Þ ðp Þ sin hði;j1Þ 2 þ dðiÞ2  dðiÞ1 ! 0 2 1 1   X X ðp Þ ðp Þ ðp ;p Þ ðp ;p Þ ðp Þ ðp Þ V ðiÞ1 V ðjÞ2 Y ði;j1Þ 2 sin hði;j1Þ 2 þ dðjÞ2  dðiÞ1 A @ j–i p2 2HðjÞ X  ðp Þ ðp Þ ðp ;p2 Þ V ðiÞ1 V ðiÞ2 Y ði;i1Þ p2 –p1   ðp ;p Þ ðp Þ ðp Þ sin hði;i1Þ 2 þ dðiÞ2  dðiÞ1 J j2I p2 –p1 ðp Þ ðp Þ ðp ;p2 Þ V ðiÞ1 V ðiÞ2 Y ði;i1Þ p2 –p1 J P ðp Þ ðp Þ g ðiÞ 1 ; dðjÞ2  ¼ X X  j–i p2 2HðjÞ   ðp ;p Þ ðp Þ ðp Þ sin hði;i1Þ 2 þ dðiÞ2  dðiÞ1 ðp Þ ðp Þ X ðp Þ ðp Þ ðp ;p2 Þ V ðiÞ1 V ðiÞ2 Y ði;i1Þ p2 –p1  j2I p2 –p1 ðp ;p2 Þ  ðp ;p2 Þ sin hði;j1Þ ðp Þ ðp Þ þ dðjÞ2  dðiÞ1  ! 8i 2 I ; p1 2 HðiÞ ð19Þ  Q  ðp Þ ðp Þ J g ðiÞ 1 ; PL11ðiÞ ¼ 1 8i 2 I; p1 2 HðiÞ ð20Þ  Q  ðp Þ ðp Þ J g ðiÞ 1 ; PGðiÞ1 ¼ 1 8i 2 I ; p1 2 HðiÞ ð21Þ  Q   bði;p Þ 1 ðp Þ ðp Þ ðp Þ ðp Þ 1 J g ðiÞ 1 ;V ðiÞ1 ¼ bði;p1 Þ Q L01ðiÞ V ðiÞ1 0 1   X X ðp Þ ðp ;p Þ ð p ;p Þ ð p Þ ð p Þ þ@ V ðjÞ2 Y ði;j1Þ 2 sin hði;j1Þ 2 þ dðjÞ2  dðiÞ1 A j–i p2 2HðjÞ ðp Þ ðp ;p Þ 2V ðiÞ2 Y ði;i1Þ 2 sin X X j2I p2 –p1 X  ðp Þ ðp ;p2 Þ V ðiÞ2 Y ði;j1Þ ðp Þ ðp ;p1 Þ 2V ðiÞ1 Y ði;j1Þ !   ðp1 ;p2 Þ ðp2 Þ ðp1 Þ hði;iÞ þ dðiÞ  dðiÞ !   ðp ;p Þ ðp Þ ðp Þ sin hði;j1Þ 2 þ dðiÞ2  dðiÞ1 !   ðp ;p Þ sin hði;j1Þ 1 8i 2 I;p1 2 HðiÞ j–i ð22Þ  Q  ðp Þ ðp Þ J g ðiÞ 1 ; V ðiÞ2 ¼ X p2 –p1  ðp Þ ðp ;p2 Þ V ðiÞ1 Y ði;i1Þ X X j2I p2 –p1 ðp Þ ðp Þ ðp Þ ðp ;p2 Þ V ðiÞ1 V ðiÞ2 Y ði;j1Þ ðp ;p Þ hði;i1Þ 2 þ ðp Þ dðiÞ2  ðp Þ dðiÞ1   ðp ;p Þ ðp Þ ðp Þ sin hði;i1Þ 2 þ dðiÞ2  dðiÞ1 ðp ;p2 Þ V ðiÞ1 Y ði;j1Þ 8i 2 I ; p1 2 HðiÞ ; p2 –p1 !   ðp ;p Þ ðp Þ ðp Þ cos hði;j1Þ 2 þ dðiÞ2  dðiÞ1  Q  ðp Þ ðp Þ J g ðiÞ 1 ; dðiÞ2 ¼ X p2 –p1  ðp Þ ðp Þ ð25Þ ðp ;p2 Þ V ðiÞ1 V ðiÞ2 Y ði;i1Þ X X ðp Þ j2I p2 –p1 ðp Þ   ðp ;p Þ ðp Þ ðp Þ cos hði;i1Þ 2 þ dðiÞ2  dðiÞ1 ðp ;p2 Þ V ðiÞ1 V ðiÞ2 Y ði;j1Þ J Q ðp Þ ðp Þ g ðiÞ 1 ; dðjÞ2  0 ¼@ X X j–i p2 2HðjÞ ðp Þ ðp Þ ðp ;p Þ V ðiÞ1 V ðjÞ2 Y ði;j1Þ 2 !   ðp ;p Þ ðp Þ ðp Þ cos hði;j1Þ 2 þ dðiÞ2  dðiÞ1 8i 2 I ; p1 2 HðiÞ ; p2 –p1  ! ! ð26Þ cos  ðp ;p Þ hði;j1Þ 2 þ ðp Þ dðjÞ2  ðp Þ dðiÞ1  1 A ð27Þ Results and discussion j–i  X X cos  8i 2 I ; p1 2 HðiÞ ; j–i !   ðp1 ;p2 Þ ðp2 Þ ðp1 Þ sin hði;jÞ þ dðiÞ  dðiÞ    ðp Þ ðp Þ ðp ;p Þ ðp ;p Þ  V ðiÞ1 V ðiÞ1 Y ði;i1Þ 1 sin hði;i1Þ 1 p2 –p1 p2 –p1 ðp Þ ðp Þ ðp ;p Þ V ðiÞ1 V ðiÞ2 Y ði;i1Þ 2 j2I p2 –p1 !   ðp ;p Þ ðp Þ ðp Þ sin hði;i1Þ 2 þ dðiÞ2  dðiÞ1   X ðp Þ ðp Þ ðp ;p Þ ðp ;p Þ  V ðiÞ1 V ðiÞ1 Y ði;j1Þ 1 sin hði;j1Þ 1 þ  ð18Þ ðp Þ ðp Þ ðp ;p Þ V ðiÞ1 V ðiÞ2 Y ði;j1Þ 2 X X 8i 2 I ; p1 2 HðiÞ    ðp Þ ðp Þ ðp ;p Þ ðp ;p Þ þ V ðiÞ1 V ðiÞ1 Y ði;i1Þ 1 sin hði;i1Þ 1 X X j–i p2 2HðjÞ ð17Þ V ðiÞ1 V ðjÞ2 Y ði;j1Þ j–i p2 2HðjÞ ð24Þ 0 1   X X ðp Þ ðp Þ ðp ;p Þ ðp1 ;p2 Þ ðp2 Þ ðp1 Þ A 1 2 1 2 @ ¼ V ðiÞ V ðjÞ Y ði;jÞ cos hði;jÞ þ dðjÞ  dðiÞ !  bði;p Þ ðp Þ ðp Þ ðp Þ ðp Þ 1 ¼ Q L01ðiÞ V ðiÞ1 þ Q L11ðiÞ  Q Gð1iÞ 0 1   X X ðp Þ ðp Þ ðp ;p Þ ðp1 ;p2 Þ ðp2 Þ ðp1 Þ A 1 2 1 2 @ þ V ðiÞ V ðjÞ Y ði;jÞ sin hði;jÞ þ dðjÞ  dðiÞ þ  þ 8i 2 I ; p1 2 HðiÞ ; j–i Q ðp Þ g ð iÞ 1 Q ðp Þ ðp Þ g ðiÞ 1 ; dðiÞ1 ! 8i 2 I ; p1 2 HðiÞ ; p2 –p1   ð16Þ  P    X X ðp Þ ðp Þ ðp ;p Þ ðp Þ ðp Þ ðp ;p Þ ðp Þ ðp Þ J g ðiÞ 1 ; dðiÞ2 ¼  V ðiÞ1 V ðiÞ2 Y ði;j1Þ 2 sin hði;j1Þ 2 þ dðiÞ2  dðiÞ1 X j–i p2 2HðjÞ 8i 2 I ; p1 2 HðiÞ ; j–i ! 8i 2 I ; p1 2 HðiÞ þ 0 1  Q    X X ðp Þ ðp ;p Þ ðp1 Þ ðp2 Þ ð p ;p Þ ð p Þ ð p Þ J g ðiÞ ; V ðjÞ ¼ @ V ðjÞ2 Y ði;j1Þ 2 sin hði;j1Þ 2 þ dðjÞ2  dðiÞ1 A For the purpose of evaluation, the proposed Jacobian was tested on two systems. The first is a simple 4-bus radial system, which is a small part of a real system located in Ontario, Canada. The system has a simple structure, as shown in Fig. 3. The second system is the IEEE 123-bus system [32], which is illustrated in Fig. 4. The load and line data for the 4-bus system and the 123-bus system are shown in Table 1 and in [32], respectively. In both cases, it is assumed that the central controller that hosts the energy management system is located at the substation or main grid connection. For each case, the exact Jacobian Jexact was generated in a MATLABÒ environment using symbolic differentiation of the power flow mismatch equations in (12) and (19). To develop the exact Jacobian, the symbolic differentiation deals with each single constraint as a whole mathematical formula while considering the derivatives w.r.t. all possible variables. On the other hand, the proposed Jacobian J proposed is customized to the structure of the distribution system and develops the derivatives as the sum of individual terms. Mathematically, both Jexact and J proposed are identical in their output. The exact Jacobian development, which is performed off-line, took 1.5 min for the simple 4-bus system and 13.5 h for the 123bus system. On the other hand, the proposed Jacobian can be implemented directly without the need for preprocessing. For large practical systems, the computer will run out of memory trying to perform symbolic differentiation. This was the case when a practical test-case of a 575-bus unbalanced system was tested. !  !   ðp ;p Þ ðp Þ ðp Þ sin hði;j1Þ 2 þ dðiÞ2  dðiÞ1 ð23Þ Fig. 3. 4-Bus distribution test system. 57 M.F. Shaaban et al. / Journal of Advanced Research 20 (2019) 51–60 Fig. 4. IEEE 123-bus distribution test system [32]. Table 1 4-Bus system data in p.u. From (i) To (j) 1 2 2 3 3 4 *ða;aÞ Y ði;jÞ 25:588\  64:25 34:118\  64:25 51:177\  64:25 *ðb;bÞ Y ði;jÞ 26:924\  62:51 35:899\  62:51 53:849\  62:51 *ðc;cÞ Y ði;jÞ 24:382\  65:82 32:510\  65:82. 48:764\  65:82 *ða;bÞ *ðb;aÞ Y ði;jÞ or Y ði;jÞ 9:578\124:31 12:771\124:31 19:157\124:31 *ða;cÞ *ðc;aÞ Y ði;jÞ or Y ði;jÞ 5:496\110:76 7:328\110:76 10:992\110:76 *ðc;bÞ *ðb;cÞ Y ði;jÞ or Y ði;jÞ 7:091\117:80 9:454\117:80 14:182\117:80 P L0ðjÞ ðaÞ 0.309 0.175 0.175 ðbÞ P L0ðjÞ ðcÞ P L0ðjÞ ðaÞ Q L0ðjÞ 0.287 0.175 0.175 0.079 0.175 0.175 0.068 0.088 0.168 ðbÞ Q L0ðjÞ ðcÞ Q L0ðjÞ 0.231 0.088 0.071 0.192 0.088 0.001 compared to the size considered in the previous section (6N bus  18N bus ) is due to the assumption that each bus has only 6 variables instead of 18. The simulations were performed on a 3.6 GHz dual core processor with 16 GB of RAM. The optimization problems for the two case studies were solved using the constrained nonlinear optimization tool based on the interior point method in the MATLABÒ environment. The proposed Jacobian J proposed was constructed for each case. Jexact and Jproposed were compared to an approximated Jacobian Japprox , which was evaluated using a finite difference approximated derivative, as in (28), for low values of Dx: In each case, a battery storage system (BSS) is controlled via an online signal from a central controller that runs the OPF problem, whereby the optimal charging/discharging power is determined and the computational time is evaluated. df f ðx þ DxÞ  f ðxÞ ¼ dx Dx ð28Þ 4-Bus test system For the 4-bus test system, the Jacobian size is 24  24. For the sake of comparison, the Jacobian is evaluated for the initial condi- Moreover, for any topology change due to fault clearing or seasonal reconfiguration, Jexact must be developed for the new topology, which is very time consuming for large systems and will cause a considerable delay in the EMS process. For simplicity, a simple OPF case was considered, in which each ða;b;cÞ bus had six variables: PG ða;b;cÞ andQ G for the slack bus and V ða;b;cÞ and dða;b;cÞ for all other buses. The slack bus is the reference bus where the voltage magnitude and the voltage angles are known. The Jacobian matrix contains the derivative of the power mismatch equations (six equations for each bus), making the Jacobian matrix a size of 6N bus  6N bus . The reduction in the Jacobian matrix size tion x0 , which corresponds to V ða;b;cÞ ¼ 1 and ða;b;cÞ ¼ 0; 2p=3; and2p=3. The root mean square error, RMSE, d which is defined in (29), is used for comparing the results: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 1 X X  RMSE ¼ Eðu;v Þ x 0ðv Þ m  n um v n ð29Þ whereEðu;v Þ is the difference between J exactðu;v Þ and either Japproxðu;v Þ or Jproposedðu;v Þ . As shown in Fig. 5, as Dx decreases, the value of the RMSE for the approximate Jacobian is further reduced. However, reducing Dx is 58 M.F. Shaaban et al. / Journal of Advanced Research 20 (2019) 51–60 Fig. 5. The RMSE of the approximate Jacobian versus Dx for the 4-bus test system. problematic because it approaches the precision limit of the software and hardware used for the computations. In this case, as Dx was reduced below 109, the value of the RMSE increased to 1.396  109 at x0 . On the other hand, for Jproposed , the error was zero, which shows that Jproposed is identical to Jexact . The three Jacobian matrices were evaluated at 1000 different values for x0 w.r.t. the computation time, t. The computational time statistics are shown in Table 2. As indicated, the proposed Jacobian computational time was much less than that for the approximate method but very close to the exact method. The mean computation time for the approximate Jacobian was 7.4 ms compared to 0.642 ms for the proposed Jacobian, which is almost 11 times faster. This discrepancy is attributable to two factors: (1) most of the derivatives are zeros, which is easily identified in Jproposed but must be evaluated for Japprox ; (2) for each term in the Jacobian matrix, Japprox must evaluate a function twice as in (28), while Jproposed evaluates only one function. A comparison of the computational time of Jproposed to that of Jexact reveals a slightly faster performance for the proposed formulation, as identified by a 2.4% reduction in the mean computational time. This decrease is due to the complexity of the form of Jexact , where symbolic differentiation is considered in comparison to the compact form of Jproposed . As the size and complexity of the system increases, this difference is expected to be significant. To test the proposed Jacobian for solving an OPF problem, a BSS is presumed to be located at bus 4. To simplify the analysis, the BSS optimal schedule (to charge and discharge) will be determined by the OPF run to minimize the overall system losses. The BSS is also assumed to operate at a unity power factor. Although the test case study is simple, the proposed Jacobian formulation is applicable (and scalable) to any OPF case study, including DER optimal management in microgrids and distribution networks, transactive energy markets modeling, simulation and analysis. To determine the optimal power amount injected from the BSS Pbat , a simple test was performed before applying the OPF, in which the injected power from the BSS was varied and the system losses were evaluated. The test outcomes showed that the system loss for the base case (Pbat ¼ 0) was 0.0229 p.u., which corresponded to 1.33% of the total system demand. On the other hand, the minimum loss was 0.0082 p.u., or 0.48%, which occurred at Pbat ¼ 1:24 p.u., as shown in Fig. 6. To test the proposed Jacobian computation method, an OPF problem was formulated as follows: min PTot loss ¼ P bat X ðp1Þ ð30Þ P loss ph1 is subject to ðp Þ PGð1iÞ þ ðp Þ P bat1 ðiÞ  ðp Þ PLði1Þ 0 1   X X ðp Þ ðp Þ ðp ;p Þ ð p ;p Þ ð p Þ ð p Þ 1 2 1 2 1 2 2 1 ¼@ V ðiÞ V ðiÞ Y ði;jÞ cos hði;jÞ þ dðiÞ  dðiÞ A 0 @ j2I p2 2HðjÞ X X j2I p2 2HðjÞ ðp Þ ðp Þ ðp ;p2 Þ V ðiÞ1 V ðjÞ2 Y ði;j1Þ 8i 2 I ; p1 2 HðiÞ 1   ðp ;p Þ ðp Þ ðp Þ cos hði;j1Þ 2 þ dðjÞ2  dðiÞ1 A ð31Þ Table 2 Computation time (ms) for evaluating the Jacobian matrix for the 4-bus test system. Mean Minimum Maximum Jexact Jproposed Japprox 0.6579 0.6401 1.0095 0.6423 0.6250 0.9900 7.4060 7.278 8.4856 Fig. 6. Total system losses for the 4-bus system versus the BSS output power. M.F. Shaaban et al. / Journal of Advanced Research 20 (2019) 51–60 ðp Þ ðp Þ Q Gð1iÞ  Q Lði1Þ ¼  X X X þ ðp Þ p2 2HðjÞ j2I X p2 2HðjÞ j2I ðp Þ ðp ;p2 Þ V ðiÞ1 V ðiÞ2 Y ði;j1Þ   ðp ;p Þ ðp Þ ðp Þ sin hði;j1Þ 2 þ dðiÞ2  dðiÞ1   ð Þ ð Þ ðp1 ;p2 Þ ð Þ ð Þ ð Þ V ðpiÞ1 V ðpjÞ2 Y ði;j sin hðpi;j1Þ;p2 þ dðpjÞ2  dðpiÞ1 Þ 8i 2 I ; p1 2 HðiÞ ðp1Þ Ploss ¼ X ðp Þ ðp Þ 59 ð32Þ ðp Þ PGð1iÞ þ P bat1 ðiÞ  PLði1Þ 8p1 ð33Þ i ðp Þ V min  V ðiÞ1  V max ðp Þ ðaÞ 8i 2 I ; p1 2 HðiÞ ðbÞ V ðiÞ1 ¼ V slack ; dðiÞ ¼ 0; dðiÞ ¼  ðp Þ ðp Þ PGð1iÞ ; Q Gð1iÞ ¼ 0 ðp Þ Pbat1 ðiÞ ¼ ð34Þ 2p ðcÞ 2p ;d ¼ 3 ð iÞ 3 8i 2 Islack 8i R Islack ð36Þ PbatðiÞ =3 8i 2 Ibat 0 ð35Þ ð37Þ otherwise where V min andV max are the minimum and maximum voltage limits, respectively, which are set to 0.95 and 1.05, respectively;V slack is the slack bus voltage, which is set to 1.05 p.u.; Islack ¼ f1g and Ibat ¼ f4g are the subsets of the slack bus and the BSS bus, respectively. The nonlinear programming (NLP) OPF problem, defined in (30)–(37), was solved using an interior-point algorithm in the MATLAB environment with default settings. The problem was first solved using finite differences. The optimal solution was P Tot loss ¼ 0:0082 p.u. at P bat ð4Þ ¼ 1:2366 p.u., which was very close to the values indicated in Fig. 6. The computational time was 100.81 ms. On the other hand, when the gradient of the objective function and the Jacobian of the constraints were provided, the optimal solution was PTot loss ¼ 0:0082 p.u. at P batð4Þ ¼ 1:2366 p.u., which was identical to the values derived from using finite differences. The process took 23.76 ms, which represents a 76.43% reduction in the computational time. 123-Bus test system The 123-bus system is assumed to be moderately sized for a practical distribution system, which may contain more than 500 buses. Closing the switches and removing zero impedance lines reduces the 123-bus system to 119 buses, in which the size of the Jacobian matrix is 714  714, which is equivalent to 509,796 entries in the matrix. The time required for evaluating the Jacobian matrices was recorded for 100 different x0 values. As shown in Table 3, the mean time for evaluating Japprox was 24% less than the mean time required for evaluating Jexact with a zero RMSE value. On the other hand, forDx = 109, the RMSE for Japprox was 4.806  106 at x0 ; this differs significantly from the RMSE for the 4-bus system, which was 1.396  109. The mean computational time for Japprox was 4.72 s, which was 78 times higher than that for Jproposed . For a moderately sized system, evaluating Japprox for a single iteration may consume more than 4 s, depending on the hardware and software used to solve for OPF. Thus, for a practical distribu- Fig. 7. The total system losses for the 123-bus system versus the BSS output power. tion system with hundreds of buses, evaluating Japprox in real time is impractical, and the error increases dramatically as the system becomes larger. In addition, developing Jexact is very time consuming and may lead to memory limitation errors, as explained previously. These considerations provide strong evidence of the need for a closed form Jacobian matrix, which could be easily coded and applied to solving real-time optimization problems in distribution systems, such as DER optimal management, transactive energy market clearing and the optimal pricing of DERs. It is assumed that a BSS unit operating at a unity power factor is located at bus 49, which is chosen arbitrarily. Changing the output power of the BSS helps us in evaluating and minimizing the distribution system losses, as shown in Fig. 7. As indicated in the figure, the loss for the base case without the BSS was 95.94 kW at point a, which was very close to the value of the system loss reported in [32] (95.611 kW). As the injected power from the BESS increases, the system losses decrease, where the minimum loss was 64.05 kW at Pbat = 1.667 MW, which is shown as point b in Fig. 7. The decrease in the system losses is attributed to the reduction in the current flowing from the grid to supply the loads, which are supplied partially by the BESS. For losses greater than 1.667 MW, increasing the injected power from the BESS causes an increase in the losses due to reverse power, where the excess current that flows from the BESS to the grid direction causes higher losses in the lines. To test the proposed Jacobian matrix computation method, the same NLP problem expressed in (30)–(37) was used, with the BSS connected to bus 49. The minimum system loss was 61.045 kW at Pbat = 1.6151 MW. Although the finite difference and the proposed Jacobian reached the same optimal solution, the computational times were 0.82158 and 25.57 s with and without applying the Jacobian method, respectively. The case study results clearly showed that applying the proposed Jacobian method could substantially reduce the OPF computational time. Moreover, as the size of the system increases, the reduction becomes even greater, which emphasizes the necessity for providing the Jacobian of the constraints in real-time DER management or transactive energy market clearing; this may require solving a large OPF problem. Conclusions Table 3 Computation time (ms) for evaluating the Jacobian matrix for the 123-Bus test system. Mean Minimum Maximum Jexact Jproposed Japprox 80.22 78.68 82.90 60.35 59.95 63.50 4718.97 4695.60 4800.86 The work presented in this paper tackles a core element of any modern EMS, which is the optimization unit. The proposed approach targets the computational process for the optimization unit for real-time unbalanced SDN operations (e.g., smart grid applications). The proposed approach could help reduce the complexity and computational time associated with the SDN applications, which usually involve real-time optimal scheduling of the 60 M.F. Shaaban et al. / Journal of Advanced Research 20 (2019) 51–60 system assets (including DERs and demand response/DR). The paper includes a detailed formulation of the proposed generalized Jacobian matrix, which can be tailored according to the real-time measurements and available system data. The formulation presented for unbalanced SDNs is based on three sets of derivatives for the power mismatch constraints: (1) the derivative w.r.t. the same bus and the same phase voltage magnitude, (2) the derivative w.r.t. the same bus and a different phase voltage magnitude, and (3) the derivative w.r.t. different bus voltage magnitudes. The proposed Jacobian matrix formulation can be applied to any distribution network after adjustments are made depending on the available equipment, such as voltage regulators and capacitor banks. The simulation results demonstrated the effectiveness and scalability of the proposed formulation for evaluating the values of the Jacobian matrix entries in a timely manner and with zero RMSE error compared to the exact and approximated Jacobian evaluation methods. The results also provided enough evidence to support the need for using the proposed Jacobian matrix formulation in realtime OPF problem solving for smart grid applications, such as DER/DR optimal management and transactive energy market clearing. Conflict of interest The authors have declared no conflict of interest. Compliance with Ethics Requirements This article does not contain any studies with human or animal subjects. References [1] Abdelaziz MM, Farag HE, El-Saadany E, Mohamed Y-R. A globally convergent trust-region method for power flow studies in active distribution systems. Paper presented at: 2012 IEEE pow ener soc ge, 2012. [2] Lightner EM, Widergren SE. An orderly transition to a transformed electricity system. IEEE Trans Smart Grid 2010;1(1):3–10. [3] Wu W, Zhang B. A three-phase power flow algorithm for distribution system power flow based on loop-analysis method. Int J Elect Power 2008;30(1):8–15. [4] Berg Jr R, Hawkins E, Pleines W. Mechanized calculation of unbalanced load flow on radial distribution circuits. IEEE Trans Power Apparat Syst 1967;1 (4):415–21. [5] Zhang F, Chen CS. A modified Newton method for radial distribution system power flow analysis. IEEE Trans Power Syst 1997;12(1):389–97. [6] Van Amerongen RA. A general-purpose version of the fast decoupled load flow. IEEE Trans Power Syst 1989;4(2):760–70. [7] Shirmohammadi D, Hong H, Semlyen A, Luo G. A compensation-based power flow method for weakly meshed distribution and transmission networks. IEEE Trans Power Syst 1988;3(2):753–62. [8] CYME user guide [Internet]. Distribution system analysis [cited 2018 October 10]. Available from: http://www.cyme.com/software/cymdist/. [9] Jasinski R, Sablerolle W, Amory M. ETAP: Scale prediction and control for the heron cluster. Paper presented at: SPE annual technical conference and exhibition, 1997. [10] Momoh JA. Electric power system applications of optimization. 2nd ed. CRC Press; 2008. [11] Dall’Anese E, Zhu H, Giannakis G. Distributed optimal power flow for smart microgrids. IEEE Trans Smart Grid 2013;4(3):1464–75. [12] Lavaei J, Low SH. Zero duality gap in optimal power flow problem. IEEE Trans Power Syst 2012;27(1):92–107. [13] Sherali HD, Adams WP. A hierarchy of relaxations and convex hull characterizations for mixed-integer zero—one programming problems. Discrete Appl Math 1994;52(1):83–106. [14] Ryoo HS, Sahinidis NV. A branch-and-reduce approach to global optimization. J Global Optim 1996;8(2):107–38. [15] Bhattarai BP, Lévesque M, Maier M, Bak-Jensen B, Radhakrishna Pillai J. Optimizing electric vehicle coordination over a heterogeneous mesh network in a scaled-down smart grid testbed. IEEE Trans Smart Grid 2015;6(2):784–94. [16] Ozturk Y, Senthilkumar D, Kumar S, Lee G. An intelligent home energy management system to improve demand response. IEEE Trans Smart Grid 2013;4(2):694–701. [17] Benetti G, Delfanti M, Facchinetti T, Falabretti D, Merlo M. Real-time modeling and control of electric vehicles charging processes. IEEE Trans Smart Grid 2015;6(3):1375–85. [18] Azzouz MA, Shaaban MF, El-Saadany EF. Real-time optimal voltage regulation for distribution networks incorporating high penetration of PEVs. IEEE Trans Power Syst 2015;30(6):3234–45. [19] Shaaban MF, Ismail M, El-Saadany EF, Zhuang W. Real-time PEV charging/ discharging coordination in smart distribution systems. IEEE Trans Smart Grid 2014;5(4):1797–807. [20] Farag HE, El-Saadany EF. A novel cooperative protocol for distributed voltage control in active distribution systems. IEEE Trans Power Syst 2013;28 (2):1645–56. [21] Shi W, Li N, Xie X, Chu C-C, Gadh R. Optimal residential demand response in distribution networks. IEEE J Sel Area Commun 2014;32(7):1441–50. [22] Soroudi A, Siano P, Keane A. Optimal DR and ESS scheduling for distribution losses payments minimization under electricity price uncertainty. IEEE Trans Smart Grid 2016;7(1):261–72. [23] Franco JF, Rider MJ, Romero R. A mixed-integer linear programming model for the electric vehicle charging coordination problem in unbalanced electrical distribution systems. IEEE Trans Smart Grid 2015;6(5):2200–10. [24] Sarker MR, Ortega-Vazquez MA, Kirschen DS. Optimal coordination and scheduling of demand response via monetary incentives. IEEE Trans Smart Grid 2015;6(3):1341–52. [25] Maffei A, Srinivasan S, Castillejo P, Martínez JF, Iannelli L, Bjerkan E, et al. A semantic middleware supported receding horizon optimal power flow in energy grids. IEEE Trans Ind Inf 2018;14(1):35–46. [26] Sharma I, Canizares C, Bhattacharya K. Smart charging of PEVs penetrating into residential distribution systems. IEEE Trans Smart Grid 2014;5(3):1196–209. [27] Sharma I, Bhattacharya K, Canizares C. Smart distribution system operations with price-responsive and controllable loads. IEEE Trans Smart Grid 2015;6 (2):795–807. [28] Zhang W, Xu Y, Dong Z, Wong P. Robust security-constrained optimal power flow using multiple microgrids for corrective control under uncertainty. IEEE Trans Ind Inf 2017;13(4):1704–13. [29] Siano P, Cecati C, Yu H, Kolbusz J. Real time operation of smart grids via FCN networks and optimal power flow. IEEE Trans Ind Inf 2012;8(4):944–52. [30] Mostafa HA, El-Shatshat R, Salama M. Multi-objective optimization for the operation of an electric distribution system with a large number of single phase solar generators. IEEE Trans Smart Grid 2013;4(2):1038–47. [31] Mathworks [Internet]. Optimization toolbox: user’s guide (R2016a) [cited 2018 October 10]. Available from: uk.mathworks.com/help/pdf_doc/optim/ optim_tb.pdf. [32] IEEE 123 node test feeder [Internet]. Available from: http://ewh.ieee.org /soc/ pes/dsacom/testfeeders/.