IEEE TRANSACTIONS ON SMART GRID, VOL. 6, NO. 6, NOVEMBER 2015 2825 Optimal Operation of Distribution Networks Considering Energy Storage Devices Leonardo H. Macedo, Student Member, IEEE, John F. Franco, Member, IEEE, Marcos J. Rider, Member, IEEE, and Rubén Romero, Senior Member, IEEE Abstract—This paper presents a mixed-integer second-order cone programing (MISOCP) model to solve the optimal opera- tion problem of radial distribution networks (DNs) with energy storage. The control variables are the active and reactive gen- erated power of dispatchable distributed generators (DGs), the number of switchable capacitor bank units in operation, the tap position of the voltage regulators and on-load tap-changers, and the operation state of the energy storage devices. The objective is to minimize the total cost of energy purchased from the dis- tribution substation and the dispatchable DGs. The steady-state operation of the DN is modeled using linear and second-order cone programing. The use of an MISOCP model guarantees convergence to optimality using existing optimization software. A mixed-integer linear programing (MILP) formulation for the original model is also presented in order to show the accuracy of the proposed MISOCP model. An 11-node test system and a 42-node real system were used to demonstrate the effectiveness of the proposed MISOCP and MILP models. Index Terms—Distributed generation, energy storage, mixed-integer linear programing (MILP), mixed-integer second-order cone programing (MISOCP), optimal operation of radial distribution networks, smart grid. ACRONYMS The abbreviations of common terms used in this paper are presented below. CB Fixed capacitor bank. DG Distributed generator. DN Distribution network. DSS Distribution substation. ESD Energy storage device. MILP Mixed-integer linear programing. MINLP Mixed-integer nonlinear programing. MISOCP Mixed-integer second-order cone programing. OLTC On-load tap-changer. OODN Optimal operation of DNs. OPF Optimal power flow. RS Renewable source. Manuscript received May 30, 2014; revised November 3, 2014 and February 9, 2015; accepted March 30, 2015. Date of publication April 27, 2015; date of current version October 17, 2015. This work was supported in part by the Conselho Nacional de Desenvolvimento Científico e Tecnológico under Grant 304234/2013-3, and in part by the Fundação de Amparo à Pesquisa do Estado de São Paulo under Grant 2012/23454-4. The authors are with the Faculdade de Engenharia de Ilha Solteira, Departamento de Engenharia Elétrica, Universidade Estadual Paulista, Ilha Solteira 15385-000, São Paulo, Brazil (e-mail: mjrider@dee.feis.unesp.br). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TSG.2015.2419134 SCB Switchable capacitor bank. VR Voltage regulator. NOMENCLATURE The notation used throughout this paper is reproduced below for quick reference. Sets �n Set of nodes. �b Set of branches. �d Set of load levels. �S Set of DSS nodes. �dg Set of dispatchable DGs. �cb Set of CBs. �scb Set of SCBs. �vr Set of VRs. �rs Set of nondispatchable RSs. �sd Set of ESDs. Function f (y, y, �) Approximation of the square value of y. Constants cdg i,d Purchase price of energy from dispatchable DG at node i at load level d (US$/kWh). cS d Purchase price of energy from DSS at load level d (US$/kWh). �t Time duration of each load level (h). PD i,d Active power load at node i at load level d (kW). QD i,d Reactive power load at node i at load level d (kvar). Rij Resistance of branch ij (m�). Xij Reactance of branch ij (m�). Zij Impedance of branch ij (m�). Iij Maximum current magnitude of branch ij (A). V Maximum voltage magnitude (kV). V Minimum voltage magnitude (kV). Vnom Nominal voltage magnitude (kV). S S i Maximum apparent power limit of substation at node i (kVA). pf dg i Lower limit of the capacitive power factor for the DG at node i. pf dg i Lower limit of the inductive power factor for the DG at node i. 1949-3053 c© 2015 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. 2826 IEEE TRANSACTIONS ON SMART GRID, VOL. 6, NO. 6, NOVEMBER 2015 S dg i Maximum apparent power limit of the DG at node i (kVA). qscb i Reactive power capacity of each module of the SCB at node i (kvar). Qcb i Reactive power injection of the CB at node i (kvar). nscb i,d Maximum integer number of SCB units at node i. � scb i Maximum variation of capacitor bank units at node i operating in the considered time period. rvr ij Regulation percentage of the VR connected to branch ij. l vr ij Maximum number of steps of the VR in branch ij. � vr ij Maximum variation of steps of the tap of the VR in branch ij in the considered time period. Prs i,d Active power generation of nondispatchable RS at node i at load level d (kW). Psd+ i Minimum power injection capacity of the storage device at node i (kW). P sd+ i,d Maximum power injection capacity of the storage device at node i (kW). Psd− i Minimum power extraction capacity of the stor- age device at node i (kW). P sd− i,d Maximum power extraction capacity of the stor- age device at node i (kW). Esd i Minimum limit of energy storage of the storage device at node i (kWh). E sd i Maximum limit of energy storage of the storage device at node i (kWh). ηsd+ i Power injection efficiency of the storage device at node i. ηsd− i Power extraction efficiency of the storage device at node i. βsd i Self-discharge rate of the storage device at node i. � sd i Maximum number of changes in the storage device operation at node i in the considered time period. � Number of discretizations used in function f . σy,γ Slope of the γ th discretization of y. y Maximum value of y. Continuous Variables Vi,d Voltage magnitude at node i at load level d (kV). Vsqr i,d Square of Vi,d. PS i,d Active power supplied by substation at node i at load level d (kW). QS i,d Reactive power supplied by substation at node i at load level d (kvar). Iij,d Current magnitude of branch ij at load level d (A). Isqr ij,d Square of Iij,d. Pij,d Active power flow of branch ij at load level d (kW). Qij,d Reactive power flow of branch ij at load level d (kvar). Pvr ij,d Active power flow of the VR between nodes i and j at load level d (kW). Qvr ij,d Reactive power flow of the VR between nodes i and j at load level d (kvar). Pdg i,d Active power generation of the DG at node i at load level d (kW). Qdg i,d Reactive power generation of the DG at node i at load level d (kvar). Qscb i,d Reactive power injection of the SCB at node i at load level d (kvar). tvr ij,d Tap of the VR in branch ij at load level d. Vc i,j,d,k Correction variable used in the linearization of the VR model. Psd+ i,d Active power production (injection) of the storage device at node i at load level d (kW). Psd− i,d Active power storage (extraction) of the storage device at node i at load level d (kW). Esd i,d State of charge of the storage device at node i at load level d (kWh). n+ i,d Auxiliary variable that indicates a positive change in the number of SCB units connected at node i. n− i,d Auxiliary variable that indicates a negative change in the number of SCB units connected at node i. l+ij,d Auxiliary variable that indicates a positive vari- ation in the steps of the tap of the VR in branch ij. l−ij,d Auxiliary variable that indicates a negative vari- ation in the steps of the tap of the VR in branch ij. e+ i,d Auxiliary variable that indicates a positive change in the operation state of the storage device at node i. e− i,d Auxiliary variable that indicates a negative change in the operation state of the storage device at node i. δy,γ Value of the γ th auxiliary variable used in the discretization of y. y Argument of function f . y+ Auxiliary variable in the calculation of |y|. y− Auxiliary variable in the calculation of |y|. Binary Variables esd i,d Operation state of the storage device at node i. bvr ij,d,k Binary variable that defines the operation region for the tap of the VR in branch ij. Integer Variables nscb i,d Number of SCB units in operation at node i at load level d. lvr ij,d Number of tap steps of the VR connected to branch ij at load level d. I. INTRODUCTION THE PROBLEM of the OODN considering ESDs aims at minimizing an objective function by determining the optimal values for a set of decision variables. The decision variables for the OODN problem are commonly the tap MACEDO et al.: OODNs CONSIDERING ESDs 2827 position in the OLTCs and VRs, the number of units in operation in each SCB, the active and reactive power injections by dispatchable DGs, and the operation state (injection or extraction of power from the network) of the ESDs. The objec- tive is to minimize the total purchase cost of energy from the DSS and the dispatchable DGs over a predetermined period of time [1]. In the specialized literature, works dealing with the OODN problem have only considered some of the devices previ- ously cited [2]–[7]; moreover, these works have not addressed energy storage. Reference [2] presented a genetic algo- rithm (GA) to obtain the optimal control of voltage magnitude in DNs with VRs, SCBs, DGs, and static var compensators. Reference [3] proposed an evolutionary algorithm to solve the OODN problem, accounting for SCBs, OLTCs in the sub- station, and the reconfiguration of the system. In [4], VRs, OLTCs, and photovoltaic generation were considered with the objective of minimizing the number of tap operations of the VRs, employing an optimal coordination strategy. A model for the OODN problem was proposed in [5] with the goal of minimizing power losses; in addition, this model included photovoltaic generation and reconfiguration of the system. Reference [6] presented an MILP model for the OODN prob- lem considering VRs, SCBs, and DGs. Reference [7] proposed a Chu-Beasley’s GA and a simplified algorithm to solve the OODN problem, also taking VRs, SCBs, and DGs into account. In [8]–[12], the optimal dispatch of ESDs was examined, but with respect to simple DN topologies. Reference [13] listed many benefits of using ESDs in DNs. Despite all of these benefits, energy storage is still expensive for high penetra- tions. However, with the development of technology, the price of ESDs has been decreasing and cost/benefit analyses have shown that the application of these devices is justified in some cases [13]. There have also been works tackling RSs and ESDs simul- taneously, along with the general network topologies in the OODN problem [14]–[18]. However, due to the complexity of the model, VRs, SCBs, and DGs have not been dealt with. For example, Levron et al. [14] presented a methodology for solving the OPF in microgrids with energy storage and RSs. The objective was to minimize the energy purchase cost at the point of common coupling. The proposed method used a load flow and dynamic search algorithm to attempt in finding an optimal solution for the OODN problem. The method worked well for systems with a small number of ESDs, but could not solve problems with more than five ESDs. Furthermore, the proposed methodology did not guarantee convergence to a global optimum. Reference [15] explicated a solution for active–reactive OPF in DNs with RSs and ESDs. Although both injected and extracted power were flexible, i.e., they could assume differ- ent values, the charging and discharging periods were fixed during the days—a consideration that led to low-quality solu- tions. Reference [16] made an improvement to the model in [15] by incorporating the flexible operation of the ESDs, i.e., both injected and extracted power by ESDs, and by allowing the duration of the injection and extraction of Fig. 1. Illustrative DN with three nodes. energy from the network to vary. The methodologies pre- sented in [15] and [16] are applicable to meshed networks. The solution technique proposed in [16] assumes only one charge/discharge cycle per day. Reference [17] proposed a dynamic OPF to solve the OODN considering only RSs, energy storage, and active network man- agement. Reference [18] presented several models with which to formulate ideal and generic storage devices. Most of these works [2]–[17] used metaheuristic or heuristic techniques to solve the OODN problem. Although these techniques are robust, flexible, and achieve good results, they present many problems, such as a high computational demand, the adjust- ment and tuning of parameters, and the definition of a stop criterion. In addition, they cannot guarantee convergence to the global optimum [19]. Also, [14]–[18] did not consider the presence of VRs, SCBs, and dispatchable DGs in the DN. This paper presents an MISOCP model for solving the OODN problem, guaranteeing convergence to the global opti- mum of the original MINLP model. The proposed model con- siders dispatchable DGs, SCBs, VRs, OLTCs, RSs, and ESDs simultaneously and is applicable to radial networks in gen- eral. In addition, the proposed formulation for ESDs permits flexible operation and modification to the number of charge cycles in the period of analysis. The use of an MISOCP model guarantees optimality by using existing classical optimization tools. An approximated MILP model is also presented in order to show the accuracy of the proposed MISOCP model. The main contribution of this paper is a novel MISOCP model for solving the OODN problem considering ESDs, which has the following benefits: 1) a flexible, realistic, and precise model; 2) efficient computational behavior with conventional MISOCP solvers; and 3) convergence to optimality that is guaranteed by using classical optimization techniques. II. MODELS FOR THE OODN PROBLEM A. Steady-State Operation of Distribution Networks The following assumptions are made in order to represent the steady-state operation of a radial DN. 1) The load is represented as constant active and reactive power. 2) The active and reactive power losses of branch ij are concentrated at node i. 3) The DN is balanced and represented by a single-phase equivalent. These three considerations are shown in Fig. 1 for each load level d. The terms RijI2 ij,d and XijI2 ij,d represent the active 2828 IEEE TRANSACTIONS ON SMART GRID, VOL. 6, NO. 6, NOVEMBER 2015 Fig. 2. Dispatchable DG. (a) Scheme and (b) capability curve. and reactive power losses of branch ij, respectively. The equa- tions representing the steady-state operation of a radial DN are shown in (1)–(4) [20]. These equations are frequently used in the load flow sweep method [21], [22] ∑ ki∈�b Pki,d − ∑ ij∈�b ( Pij,d + RijI 2 ij,d ) + PS i,d = PD i,d ∀i ∈ �n,∀d ∈ �d (1)∑ ki∈�b Qki,d − ∑ ij∈�b ( Qij,d + XijI 2 ij,d ) + QS i,d = QD i,d ∀i ∈ �n,∀d ∈ �d (2) V2 i,d − V2 j,d = 2 ( RijPij,d + XijQij,d )+ Z2 ijI 2 ij,d ∀ij ∈ �b,∀d ∈ �d (3) I2 ij,dV2 j,d = P2 ij,d + Q2 ij,d ∀ij ∈ �b,∀d ∈ �d. (4) The active and reactive power balance equations are repre- sented by (1) and (2), respectively. The magnitude of voltage drop is calculated in terms of the active and reactive power flows, the current magnitude, and the electrical parameters of branch ij, as expressed in (3), while (4) is the square current magnitude calculation. B. Modeling of the Equipment 1) Dispatchable Distributed Generators: In this paper, it is assumed that dispatchable DGs (e.g., biomass genera- tors) are synchronous machines that are sources of both active and reactive power. It is also assumed that DGs operate with a restricted power factor and free terminal volt- age. Fig. 2(a) shows the scheme of a dispatchable DG and Fig. 2(b) shows its capability curve. Equations (5)–(7) model the operation of dispatchable DGs ( Pdg i,d )2 + ( Qdg i,d )2 ≤ ( S̄dg i )2 ∀i ∈ �dg,∀d ∈ �d (5) Pdg i,d ≥ 0 ∀i ∈ �dg,∀d ∈ �d (6) −Pdg i,d tan ( cos−1 ( pf dg i )) ≤ Qdg i,d ≤ Pdg i,d tan ( cos−1 ( pf dg i )) ∀i ∈ �dg,∀d ∈ �d. (7) Constraint (5) is the power generation capacity, while (6) establishes that the active power generated by the dispatch- able DG connected at node i is always positive. Constraint (7) limits the reactive power generation to within capacitive and inductive minimum power factors. 2) Capacitor Banks: There are two types of capacitor banks: 1) CBs and 2) SCBs. CBs are formed by units that are always connected to the DN, whereas SCBs have units that may be fully or partially connected to the DN. Fig. 3 shows Fig. 3. SCB. Fig. 4. VR. a scheme of an SCB. Expressions (8)–(11) model the opera- tion of SCBs Qscb i,d = nscb i,d qscb i ∀i ∈ �scb,∀d ∈ �d (8) ∑ d∈�d ∣∣∣nscb i,d − nscb i,d−1 ∣∣∣ ≤ �̄scb i ∀i ∈ �scb (9) 0 ≤ nscb i,d ≤ n̄scb i,d ∀i ∈ �scb,∀d ∈ �d (10) nscb i,d integer ∀i ∈ �scb,∀d ∈ �d. (11) Equation (8) imposes that the total reactive power injected by the SCB connected at node i at load level d is equal to the number of units connected at the node times the reactive power of each unit. Constraint (9) limits the number of changes in the operation of switchable units during the period of analysis, where for d = 0, nscb i,0 is the initial value of nscb i,d ; (10) limits the number of units connected at a node. The integer condition of the number of units connected is established by (11). 3) Voltage Regulators and OLTCs: VRs and OLTCs aim at improving voltage levels in DNs. They have a mechanism that allows the adjustment of the tap, thus adjusting the value of regulation. In this paper, the switching mechanisms of VRs and OLTCs are remotely controlled. Fig. 4 shows a VR connected between nodes i and j. Equations (12)–(16) show the model of the VR Vj,d = tvr ij,dVi,d ∀ij ∈ �vr,∀d ∈ �d (12) tvr ij,d = 1 + rvr ij lvr ij,d l̄vr ij ∀ij ∈ �vr,∀d ∈ �d (13) ∑ d∈�d ∣∣∣lvr ij,d − lvr ij,d−1 ∣∣∣ ≤ �̄vr ij ∀ij ∈ �vr (14) −l̄vr ij ≤ lvr ij,d ≤ l̄vr ij ∀ij ∈ �vr,∀d ∈ �d (15) lvr ij,d integer ∀ij ∈ �vr,∀d ∈ �d. (16) Equation (12) calculates the voltage at node j as a func- tion of the nonregulated voltage at node i and the value of MACEDO et al.: OODNs CONSIDERING ESDs 2829 Fig. 5. ESD. the tap. From (13), it can be shown that the value of the tap varies in the range of regulation [1 − rvr ij , 1 + rvr ij ] in intervals of rvr ij /l̄vr ij . Constraint (14) defines the maximum variation of the VR/OLTC tap steps during the period of analysis, where for d = 0, lvr ij,0 is the initial value of lvr ij,d. Expression (15) represents the interval of variation of the discrete steps of the tap. Expression (16) determines that the steps of the VRs and OLTCs must be integers [6]. 4) Intermittent Renewable Sources: Intermittent RSs of energy (wind turbines and photovoltaic panels) are considered in this paper as nondispatchable sources operating with unity power factor, i.e., they can only produce active power. Because they are nondispatchable, they are represented as injections of active power in the node at which they are installed. Thus, the value of the injected power is a parameter, and not a variable, that depends on the weather forecast, as shown in [14]. 5) Energy Storage Devices: It is appropriate to store energy in storage devices in two situations: 1) when energy from nondispatchable RSs is in excess in the system and 2) when the purchase costs of energy from DSS or dispatchable DGs are low. The stored energy can then be consumed during peak hours, when energy is expensive and load is high [10]. Fig. 5 shows a diagram of the model of an ESD connected at node i. Equations (17)–(22) model the operation of the ESDs Psd+ i esd i,d ≤ Psd+ i,d ≤ P̄sd+ i,d esd i,d ∀i ∈ �sd,∀d ∈ �d (17) Psd− i ( 1 − esd i,d ) ≤ Psd− i,d ≤ P̄sd− i,d ( 1 − esd i,d ) ∀i ∈ �sd,∀d ∈ �d (18) Esd i,d = Esd i,d−1 + ηsd− i �tPsd− i,d − 1 ηsd+ i �tPsd+ i,d − βsd i �tEsd i,d ∀i ∈ �sd,∀d ∈ �d (19) Esd i ≤ Esd i,d ≤ Ēsd i ∀i ∈ �sd,∀d ∈ �d (20) ∑ d∈�d ∣∣∣esd i,d − esd i,d−1 ∣∣∣ ≤ �̄sd ij ∀i ∈ �sd (21) esd i,d ∈ {0, 1} ∀i ∈ �sd,∀d ∈ �d (22) where if esd i,d = 0, the ESD is extracting power from the DN and if esd i,d = 1, the ESD is injecting power into the DN. Constraint (17) limits power injection and (18) restricts power extraction from the ESDs. Equation (19) determines that the stored energy (charge state) of the ESD connected at node i at load level d depends on the previous state of charge, the injected and extracted power for the time interval (duration of load level) multiplied by their respective efficiencies, and the self-discharge rate of the ESD. For d = 0, Esd i,0 is the ini- tial value of Esd i,d. Expression (20) represents the extent of the charge state, i.e., the maximum and minimum energy that must remain in the ESD. Equation (21), together with (17) and (18), limits the number of changes in the operation state of the ESDs during the period of analysis, where for d = 0, esd i,0 is the initial value of esd i,d. Equation (22) imposes that esd i,d must be binary. 6) Additional Operational Constraints: Const- raints (23)–(25) are also part of the model ( PS i,d )2 + ( QS i,d )2 ≤ ( S̄S i )2 ∀i ∈ �S,∀d ∈ �d (23) V ≤ Vi,d ≤ V̄ ∀i ∈ �n,∀d ∈ �d (24) 0 ≤ Iij,d ≤ Īij ∀ij ∈ �b,∀d ∈ �d. (25) Constraint (23) models the generation capacity of the substation, (24) defines the voltage limits of the nodes, and (25) confines the current in the conductors to its maximum supported current. C. Mixed-Integer Nonlinear Programing Model The objective in this paper is to minimize the energy pur- chase cost of the DSS and dispatchable DGs, as shown in (26). Equations (27) and (28) are obtained by adding the contribu- tions of all of the devices cited previously to the active and reactive power balance equations (1) and (2). The model for the OODN problem is described min ∑ d∈�d ∑ i∈�S cS d�tPS i,d + ∑ d∈�d ∑ i∈�dg cdg d �tPdg i,d (26) subject to constraints (3)–(25) ∑ ki∈�b Pki,d − ∑ ij∈�b ( Pij,d + RijI 2 ij,d ) + ∑ ki∈�vr Pvr ki,d − ∑ ij∈�vr Pvr ij,d + PS i,d + Pdg i,d + Prs i,d + Psd+ i,d − Psd− i,d = PD i,d ∀i ∈ �n,∀d ∈ �d (27)∑ ki∈�b Qki,d − ∑ ij∈�b ( Qij,d + XijI 2 ij,d ) + ∑ ki∈�vr Qvr ki,d − ∑ ij∈�vr Qvr ij,d + QS i,d + Qdg i,d + Qcb i + Qscb i,d = QD i,d ∀i ∈ �n,∀d ∈ �d. (28) The proposed model (26)–(28) is an MINLP problem due to the following: (3)–(5), (23), (27), and (28) have squared variables and the product of squared variables; (12) presents the product of two variables; moreover, (9), (14), and (21) are nonlinear because they present the sum of the absolute value of a difference in variables. The presence of integer and binary variables makes this problem even more difficult to solve. Existing methodologies are able to find only low-quality solutions, and convergence, even for a feasible solution, is not guaranteed. D. Mixed-Integer Second-Order Cone Programing Model This section presents an MISOCP model for the OODN problem considering ESDs. The process to obtain this model from (26) to (28) involves variable substitutions, convexification of constraints [23], and an equivalent dis- junctive formulation for the model of VRs and OLTCs. In (3), (4), (27), and (28), the voltage and current magnitudes 2830 IEEE TRANSACTIONS ON SMART GRID, VOL. 6, NO. 6, NOVEMBER 2015 are squared, so that the change in variables can be accomplished Vsqr i,d = V2 i,d (29) Isqr ij,d = I2 ij,d. (30) After this, (3), (27), and (28) become linear ∑ ki∈�b Pki,d − ∑ ij∈�b ( Pij,d + RijI sqr ij,d ) + ∑ ki∈�vr Pvr ki,d − ∑ ij∈�vr Pvr ij,d + PS i,d + Pdg i,d + Prs i,d + Psd+ i,d − Psd− i,d = PD i,d ∀i ∈ �n,∀d ∈ �d (31)∑ ki∈�b Qki,d − ∑ ij∈�b ( Qij,d + XijI sqr ij,d ) + ∑ ki∈�vr Qvr ki,d − ∑ ij∈�vr Qvr ij,d + QS i,d + Qdg i,d + Qcb i + Qscb i,d = QD i,d ∀i ∈ �n,∀d ∈ �d (32) Vsqr i,d − Vsqr j,d = 2 ( RijPij,d + XijQij,d )+ Z2 ijI sqr ij,d ∀ij ∈ �b,∀d ∈ �d. (33) Additionally, (24) and (25) are expressed in a linear form V2 ≤ Vsqr i,d ≤ V̄2 ∀i ∈ �n,∀d ∈ �d (34) 0 ≤ Isqr ij,d ≤ Ī2 ij ∀ij ∈ �b,∀d ∈ �d. (35) Conversely, (4) continues to be nonlinear because on its left- hand side there is a product of variables and on its right-hand side there are squared variables Isqr ij,dVsqr j,d = P2 ij,d + Q2 ij,d ∀ij ∈ �b,∀d ∈ �d. (36) However, convexity and optimality can be achieved by relaxing the equalities in (36) into the second-order conic constraints Isqr ij,dVsqr j,d ≥ P2 ij,d + Q2 ij,d ∀ij ∈ �b,∀d ∈ �d. (37) If the dual variables related to constraints (37) in the OODN problem are greater than zero, then the conic con- straint (37) is active and, therefore, is equivalent to the nonlinear constraint (36) [24]. A set of sufficient conditions for constraint (37) to be active is provided in [25]. As a result of the variable changes (29) and (30), (12) is squared on both sides, and the squared value of the tap is needed in (13). Also, (12) presents a product of two variables: one continuous and one integer. Expressions (12) and (13) are (38) and (39) after altering the variable Vsqr j,d = tsqr ij,dVsqr i,d ∀ij ∈ �vr,∀d ∈ �d (38) tsqr ij,d = 1 + 2rvr ij lvr ij,d l̄vr ij + ( rvr ij )2 ( lvr ij,d l̄vr ij )2 ∀ij ∈ �vr,∀d ∈ �d. (39) The integer variable lvr ij,d can be represented by a set of binary variables lvr ij,d = 2l̄vr ij∑ k=0 [( k − l̄vr ij ) bvr ij,d,k ] ∀ij ∈ �vr,∀d ∈ �d (40) 2l̄vr ij∑ k=0 bvr ij,d,k = 1 ∀ij ∈ �vr,∀d ∈ �d (41) bvr ij,d,k ∈ {0, 1} ∀ij ∈ �vr,∀d ∈ �d,∀k = 0 . . . 2l̄vr ij . (42) Using the representation in (40)–(42), the maximum varia- tion of taps during the analyzed period (14) is modeled ∑ d∈�d ∣∣∣∣∣∣∣ 2l̄vr ij∑ k=0 [( k − l̄vr ij ) bvr ij,d,k ] − 2l̄vr ij∑ k=0 [( k − l̄vr ij ) bvr ij,d−1,k ] ∣∣∣∣∣∣∣ ≤ �̄vr ij ∀ij ∈ �vr. (43) The squared value of the tap is shown tsqr ij,d = 2l̄vr ij∑ k=0 ⎡ ⎢⎣ ⎛ ⎝1 + rvr ij ( k − l̄vr ij ) l̄vr ij ⎞ ⎠ 2 bvr ij,d,k ⎤ ⎥⎦ ∀ij ∈ �vr,∀d ∈ �d. (44) The squared value of the regulated voltage can be calculated Vsqr j,d = 2l̄vr ij∑ k=0 ⎡ ⎢⎣ ⎛ ⎝1 + rvr ij ( k − l̄vr ij ) l̄vr ij ⎞ ⎠ 2 Vsqr i,d bvr ij,d,k ⎤ ⎥⎦ ∀ij ∈ �vr,∀d ∈ �d. (45) The nonlinear product Vsqr i,d bvr ij,d,k is represented by using the variable Vc i,j,d,k Vsqr j,d = 2l̄vr ij∑ k=0 ⎡ ⎢⎣ ⎛ ⎝1 + rvr ij ( k − l̄vr ij ) l̄vr ij ⎞ ⎠ 2 Vc i,j,d,k ⎤ ⎥⎦ ∀ij ∈ �vr,∀d ∈ �d. (46) Vc i,j,d,k = Vsqr i,d bvr ij,d,k ∀ij ∈ �vr,∀d ∈ �d,∀k = 0 . . . 2l̄vr ij . (47) Equation (46) is linear, but (47) presents the product of a continuous variable and a binary variable that can be sub- stituted by two equivalent linear constraints, according to the disjunctive formulation shown V2bvr ij,d,k ≤ Vc i,j,d,k ≤ V̄2bvr ij,d,k ∀ij ∈ �vr,∀d ∈ �d,∀k = 0 . . . 2l̄vr ij (48) V2 ( 1 − bvr ij,d,k ) ≤ Vsqr i,d − Vc i,j,d,k ≤ V̄2 ( 1 − bvr ij,d,k ) ∀ij ∈ �vr,∀d ∈ �d,∀k = 0 . . . 2l̄vr ij (49) bvr ij,d,k ∈ {0, 1} ∀ij ∈ �vr,∀d ∈ �d,∀k = 0 . . . 2l̄vr ij . (50) After all of the aforementioned changes, the model still presents three constraints that are nonlinear, i.e. (9), (14) and (21). However, note that (14) was replaced by (43). These constraints limit the number of changes in MACEDO et al.: OODNs CONSIDERING ESDs 2831 operation of the SCBs, VRs, and ESDs, respectively, during the period of analysis. A linear equivalent to constraint (9) is shown ∑ d∈�d ( n+ i,d + n− i,d ) ≤ �̄scb i ∀i ∈ �scb (51) nscb i,d − nscb i,d−1 = n+ i,d − n− i,d ∀i ∈ �scb,∀d ∈ �d (52) n+ i,d ≥ 0 ∀i ∈ �scb,∀d ∈ �d (53) n− i,d ≥ 0 ∀i ∈ �scb,∀d ∈ �d. (54) A similar linear equivalent can be obtained in (43). Constraints (55)–(58) represent a linear equivalent to (43) ∑ d∈�d ( l+ij,d + l−ij,d ) ≤ �̄vr ij ∀ij ∈ �vr (55) 2l̄vr ij∑ k=0 [( k − l̄vr ij ) bvr ij,d,k ] − 2l̄vr ij∑ k=0 [( k − l̄vr ij ) bvr ij,d−1,k ] = l+ij,d − l−ij,d ∀ij ∈ �vr,∀d ∈ �d (56) l+ij,d ≥ 0 ∀ij ∈ �vr,∀d ∈ �d (57) l−ij,d ≥ 0 ∀ij ∈ �vr,∀d ∈ �d. (58) For constraint (21), the linear equivalent is ∑ d∈�d ( e+ i,d + e− i,d ) ≤ �̄sd i ∀i ∈ �sd (59) esd i,d − esd i,d−1 = e+ i,d − e− i,d ∀i ∈ �sd,∀d ∈ �d (60) 0 ≤ e+ i,d ≤ 1 ∀i ∈ �sd,∀d ∈ �d (61) 0 ≤ e− i,d ≤ 1 ∀i ∈ �sd,∀d ∈ �d. (62) Now, after these substitutions and linearizations, it is possi- ble to write an MISOCP model for solving the OODN problem considering ESDs min ∑ d∈�d ∑ i∈�S cS d�tPS i,d + ∑ d∈�d ∑ i∈�dg cdg d �tPdg i,d (63) subject to (5)–(7), (8), (10), (11), (17)–(20), (22), (23), (31)–(35), (37), (46), (48)–(50), and (51)–(62). (64) The model (63) and (64) has linear, quadratic, and second- order cone constraints, all of them convex. It also has continu- ous and integer variables. The commercial solver CPLEX [26] ensures global convergence for the proposed model. E. Mixed-Integer Linear Programing Model In order to compare the results of the MISOCP model for the OODN problem considering ESDs, an MILP model is pre- sented below. Recall that (5), (23), and (37) are the nonlinear constraints of the model and (37) was obtained from (36) using conic relaxation. Instead of the presented approach, an MILP model can be obtained by linearizing these constraints. First, assume for (36) that the voltage magnitude interval [V2, V̄2] is small, so that (65) is valid to approximate the left-hand side of (36) Isqr ij,dVsqr j,d ≈ Isqr ij,d ( Vnom)2 ∀ij ∈ �b,∀d ∈ �d. (65) Quadratic terms from (5), (23), and (36) can be approxi- mated using piecewise linear functions ( Pdg i,d )2 + ( Qdg i,d )2 ≈ f ( Pdg i,d, S̄dg i , � ) + f ( Qdg i,d, S̄dg i , � ) ∀i ∈ �dg,∀d ∈ �d (66) ( PS i,d )2 + ( QS i,d )2 ≈ f ( PS i,d, S̄S i , � ) + f ( QS i,d, S̄S i , � ) ∀i ∈ �S,∀d ∈ �d (67) P2 ij,d + Q2 ij,d ≈ f ( Pij,d, V̄ Īij, � )+ f ( Qij,d, V̄ Īij, � ) ∀ij ∈ �b,∀d ∈ �d. (68) The piecewise linear approximation function f (y, ȳ, �) used in (66)–(68) is defined f (y, ȳ, �) = �∑ γ=1 σy,γ δy,γ (69) y = y+ + y− (70) y+ + y− = �∑ γ=1 δy,γ (71) 0 ≤ δy,γ ≤ ȳ/� ∀γ = 1, . . . , � (72) σy,γ = (2γ − 1)ȳ/� ∀γ = 1, . . . , � (73) y+, y− ≥ 0. (74) Equation (69) approximates the square value of variable y by using variables δy,γ in order to discretize the absolute value of y. The parameter σy,γ represents the slope of the γ th line segment of the linearization method and is used to evaluate the contribution of δy,γ in each step of the discretization. The bound ȳ/� in (72) and (73) represents the length of the � discretized segments. δy,γ is a set of variables that denotes the length of the discretized segments, whose sum is equal to |y| according to (71). y+ and y− are nonnegative variables used to model |y| and the parameter ȳ is the maximum value of y. The nonlinear constraints (5), (23), and (36) can be replaced by the linear constraints f ( Pdg i,d, S̄dg i , � ) + f ( Qdg i,d, S̄dg i , � ) ≤ ( S̄dg i )2 ∀i ∈ �dg,∀d ∈ �d (75) f ( PS i,d, S̄S i , � ) + f ( QS i,d, S̄S i , � ) ≤ ( S̄S i )2 ∀i ∈ �S,∀d ∈ �d (76) Isqr ij,d ( Vnom)2 = f ( Pij,d, V̄ Īij, � )+ f ( Qij,d, V̄ Īij, � ) ∀ij ∈ �b,∀d ∈ �d. (77) Finally, the MILP model for the problem of OODN consid- ering ESDs is given min ∑ d∈�d ∑ i∈�S cS d�tPS i,d + ∑ d∈�d ∑ i∈�dg cdg d �tPdg i,d (78) subject to (6), (7), (8), (10), (11), (17)–(20), (22), (31)–(35), (46), (48)–(50), (51)–(62), and (75)–(77). (79) 2832 IEEE TRANSACTIONS ON SMART GRID, VOL. 6, NO. 6, NOVEMBER 2015 Fig. 6. Eleven-node test system of Levron et al. [14]. Notice that the linearization of the quadratic and conic constraints increases the number of continuous variables, but the number of integer variables remains unchanged. Also, greater values of � improve the fitting of the piecewise lin- ear approximation function, but also increases the size of the problem. III. TESTS AND RESULTS An 11-node test system and a real 42-node system were used to demonstrate the performance and robustness of the proposed methodologies. For both systems, �t = 0.5 h, V = 0.95 p.u., and V̄ = 1.05 p.u. The proposed MISOCP and MILP models, given by (63)–(64) and (78)–(79), respectively, were imple- mented in the mathematical programing language AMPL [27] and solved using the commercial solver for optimization prob- lems CPLEX (version 12.6.0.0, 64 bits using default settings) on a computer with a 3.40 GHz Intel core i7 processor and 16 GB of RAM. Complete results are shown for the MISOCP model, while only numerical results are shown for the MILP model for comparison. A. 11-Node Test System In order to validate the proposed model, it was first analyzed in a slightly modified version of the microgrid proposed by Levron et al. [14]. Fig. 6 shows the test system containing two photovoltaic RSs with power peaks of 1 and 0.5 MW at nodes 8 and 5, respectively; two ESDs at nodes 4 and 11 with 0.4 MWh capacities each and a self-discharge rate of βsd i = 0.021; a CB of 0.3 MVAr at node 2; and five loads. The objective was to minimize the cost of energy purchase from node 1. Impedances are in percent with a base (power and voltage magnitude) equal to the transformer’s nominal ratings. It was also assumed that both ESDs had the same value of efficiency ηsd+ i = ηsd− i = 0.95. The maximum daily change in operation of the ESDs was equal to 2. The cost of energy cS d, the active and reactive total load power, and the active power of photovoltaic RSs (shown in Fig. 7) were the same as those utilized in [14]. They were generated over a 72-h period. Seven cases were analyzed with different parameter values. For each case, the values of the modified parameters were specified, while all other values remained the same as previously cited. Complete results are shown for Cases A1 and A2, while for Cases A3–A7 the Fig. 7. Eleven-node test system of Levron et al. [14]. (a) Cost of energy purchase from node 1. (b) Total active load power. (c) Total reactive load power. (d) Active power of photovoltaic RSs. TABLE I RESULTS FOR CASES USING THE 11-NODE TEST SYSTEM results are reported in Table I. The cases are described as follows. Case A1: Without ESDs in the system and consid- ering a transformer with increased capacity: S̄S 1 = 10 MVA. Case A2: With both ESDs in the system and considering the original transformer capacity: S̄S 1 = 5 MVA. In order to analyze the influence of different characteris- tics of the ESDs in the system, the following cases consider a transformer with increased capacity: S̄S 1 = 10 MVA. Case A3: With the ESDs and considering a transformer with increased capacity (S̄S 1 = 10 MVA). Case A4: Given ideal ESDs, i.e., ηsd+ i = ηsd− i = 1 and βsd i = 0, in this case, the ESDs will not present losses. Case A5: ESDs at nodes 4 and 11 with 0.8 MWh each. This case considers ESDs with twice the capacity of the ESDs of Case A3. Case A6: Only two changes in operation for each ESD during the analyzed period. Case A7: Double the self-discharge for each ESD, i.e., βsd i = 0.042. The optimal solutions for Cases A1 and A2 are shown in Fig. 8. Fig. 8(a) shows the total energy stored by the ESDs in Case A2. As expected, the energy was stored when the MACEDO et al.: OODNs CONSIDERING ESDs 2833 Fig. 8. Optimal operation of the 11-node test system. (a) Total stored energy (Case A2). (b) Total active power extracted and injected by the ESDs (Case A2). (c) Apparent power generation at node 1, SS d , without and with the ESDs in the system (Cases A1 and A2). (d) Active power losses in the lines without and with the ESDs in the system (Cases A1 and A2). (e) Active power losses in the ESDs (Case A2). (f) Voltage profile on the load nodes with the highest and lowest voltages during the operation period, Cases A1 and A2 (same profiles). energy cost was low and injected into the DN when the energy cost was high. Moreover, it is easy to see the rela- tionship between Fig. 8(a) and (b): when power is extracted from the system [negative values in Fig. 8(b) between 5:30 and 6:30 A.M. on all three days of Case A2], the amount of stored energy increases, meeting the maximum value of 0.8 MWh; when power is injected into the system [positive val- ues in Fig. 8(b) between 8:30 and 9:30 A.M. on days 1 and 3, and between 7:30 and 10:00 A.M. on day 2], the stored energy decreases. Because the maximum number of daily changes in operation for the ESD is two, there is a total of six periods of operation for the ESDs—three of them extract- ing and three of them injecting power. Fig. 8(c) shows the apparent power generation of the DSS with (Case A2) and without (Case A1) the ESDs in the system. Without the ESDs, the system operation violates the capacity of the transformer considered in Case A2, and the problem is infeasible (for Case A2). Fig. 8(d) and (e) shows the power losses in the lines (for Cases A1 and A2) and in the ESDs (Case A2), respectively. Note that when the ESDs are extracting/injecting power in the system, the losses increase. When they are charged, due to self-discharge, there is a small amount of loss. Fig. 8(f) shows that all voltages in the system are within their limits (Case A2). Constraint (37) was also active for all branches in both cases. The results for Cases A1–A7 are shown in Table I, which contains information about total energy losses in lines and ESDs, and the cost of energy purchase Fig. 9. Forty-two-node real system. (value of the objective function). Comparing the results of Cases A1 and A2, as shown in Table I, an economy of energy purchase of 0.82% is obtained when the ESDs are present in the system. With a transformer with increased capacity (Case A3), the economy is 0.85%, when compared to Case A1. If the ESDs were ideal (Case A4), the economy would be 1.02%. With ESDs possessing twice the capacity (Case A5), the economy is 1.69%. Limiting the changes, the operation with two changes in the period (Case A6) reduces the econ- omy to 0.29%. Finally, doubling the self-discharge rate of both ESDs (Case A7) reduces the economy to 0.80%. Also, note that the results obtained with the MISOCP and MILP models are almost identical. The minimum and the maximum voltage magnitudes in the system (on load nodes) for all cases were 0.9621 and 0.9985 p.u., respectively. B. 42-Node Real System Fig. 9 shows the second system analyzed. It is a 42-node real DN [28] containing two OLTCs with five tap positions (−2; −1; 0; +1; +2) and a regulation of 2.5% per tap (±5% of maximum regulation); two VRs with 33 tap posi- tions (−16;−15;−14; . . . ; 0;+1;+2; . . . ;+16) and a regu- lation of 0.625% per tap (±10% of maximum regulation); five SCBs of 480 kvar with four units (120 kvar/unit) at nodes 38, 39, and 40, and of 1200 kvar (300 kvar/unit) at nodes 41 and 42; two biomass dispatchable DGs, both with a nominal capacity of 1 MVA and capacitive and inductive minimum power factors of 1 and 0.95, respectively; two wind turbines, three photovoltaic generators, and eight ESDs with a capacity of 0.8 MWh each. It was assumed that �̄scb i = 4 per day for all SCBs and �̄vr ij = 4 for all OLTCs, and �̄vr ij = 11 for all VRs. The objective is to minimize the cost of energy purchased from the DSS at node 1 and from dispatchable DGs 2 and 3. The voltage at the DSS is 13.8 kV, and its capacity is 50 MVA. Fig. 10(a) shows the cost of energy, cS d. Figs. 10(b)-(c) show the active and reactive total load power. Figs. 10(d)-(e) show the active power of RSs. These were generated over a 72-h period. We also assumed an efficiency of ηsd+ i = ηsd− i = 0.95, βsd i = 0.01, Psd+ i = Psd− i = Esd i = 0, and P̄sd+ i,d = P̄sd− i,d = 800 kW for all ESDs. Two cases were analyzed: without the presence of the ESDs in the system (Case B1) and considering the ESDs (Case B2). Results for both cases are presented. Case B1: The global optimal solution for the model (63) and (64) for the 42-node DN without the 2834 IEEE TRANSACTIONS ON SMART GRID, VOL. 6, NO. 6, NOVEMBER 2015 Fig. 10. Forty-two-node real system. (a) Cost of energy purchase from node 1 and dispatchable DGs. (b) Total active load power. (c) Total reactive load power. (d) Active power of photovoltaic RSs. (e) Active power of wind turbines RSs. ESDs is shown in Fig. 11. Constraint (37) was also active for all branches. Fig. 11 shows the sum of reactive power injected by all the SCBs, the tap movements of the VRs, the total active and reactive power injected by the dispatchable DGs, and the voltage profile of the system. Fig. 11(a) shows that more units of SCBs were connected when the load was heavy. The taps of both OLTCs remained at position +2 for all load levels. Note that the voltages were always within their limits. Case B2: Figs. 12 and 13 show the optimal solution for the 42-node DN considering the ESDs. Constraint (37) was active for all branches. Figs. 11 and 12 contain the same information, from which one can draw the same conclusions. Fig. 13 shows the total energy stored in the ESDs, total injected and extracted power by the ESDs, the apparent power injected by the DSS, the active power losses in the lines, the active power losses in the ESDs, and the voltage profile of the system. Again, as in the previous test system, energy was stored when the energy price was low and injected into the system when energy was expensive. Due to the self-discharge rate, there was a small decrease in stored energy over the three days, as can be seen in Fig. 13(a). Likewise, all voltages in the system were within their limits, as can be seen in Fig. 13(f). In Fig. 13(a), the numbered arrows indicate the phases of operation of the ESDs. During inter- vals 1 and 2, the ESDs were discharged. From 2, they started charging, reaching full charge at 3. In intervals 3 and 4, the ESDs neither extracted nor injected power in the system; so, due to self-discharge effect, they lost energy. At 4, they began extracting energy from the system, acquiring again full charge. During intervals 5 and 6, power was injected into the system, so the ESDs discharged. Afterwards, they remained uncharged Fig. 11. Optimal operation of the 42-node real system without ESDs (Case B1). (a) Total reactive power injected by the SCBs. (b) Tap move- ments of the VRs. (c) Total active power injected by the dispatchable DGs. (d) Total reactive power injected by the dispatchable DGs. (e) Voltage profile of the load nodes with the highest and lowest voltage during the operation period. Fig. 12. Optimal operation of the 42-node real system with ESDs (Case B2). (a) Total reactive power injected by the SCBs. (b) Tap movements of the VRs. (c) Total active power injected by the dispatchable DGs. (d) Total reactive power injected by the dispatchable DGs. during the periods 6 and 7; subsequently, they repeated a similar cycle. From Fig. 13(a), it can be seen that this cycle was similar for all three days. MACEDO et al.: OODNs CONSIDERING ESDs 2835 Fig. 13. Optimal operation of the 42-node real system with ESDs. (a) Total stored energy (Case B2). (b) Total active power extracted and injected by the ESDs (Case B2). (c) Apparent power generated at node 1, SS d , without and with the ESDs in the system (Cases B1 and B2). (d) Active power losses in the lines without and with the ESDs in the system (Cases B1 and B2). (e) Active power losses in the ESDs (Case B2). (f) Voltage profile of the load nodes with the highest and lowest voltage during the operation period (Case B2). TABLE II RESULTS FOR CASES USING THE 42-NODE REAL SYSTEM It must be clarified that the ESDs do not necessarily inject or extract power at the same time, i.e., at any given time, some ESDs may inject power, while others extract power. The mathematical model defines the optimal operation for each device. Further results for Cases B1 and B2 are shown in Table II. The values of total energy losses in lines and ESDs, and the cost of energy purchase are presented. The results in Table II show that an economy in energy pur- chase cost of 1.85% is obtained with the ESDs in the system (using the MISOCP model). The operation cost for the results of the MISOCP model and the MILP model for this case are also close. Finally, it is emphasized that the model is flexi- ble since the maximum number of changes in operation for TABLE III COMPUTATIONAL TIME NEEDED TO SOLVE THE MODELS FOR CASES USING THE 11-NODE TEST SYSTEM TABLE IV COMPUTATIONAL TIME NEEDED TO SOLVE THE MODELS FOR CASES USING THE 42-NODE REAL SYSTEM every piece of equipment is adjustable. Thus, if a solution that requires all controls to be adjusted only once over a period is desired, the proposed model can accommodate it. C. Computational Performance Below, further information is presented regarding the time needed to solve the problems presented in this paper, con- sidering both MISOCP and MILP models. The computational time needed to solve Cases A2–A7 using the 11-node system is shown in Table III. Since Case A1 was a simple power flow solution for each load level, the time to solve it has been omitted. Table IV shows the time needed to solve Cases B1 and B2 for the 42-node system. The results shown in Tables I–IV indicate that the MISOCP makes it possible to find the optimal solution of the problem in a reasonable amount of time. However, an MILP is an interest- ing alternative if a good quality solution is needed in a shorter amount of time. For the tests carried out, the time expended by the MILP was about half that required by the MISOCP; in real applications, this difference in the processing time can be taken into account. IV. CONCLUSION This paper presented MISOCP and MILP models for solv- ing the problem of OODNs considering ESDs. The proposed models considered as control variables the generation of active and reactive power of dispatchable DGs, the number of units in operation in each SCB, the tap position of the VRs and OLTCs, and the operation state of the ESDs. The objective was to minimize the total cost of energy purchase from the DSS and the dispatchable DGs. The steady- state operation of the radial DN was modeled through the use of linear and second-order conic expressions. The use of an MISOCP model guarantees convergence to optimality using classical optimization techniques. An 11-node test system and a 42-node real distribution sys- tem were used to demonstrate the accuracy of the mathemati- cal models, as well as the efficiency of the proposed solution 2836 IEEE TRANSACTIONS ON SMART GRID, VOL. 6, NO. 6, NOVEMBER 2015 techniques. The results showed that it is possible to determine the OODNs with ESDs using a flexible, realistic, and pre- cise MISOCP model. The approximated MILP model provides good-quality solutions with shorter computational time. ACKNOWLEDGMENT The authors would like to thank Y. Levron for his help in providing the test system data and results. REFERENCES [1] A. B. M. Shawkat Ali, Ed., Smart Grids: Opportunities, Developments, and Trends, Green Energy and Technology. New York, NY, USA: Springer, 2013. [2] T. Senjyu, Y. Miyazato, A. Yona, N. Urasaki, and T. Funabashi, “Optimal distribution voltage control and coordination with distributed generation,” IEEE Trans. Power Del., vol. 23, no. 2, pp. 1236–1242, Apr. 2008. [3] A. Augugliaro, L. Dusonchet, S. Favuzza, and E. R. Sanseverino, “Voltage regulation and power losses minimization in automated dis- tribution networks by an evolutionary multiobjective approach,” IEEE Trans. Power Syst., vol. 19, no. 3, pp. 1516–1527, Aug. 2004. [4] Y. P. Agalgaonkar, B. C. Pal, and R. A. Jabr, “Distribution volt- age control considering the impact of PV generation on tap changers and autonomous regulators,” IEEE Trans. Power Syst., vol. 29, no. 1, pp. 182–192, Jan. 2014. [5] R. A. Jabr, “Minimum loss operation of distribution networks with pho- tovoltaic generation,” IET Renew. Power Gen., vol. 8, no. 1, pp. 33–44, Jan. 2014. [6] R. R. Gonçalves, R. P. Alves, J. F. Franco, and M. J. Rider, “Operation planning of electrical distribution systems using a mixed integer linear model,” J. Control Autom. Elect. Syst., vol. 24, no. 5, pp. 668–679, Jun. 2013. [7] R. A. Araujo, P. C. M. Meira, and M. C. de Almeida, “Algorithms for operation planning of electric distribution networks,” J. Control Autom. Elect. Syst., vol. 24, no. 3, pp. 377–387, Jun. 2013. [8] F. A. Chacra, P. Bastard, G. Fleury, and R. Clavreul, “Impact of energy storage costs on economical performance in a distribution substation,” IEEE Trans. Power Syst., vol. 20, no. 2, pp. 684–691, May 2005. [9] Y. Levron and D. Shmilovitz, “Power systems’ optimal peak-shaving applying secondary storage,” Elect. Power Syst. Res., vol. 89, pp. 80–84, Aug. 2012. [10] Y. Levron and D. Shmilovitz, “Optimal power management in fueled systems with finite storage capacity,” IEEE Trans. Circuits Syst. I, Reg. Papers, vol. 57, no. 8, pp. 2221–2231, Aug. 2010. [11] J. P. Barton and D. G. Infield, “Energy storage and its use with inter- mittent renewable energy,” IEEE Trans. Energy Convers., vol. 19, no. 2, pp. 441–448, Jun. 2004. [12] T. K. A. Brekken et al., “Optimal energy storage sizing and control for wind power applications,” IEEE Trans. Sustain. Energy, vol. 2, no. 1, pp. 69–77, Jan. 2011. [13] P. Poonpun and W.T. Jewell, “Analysis of the cost per kilowatt hour to store electricity,” IEEE Trans. Energy Convers., vol. 23, no. 2, pp. 529–534, Jun. 2008. [14] Y. Levron, J. M. Guerrero, and Y. Beck, “Optimal power flow in micro- grids with energy storage,” IEEE Trans. Power Syst., vol. 28, no. 3, pp. 3226–3234, Aug. 2013. [15] A. Gabash and P. Li, “Active-reactive optimal power flow in distribution networks with embedded generation and battery storage,” IEEE Trans. Power Syst., vol. 27, no. 4, pp. 2026–2035, Nov. 2012. [16] A. Gabash and P. Li, “Flexible optimal operation of battery storage systems for energy supply networks,” IEEE Trans. Power Syst., vol. 28, no. 3, pp. 2788–2797, Aug. 2013. [17] S. Gill, I. Kockar, and G. W. Ault, “Dynamic optimal power flow for active distribution networks,” IEEE Trans. Power Syst., vol. 29, no. 1, pp. 121–131, Jan. 2014. [18] D. Pozo, J. Contreras, and E. E. Sauma, “Unit commitment with ideal and generic energy storage units,” IEEE Trans. Power Syst., vol. 29, no. 6, pp. 2974–2984, Nov. 2014. [19] F. Glover and G. Kochenberger, Handbook of Metaheuristics. Norwell, MA, USA: Kluwer Academic, 2003. [20] J. F. Franco, M. J. Rider, M. Lavorato, and R. Romero, “Optimal con- ductor size selection and reconductoring in radial distribution systems using a mixed-integer LP approach,” IEEE Trans. Power Syst., vol. 28, no. 1, pp. 10–20, Feb. 2013. [21] D. Shirmohammadi, H. W. Hong, A. Semlyen, and G. X. Luo, “A compensation-based power flow method for weakly meshed distribu- tion and transmission networks,” IEEE Trans. Power Syst., vol. 3, no. 2, pp. 753–762, May 1988. [22] R. G. Cespedes, “New method for the analysis of distribution networks,” IEEE Trans. Power Del., vol. 5, no. 1, pp. 391–396, Jan. 1990. [23] R. A. Jabr, “Radial distribution load flow using conic programming,” IEEE Trans. Power Syst., vol. 21, no. 3, pp. 1458–1459, Aug. 2006. [24] J. F. Franco, M. J. Rider, and R. Romero, “A mixed-integer quadratically- constrained programming model for the distribution system expansion planning,” Int. J. Elect. Power Energy Syst., vol. 62, pp. 265–272, Nov. 2014. [25] M. Farivar and S. H. Low, “Branch flow model: Relaxations and convexification—Part I,” IEEE Trans. Power Syst., vol. 28, no. 3, pp. 2554–2564, Aug. 2013. [26] (Apr. 10, 2014). IBM ILOG CPLEX. IBM. [Online]. Available: http://www-01.ibm.com/software/integration/optimization/ cplex-optimization-studio [27] R. Fourer, D. M. Gay, and B. W. Kernighan, AMPL: A Modeling Language for Mathematical Programming. 2nd ed. Pacific Grove, CA, USA: Brooks/Cole-Thomson Learning, 2003. [28] LaPSEE Power System Test Cases Repository. (2014, Oct. 27). [Online]. Available: http://www.feis.unesp.br/#!/lapsee Leonardo H. Macedo (S’14) received the B.Sc. and M.Sc. degrees from the Universidade Estadual Paulista, Ilha Solteira, São Paulo, Brazil, in 2012 and 2015, respectively, where he is currently pursuing the Ph.D. degree, all in electrical engineering. His current research interests include the development of methodologies for the optimization, planning, and control of electrical power systems. John F. Franco (S’11–M’13) received the B.Sc. and M.Sc. degrees from the Universidad Tecnológica de Pereira, Pereira, Colombia, in 2004 and 2006, respectively, and the Ph.D. degree from the Universidade Estadual Paulista, Ilha Solteira, São Paulo, Brazil, in 2012, all in electrical engineering. He is currently a Post-Doctoral Research Fellow with the Universidade Estadual Paulista. His current research interests include the development of methodologies for the optimization and planning of distribution systems. Marcos J. Rider (S’97–M’06) received the B.Sc. (Hons.) and P.E. degrees from the National University of Engineering, Lima, Perú, in 1999 and 2000, respectively; the M.Sc. degree from the Federal University of Maranhão, Maranhão, Brazil, in 2002; and the Ph.D. degree from the University of Campinas, Campinas, Brazil, in 2006, all in electrical engineering. He is currently a Professor with the Electrical Engineering Department, Universidade Estadual Paulista, Ilha Solteira, Brazil. His current research interests include the development of methodologies for the optimization, planning, and control of electrical power systems and applications of artificial intelligence in power systems. Rubén Romero (M’93–SM’08) received the B.Sc. and P.E. degrees from the National University of Engineering, Lima, Perú, in 1978 and 1984, respec- tively, and the M.Sc. and Ph.D. degrees from the University of Campinas, Campinas, Brazil, in 1990 and 1993, respectively, all in electrical engineering. He is currently a Professor with the Electrical Engineering Department, Universidade Estadual Paulista, Ilha Solteira, Brazil. His current research interests include the area of electrical power systems planning. << /ASCII85EncodePages false /AllowTransparency false /AutoPositionEPSFiles false /AutoRotatePages /None /Binding /Left /CalGrayProfile (Gray Gamma 2.2) /CalRGBProfile (sRGB IEC61966-2.1) /CalCMYKProfile (U.S. Web Coated \050SWOP\051 v2) /sRGBProfile (sRGB IEC61966-2.1) /CannotEmbedFontPolicy /Warning /CompatibilityLevel 1.4 /CompressObjects /Off /CompressPages true /ConvertImagesToIndexed true /PassThroughJPEGImages true /CreateJobTicket false /DefaultRenderingIntent /Default /DetectBlends true /DetectCurves 0.0000 /ColorConversionStrategy /LeaveColorUnchanged /DoThumbnails false /EmbedAllFonts true /EmbedOpenType false /ParseICCProfilesInComments true /EmbedJobOptions true /DSCReportingLevel 0 /EmitDSCWarnings false /EndPage -1 /ImageMemory 1048576 /LockDistillerParams true /MaxSubsetPct 100 /Optimize true /OPM 0 /ParseDSCComments false /ParseDSCCommentsForDocInfo false /PreserveCopyPage true /PreserveDICMYKValues true /PreserveEPSInfo false /PreserveFlatness true /PreserveHalftoneInfo true /PreserveOPIComments false /PreserveOverprintSettings true /StartPage 1 /SubsetFonts false /TransferFunctionInfo /Remove /UCRandBGInfo /Preserve /UsePrologue false /ColorSettingsFile () /AlwaysEmbed [ true /Arial-Black /Arial-BoldItalicMT /Arial-BoldMT /Arial-ItalicMT /ArialMT /ArialNarrow /ArialNarrow-Bold /ArialNarrow-BoldItalic /ArialNarrow-Italic /ArialUnicodeMS /BookAntiqua /BookAntiqua-Bold /BookAntiqua-BoldItalic /BookAntiqua-Italic /BookmanOldStyle /BookmanOldStyle-Bold /BookmanOldStyle-BoldItalic /BookmanOldStyle-Italic /BookshelfSymbolSeven /Century /CenturyGothic /CenturyGothic-Bold /CenturyGothic-BoldItalic /CenturyGothic-Italic /CenturySchoolbook /CenturySchoolbook-Bold /CenturySchoolbook-BoldItalic /CenturySchoolbook-Italic /ComicSansMS /ComicSansMS-Bold /CourierNewPS-BoldItalicMT /CourierNewPS-BoldMT /CourierNewPS-ItalicMT /CourierNewPSMT /EstrangeloEdessa /FranklinGothic-Medium /FranklinGothic-MediumItalic /Garamond /Garamond-Bold /Garamond-Italic /Gautami /Georgia /Georgia-Bold /Georgia-BoldItalic /Georgia-Italic /Haettenschweiler /Helvetica /Helvetica-Bold /HelveticaBolditalic-BoldOblique /Helvetica-BoldOblique /Impact /Kartika /Latha /LetterGothicMT /LetterGothicMT-Bold /LetterGothicMT-BoldOblique /LetterGothicMT-Oblique /LucidaConsole /LucidaSans /LucidaSans-Demi /LucidaSans-DemiItalic /LucidaSans-Italic /LucidaSansUnicode /Mangal-Regular /MicrosoftSansSerif /MonotypeCorsiva /MSReferenceSansSerif /MSReferenceSpecialty /MVBoli /PalatinoLinotype-Bold /PalatinoLinotype-BoldItalic /PalatinoLinotype-Italic /PalatinoLinotype-Roman /Raavi /Shruti /Sylfaen /SymbolMT /Tahoma /Tahoma-Bold /Times-Bold /Times-BoldItalic /Times-Italic /TimesNewRomanMT-ExtraBold /TimesNewRomanPS-BoldItalicMT /TimesNewRomanPS-BoldMT /TimesNewRomanPS-ItalicMT /TimesNewRomanPSMT /Times-Roman /Trebuchet-BoldItalic /TrebuchetMS /TrebuchetMS-Bold /TrebuchetMS-Italic /Tunga-Regular /Verdana /Verdana-Bold /Verdana-BoldItalic /Verdana-Italic /Vrinda /Webdings /Wingdings2 /Wingdings3 /Wingdings-Regular /ZapfChanceryITCbyBT-MediumItal /ZWAdobeF ] /NeverEmbed [ true ] /AntiAliasColorImages false /CropColorImages true /ColorImageMinResolution 200 /ColorImageMinResolutionPolicy /OK /DownsampleColorImages true /ColorImageDownsampleType /Bicubic /ColorImageResolution 300 /ColorImageDepth -1 /ColorImageMinDownsampleDepth 1 /ColorImageDownsampleThreshold 1.50000 /EncodeColorImages true /ColorImageFilter /DCTEncode /AutoFilterColorImages false /ColorImageAutoFilterStrategy /JPEG /ColorACSImageDict << /QFactor 0.76 /HSamples [2 1 1 2] /VSamples [2 1 1 2] >> /ColorImageDict << /QFactor 0.76 /HSamples [2 1 1 2] /VSamples [2 1 1 2] >> /JPEG2000ColorACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 15 >> /JPEG2000ColorImageDict << /TileWidth 256 /TileHeight 256 /Quality 15 >> /AntiAliasGrayImages false /CropGrayImages true /GrayImageMinResolution 200 /GrayImageMinResolutionPolicy /OK /DownsampleGrayImages true /GrayImageDownsampleType /Bicubic /GrayImageResolution 300 /GrayImageDepth -1 /GrayImageMinDownsampleDepth 2 /GrayImageDownsampleThreshold 1.50000 /EncodeGrayImages true /GrayImageFilter /DCTEncode /AutoFilterGrayImages false /GrayImageAutoFilterStrategy /JPEG /GrayACSImageDict << /QFactor 0.76 /HSamples [2 1 1 2] /VSamples [2 1 1 2] >> /GrayImageDict << /QFactor 0.76 /HSamples [2 1 1 2] /VSamples [2 1 1 2] >> /JPEG2000GrayACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 15 >> /JPEG2000GrayImageDict << /TileWidth 256 /TileHeight 256 /Quality 15 >> /AntiAliasMonoImages false /CropMonoImages true /MonoImageMinResolution 400 /MonoImageMinResolutionPolicy /OK /DownsampleMonoImages true /MonoImageDownsampleType /Bicubic /MonoImageResolution 600 /MonoImageDepth -1 /MonoImageDownsampleThreshold 1.50000 /EncodeMonoImages true /MonoImageFilter /CCITTFaxEncode /MonoImageDict << /K -1 >> /AllowPSXObjects false /CheckCompliance [ /None ] /PDFX1aCheck false /PDFX3Check false /PDFXCompliantPDFOnly false /PDFXNoTrimBoxError true /PDFXTrimBoxToMediaBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXSetBleedBoxToMediaBox true /PDFXBleedBoxToTrimBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXOutputIntentProfile (None) /PDFXOutputConditionIdentifier () /PDFXOutputCondition () /PDFXRegistryName () /PDFXTrapped /False /CreateJDFFile false /Description << /CHS /CHT /DAN /DEU /ESP /FRA /ITA (Utilizzare queste impostazioni per creare documenti Adobe PDF adatti per visualizzare e stampare documenti aziendali in modo affidabile. I documenti PDF creati possono essere aperti con Acrobat e Adobe Reader 5.0 e versioni successive.) /JPN /KOR /NLD (Gebruik deze instellingen om Adobe PDF-documenten te maken waarmee zakelijke documenten betrouwbaar kunnen worden weergegeven en afgedrukt. De gemaakte PDF-documenten kunnen worden geopend met Acrobat en Adobe Reader 5.0 en hoger.) /NOR /PTB /SUO /SVE /ENU (Use these settings to create PDFs that match the "Recommended" settings for PDF Specification 4.01) >> >> setdistillerparams << /HWResolution [600 600] /PageSize [612.000 792.000] >> setpagedevice