Distribution Systems Analysis andOptimizationbyHamed AhmadiM.A.Sc., The University of Tehran, 2011a thesis submitted in partial fulfillmentof the requirements for the degree ofDoctor of Philosophyinthe faculty of graduate and postdoctoral studies(Electrical and Computer Engineering )The University of British Columbia(Vancouver)April 2015c© Hamed Ahmadi, 2015AbstractDistribution systems (DS) are the last stage of any large power system, delivering elec-tricity to the end-users. Conventionally, simplicity of DS operation has been a priorityover its optimality. However, with the recent advancements in the automation and mea-surement infrastructures, it is now possible to improve the efficiency of DS operation. Inthis dissertation, a load modeling procedure is proposed which takes advantage of thedata available at the smart meters. An algorithm is proposed to decompose the loadat each customer level using the smart meter measurements. The proposed load modelrepresents the voltage dependence of loads according to the load composition. Based onthe voltage-dependent load model, a linear power flow formulation is developed for DSanalysis. The linear current flow equations are then proposed which calculate the branchflows directly without requiring the nodal voltages. Sensitivity factors in terms of currenttransfer and branch outage distribution factors are also derived using the linear powerflow concept. The advantages of having a set of linear equations describing the systemstatics are demonstrated in a variety of DS optimization problems, such as topological re-configuration, capacitor placement, and volt-VAR optimization. Using the linear currentflow equations, the mixed-integer nonlinear programming problem of DS reconfiguration isreformulated into a mixed-integer quadratic/linear programming problem, which substan-tially reduces the computational burden of the nonlinear combinatorial problem. Besidesdeveloping a direct mathematical optimization approach, a fast heuristic method is alsodeveloped here for the minimum-loss network reconfiguration based on the minimum span-ning tree problem. This heuristic method provides a good suboptimal solution to initializethe direct mathematical optimization approaches such as branch-and-cut algorithm usedfor solving combinatorial problems. Based on planar graph theory, an efficient mathe-matical formulation for the representation of the radiality constraint in reconfigurationproblems is introduced. It is shown that this formulation is advantageous over the avail-able methods in terms of accuracy and computational efficiency. The proposed algorithmsare tested using a variety of DS benchmarks and promising results are achieved.iiPrefaceThe contributions pointed in this dissertation have led to a number of already published,accepted for publication, or currently under review publications in well-known journalsand conferences. My research work and all my publications have been done by me underthe supervision of Prof. Jose´ R. Mart´ı. The co-authors of publications have provided uswith constructive feedback. The outcomes of each chapter in terms of publications are asfollows.Chapter 2• H. Ahmadi, J. R. Mart´ı, “Load decomposition at smart meters level using Eigenloadsapproach,” IEEE Transactions on Power Systems, Accepted on Nov. 2014.Chapter 3• J. R. Mart´ı, H. Ahmadi, L. Bashualdo “Linear power flow formulation based on avoltage-dependent load model,” IEEE Transactions on Power Delivery, vol. 28, no.3, pp. 1682–1690, July 2013.• H. Ahmadi, J. R. Mart´ı, “Linear current flow equations with application to distri-bution systems reconfiguration,” IEEE Transactions on Power Systems, Acceptedon Sep. 2014.• H. Ahmadi, J. R. Mart´ı, “Sensitivity factors for distribution systems,” IEEE Power& Energy General Meeting, Denver, CO, Jul. 2015.• H. Ahmadi, J. R. Mart´ı, “Power flow formulation based on a mixed-linear and nonlin-ear system of equations,” in Proc. 13th International Conference on Environmentaland Electrical Engineering (EEEIC), Wroclaw, Poland, Nov. 2013.iiiChapter 4• H. Ahmadi, J. R. Mart´ı, “Distribution system optimization based on a linear powerflow formulation,” IEEE Transactions on Power Delivery, vol. 30, no. 1, pp. 25–33,Feb. 2015.• H. Ahmadi, J. R. Mart´ı, “Mathematical representation of radiality constraint in dis-tribution system reconfiguration problem,” International Journal of Electrical Power& Energy Systems, vol. 16, pp. 293–299, Jan. 2015.• H. Ahmadi, J. R. Mart´ı, “Minimum-loss network reconfiguration: a minimum span-ning tree problem,” Sustainable Energy, Grids & Networks, vol. 1, pp. 1–9, Mar.2015.• H. Ahmadi, A. Alsubaie, J. R. Mart´ı, “Distribution system restoration consideringcritical infrastructures interdependencies,” in Proc. IEEE Power & Energy GeneralMeeting, Washington, DC, Jul. 2014.• H. Ahmadi, J. R. Mart´ı, A. Moshref, “Piecewise linear approximation of generatorscost functions using max-affine functions,” in Proc. IEEE Power & Energy GeneralMeeting, Vancouver, Canada, Jul. 2013.Chapter 5• H. Ahmadi, J. R. Mart´ı, H. W. Dommel, “A framework for volt-VAR optimization indistribution systems,” IEEE Transactions on Smart Grid, Accepted on Nov. 2014.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixGlossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xivDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2.1 Load Disaggregation . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2.2 Power Flow Formulation for Distribution Systems . . . . . . . . . . 41.2.3 Distribution System Optimization . . . . . . . . . . . . . . . . . . . 61.3 Research Objectives and Anticipated Impacts . . . . . . . . . . . . . . . . . 102 Load Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.1 Load Decomposition at Smart Meters . . . . . . . . . . . . . . . . . . . . . 132.1.1 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.1.2 Load Classification . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.1.3 Laboratory Measurements . . . . . . . . . . . . . . . . . . . . . . . . 212.1.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 292.2 Voltage-Dependent Load Modeling . . . . . . . . . . . . . . . . . . . . . . . 34v2.2.1 Conventional Load Models . . . . . . . . . . . . . . . . . . . . . . . 342.2.2 Proposed Load Model . . . . . . . . . . . . . . . . . . . . . . . . . . 343 Linear Power Flow Analysis . . . . . . . . . . . . . . . . . . . . . . . . . 383.1 The Linear Power Flow Equations . . . . . . . . . . . . . . . . . . . . . . . 383.1.1 Voltage Variation in a Distribution Feeder . . . . . . . . . . . . . . . 403.1.2 A Norton Equivalent of Load Synthesis . . . . . . . . . . . . . . . . 413.1.3 The Linear Power Flow Formulation . . . . . . . . . . . . . . . . . . 413.1.4 Distributed Generation Modeling within the Linear Power Flow (LPF) 423.1.5 Multi-Region Linear Power Flow . . . . . . . . . . . . . . . . . . . . 433.2 Simulation Results-LPF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 443.3 The Linear Current Flow Equations . . . . . . . . . . . . . . . . . . . . . . 483.4 Sensitivity Factors for Distribution Systems . . . . . . . . . . . . . . . . . . 513.4.1 Current Transfer Distribution Factors . . . . . . . . . . . . . . . . . 523.4.2 Branch Outage Distribution Factors . . . . . . . . . . . . . . . . . . 533.4.3 Active/Reactive Power Injection . . . . . . . . . . . . . . . . . . . . 543.4.4 Performance Indices . . . . . . . . . . . . . . . . . . . . . . . . . . . 553.4.5 Simulation Results-Sensitivity Factors . . . . . . . . . . . . . . . . . 574 Optimal Network Reconfiguration . . . . . . . . . . . . . . . . . . . . . . 634.1 Network Reconfiguration Using Mixed-Integer Programming . . . . . . . . . 634.1.1 LPF-based Mixed-Integer Programming . . . . . . . . . . . . . . . . 634.1.2 LCF-based Mixed-Integer Programming . . . . . . . . . . . . . . . . 724.1.3 A MILP Formulation of the Network Reconfiguration Problem . . . 754.2 Mathematical Representation of Radiality Constraint . . . . . . . . . . . . . 774.2.1 Planar Graph . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 774.2.2 Dual Graph . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 784.2.3 The Spanning Tree Constraint . . . . . . . . . . . . . . . . . . . . . 794.2.4 Formulation of Radiality Constraints . . . . . . . . . . . . . . . . . . 804.2.5 Inadequacy of Existing Methods in Representing the Radiality Con-straint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 824.2.6 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 834.3 Network Reconfiguration Using Minimum Spanning Tree . . . . . . . . . . . 854.3.1 Background: Minimum Spanning Tree . . . . . . . . . . . . . . . . . 864.3.2 Reconfiguration as A Minimum Spanning Tree Problem . . . . . . . 874.3.3 Series Neighborhood Search . . . . . . . . . . . . . . . . . . . . . . . 894.3.4 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 914.3.5 Feeders Loading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94vi4.3.6 Computational Complexity . . . . . . . . . . . . . . . . . . . . . . . 955 Volt-VAR Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 985.1 Load Characterization For Volt-VAR Optimization . . . . . . . . . . . . . . 995.2 Volt-VAR Optimization Algorithm . . . . . . . . . . . . . . . . . . . . . . . 1005.2.1 System Modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1035.3 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1065.3.1 Test System Description . . . . . . . . . . . . . . . . . . . . . . . . . 1065.3.2 Loss Reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1075.3.3 Conservation Voltage Reduction . . . . . . . . . . . . . . . . . . . . 1096 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1156.1 Accomplished Research Objectives . . . . . . . . . . . . . . . . . . . . . . . 1156.1.1 Objectives 1 & 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1156.1.2 Objectives 3–5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1166.1.3 Objective 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1166.2 Potential Impacts of the Presented Contributions . . . . . . . . . . . . . . . 1176.3 Future Developments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1186.3.1 Load Composition Forecasting . . . . . . . . . . . . . . . . . . . . . 1186.3.2 Three-Phase Unbalanced Robust Power Flow Solution . . . . . . . . 1196.3.3 Distribution System Reconfiguration for Phase Balancing . . . . . . 1196.3.4 Considering Reconfiguration within the VVO Framework . . . . . . 119Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121A Piecewise Linearization of a Quadratic Function . . . . . . . . . . . . . 131B Linearization of a Convex Quadratic Constraint . . . . . . . . . . . . . 133C The 69-Node System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136viiList of TablesTable 2.1 Load Classification for a Typical Residential Customer . . . . . . . . . . 21Table 2.2 Load Parameters in (2.16) and (2.17) Obtained by Measurements . . . . 36Table 3.1 Error of Approximating a Constant-Power Load Model by (2.16) . . . . 43Table 3.2 Standard Voltage Ranges for Distribution Systems [90] . . . . . . . . . . 44Table 3.3 Dimensions of the Employed Test Systems . . . . . . . . . . . . . . . . . 45Table 3.4 CPU Times for 20 Runs of the LPF and the Implicit Z-Bus Method . . . 47Table 3.5 Changes in the Total Power Losses (kW) Due to a Branch Outage . . . . 62Table 4.1 Dimensions of the Test Systems . . . . . . . . . . . . . . . . . . . . . . . 68Table 4.2 Initial Configurations and Losses for Test Systems . . . . . . . . . . . . . 69Table 4.3 Optimal Solutions For System Reconfiguration . . . . . . . . . . . . . . . 71Table 4.4 Optimal Configurations for Test Systems . . . . . . . . . . . . . . . . . . 72Table 4.5 Results Obtained By Assuming A Limit On the Number of SwitchingActions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73Table 4.6 Computation Times (s) for the LCF-based Formulations . . . . . . . . . 77Table 4.7 Values for βij According to (4.15) for the Network in Fig. 4.6 . . . . . . 82Table 4.8 CPU Times When Different Methods are Used To Represent RadialityConstraint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84Table 4.9 MST Search For The 14-Node Graph Using Prim’s Algorithm . . . . . . 87Table 4.10 Comparison of The Best Possible Configuration and Meshed Network . . 93Table 4.11 Comparison of Network Losses Obtained by The Proposed Algorithmsand Other References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93Table 4.12 Radial Configurations Obtained by The LS Algorithm . . . . . . . . . . 94Table 4.13 Comparison of CPU Times for Different Methods . . . . . . . . . . . . . 97Table 5.1 List of Typical Appliances for a Residential Area . . . . . . . . . . . . . 100Table C.1 Power Flow Data for the 69-Node System. P in kW and Q in kVAR. . . 137viiiList of FiguresFigure 1.1 Overall view of the proposed research plan. . . . . . . . . . . . . . . . . 11Figure 2.1 Illustration of Principal Component Analysis (PCA). V1 and V2 arescaled for proper representation. . . . . . . . . . . . . . . . . . . . . . . 17Figure 2.2 Flowchart of near-real-time load decomposition algorithm. . . . . . . . . 19Figure 2.3 Current waveform of a vacuum cleaner and the ST of its instantaneouspower. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23Figure 2.4 Current waveform of a MAC laptop charger and the ST of its instanta-neous power. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24Figure 2.5 Current waveform of a T5 fluorescent lamp with electronic ballast andthe ST of its instantaneous power. . . . . . . . . . . . . . . . . . . . . . 25Figure 2.6 Results of recognition tests on devices not present in the training set (θis defined by (2.8)). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26Figure 2.7 Results of recognition tests on devices not present in the training set (θis defined by (2.8)). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27Figure 2.8 The S Transform of the vacuum cleaner reconstructed by a reducednumber of eigenloads (r is the number of used eigenloads). . . . . . . . 28Figure 2.9 Differences between the original and reconstructed images with reducednumber of eigenloads. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28Figure 2.10 Rate of success in the recognition process versus the number of discardedeigenloads. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29Figure 2.11 Current waveforms of a LCD; bottom: original; top: noisy by a randomGaussian noise with a zero-mean and 30% of RMS value of current asits variance. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30Figure 2.12 Effect of noise (as in Fig. 2.11) on recognition accuracy. . . . . . . . . . 30Figure 2.13 Comparison of load models for a 3-phase induction motor (460 V, 3-phase, 1.4 HP, 1725 RPM). . . . . . . . . . . . . . . . . . . . . . . . . . 37Figure 3.1 A simple 2-node distribution feeder. . . . . . . . . . . . . . . . . . . . . 40ixFigure 3.2 Phasor diagram for voltages in the 2-node distribution feeder. . . . . . . 40Figure 3.3 A generic part of distribution system: Norton equivalent. . . . . . . . . 41Figure 3.4 Voltage regions, divided by vertical dashed-lines, and the correspondingapproximate constant-power loads (parabolas). CZ and CI are shownin single-column matrices for each region. ξ+ and ξ− are the relativeapproximation errors (%) above and below 1 p.u., respectively. . . . . . 44Figure 3.5 Histograms of relative errors for the LPF solution for various CZ . . . . . 46Figure 3.6 Average relative errors for the LPF solutions as a function of CZ . . . . . 47Figure 3.7 Voltage magnitudes and angles for the 135-node test system obtainedby the LPF method and the implicit Z-bus method. . . . . . . . . . . . . 48Figure 3.8 Maximum relative error in nodal voltages (µ) obtained using the MR-LPF and LPF versus loading factor (λ) in the 873-node and 10476-nodesystems. The minimum nodal voltage (Vmin) is also shown on the rightaxis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49Figure 3.9 A generic part of distribution system: The´venin equivalent. . . . . . . . 50Figure 3.10 Modeling of a branch outage using fictitious nodal current injections. . . 54Figure 3.11 Initial configuration of the 33-node system [95]. The dotted-lines indi-cate open switches. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58Figure 3.12 Voltage profile indices obtained using (3.41) (denoted as LPF) and Newton-Raphson power flow solutions (denoted as NR-PF) when a 1 MW/M-VAR DG/capacitor is installed at each node. . . . . . . . . . . . . . . . 59Figure 3.13 Total power losses obtained using (3.43) (denoted as LPF) and Newton-Raphson power flow solutions (denoted as NR-PF) when a 1 MW/M-VAR DG/capacitor is installed at each node. . . . . . . . . . . . . . . . 60Figure 4.1 Active and reactive power losses versus load voltage dependence in the70-node test system (optimal configuration). . . . . . . . . . . . . . . . 70Figure 4.2 Total active and reactive power consumptions versus load voltage de-pendence in the 70-node test system (optimal configuration). . . . . . . 70Figure 4.3 The 9-node network. The letters show the faces of the graph. . . . . . . 78Figure 4.4 The Kuratowski’s two graphs. . . . . . . . . . . . . . . . . . . . . . . . . 79Figure 4.5 Dual graph (dotted lines) of the 9-node network. . . . . . . . . . . . . . 79Figure 4.6 The network obtained by applying the conventional radiality constraints. 83Figure 4.7 The 16-node graph and its MST. . . . . . . . . . . . . . . . . . . . . . . 88Figure 4.8 Flowchart of the MST+LS algorithm for network reconfiguration. . . . 90Figure 4.9 Visual examples of series and non-series neighbors. . . . . . . . . . . . . 90Figure 4.10 Voltage profiles of the 119-node and 873-node systems under initial andobtained configurations. . . . . . . . . . . . . . . . . . . . . . . . . . . . 95xFigure 4.11 Feeders currents for the MST+LS and the optimal solutions. . . . . . . 96Figure 5.1 Typical load profile for a residential area in a weekday. . . . . . . . . . . 101Figure 5.2 Voltage dependence of typical appliances [83]. . . . . . . . . . . . . . . . 102Figure 5.3 The pi-equivalent circuit of a ULTC-equipped transformer. . . . . . . . . 104Figure 5.4 Daily load curve for the 69-Node feeder. . . . . . . . . . . . . . . . . . . 107Figure 5.5 Total power losses at peak load. . . . . . . . . . . . . . . . . . . . . . . 108Figure 5.6 Tap ratio and capacitor banks statuses under different objectives forthe radial network. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110Figure 5.7 Total active power drawn from the substation. . . . . . . . . . . . . . . 111Figure 5.8 Total reactive power drawn from the substation. . . . . . . . . . . . . . 111Figure 5.9 Total active power losses. . . . . . . . . . . . . . . . . . . . . . . . . . . 111Figure 5.10 Minimum nodal voltage. . . . . . . . . . . . . . . . . . . . . . . . . . . . 112Figure 5.11 Reduction in total active power demand w.r.t. “No Control” case. . . . 113Figure 5.12 Reduction in total active power losses w.r.t. “No Control” case. . . . . . 113Figure 5.13 Tap ratio and capacitor banks statuses under different objectives forlooped network (Switch 50-59 is closed). . . . . . . . . . . . . . . . . . . 114Figure A.1 A standard 2-D quadratic function and its PWL approximation. . . . . 132Figure B.1 Hexagon approximations of a circle. . . . . . . . . . . . . . . . . . . . . 135Figure B.2 Error of polygon approximation of a circle. . . . . . . . . . . . . . . . . 135Figure C.1 The 69-node system [121]. Dotted line: normally-open switch. . . . . . . 138xiGlossaryAMI Advanced Metering Infrastructure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1ANN Artificial Neural Network . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3BODF Branch Outage Distribution Factor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6CPU Central Processing Unit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8CTDF Current Transfer Distribution Factor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51CVR Conservation Voltage Reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9CWT Continuous Wavelet Transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14DG Distributed Generation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5DMS Distribution Management System. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .8DS Distribution Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1EMS Energy Management System. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .8EPRI Electric Power Research Institution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98FPGA Field-Programmable Gate Array. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .33IM Induction Motor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .21KCL Kirchhoff’s Current Law. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .53LCF Linear Current Flow. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .10LODF Line Outage Distribution Factor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5LPF Linear Power Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10LS Local Search . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89MICP Mixed-Integer Conic Programming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7MILP Mixed-Integer Linear Programming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7MINLP Mixed-Integer Nonlinear Programming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7MIP Mixed-Integer Programming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71xiiMIQCP Mixed-Integer Quadratically-Constrained Programming . . . . . . . . . . . . . . . . . . . . . . . 8MR-LPF Multi-Region Linear Power Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43MST Minimum Spanning Tree . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85PCA Principal Component Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .15PF Power Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4PTDF Power Transfer Distribution Factor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5PWL Piecewise Linearization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75RAM Random Access Memory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33SOS Special-Ordered Set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66SS Sectionalizing Switch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7ST S Transform. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13STFT Short-Time Fourier Transform. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14TS Tie Switch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7ULTC Under-Load Tap-Changing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6VPI Voltage Profile Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55VVO Volt-VAR Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2xiiiAcknowledgmentsThere are countless people who helped me during the course of my Ph.D. studies and itwill not be feasible to mention all of them. I try my best to take advantage of this shortspace and thank as many as I can.I would like to thank my family who supported me by all means during the hard yearsof being away from them. My dear girlfriend, Jing, was always with me and cheered meup with her unconditional love. It was not possible to reach this point in my life withouttheir unhesitant support.The continuous support of Prof. Jose´ Mart´ı is priceless. As my supervisor, he neverhesitated to dedicate his time to hear me and guide me through the tough path of research.This work was supported, financially, by the Natural Sciences and Engineering ResearchCouncil of Canada (NSERC).I need to acknowledge the scientific supports and feedback provided by Dr. Ali Moshref(BBA) and Dr. Hermann W. Dommel (Professor Emeritus) in distribution systems mod-eling and analysis. Also, Dr. Mukesh Nagpal, Cheong K. Siew, and Lincol Bashualdofrom BC Hydro, who provided information about the present practices at the distributionsystem industry, deserve to be acknowledged.The staffs at the electrical and computer engineering department at UBC have beenhelpful in the course of my studies by providing guidelines and advice, when needed, toresolve my issues promptly. Also, I should thank the Faculty of Graduate and PostdoctoralStudies at UBC for offering me awards for academic excellent. The UBC Let’s Talk Scienceis a nonprofit organization for connecting university and schools. I should thank them,especially Sherie Duncan, for helping us reach out as volunteers.I would like to express my sincere gratitude toward my colleagues, Abdullah Alsubaieand Francis Therrien, who helped me look at subjects from different point of views. I alsoacknowledge the comments provided by the committee members, Dr. Srivastava and Dr.Vaahedi, and University Examiners, Dr. Sassani, Dr. Mirabbasi, and Dr. Romero-Ramos,who were instructive in making this dissertation a more complete and clear.xivDedicationThis dissertation is dedicated to my beloved family, my mother, Rababeh Ali-Ahmadi, myfather, Mahmoud Ahmadi, and my sisters, Fatemeh, Zahra, and Somayeh. None of thesecould be done without their continuous and unhesitant support.xvChapter 1Introduction1.1 MotivationDistribution Systems (DS) are the last part of the bulk power system, delivering electricityto the end-users. Conventionally, DS are built and operated with minimal automation andoptimization in order to keep the operation as simple as possible. With the recent progressin the area of Advanced Metering Infrastructure (AMI) and the smart grid initiatives, it istime to revisit the traditional operation strategies and look for opportunities to improvethe efficiency of DS.The smart grid movement has made numerous improvements in transmission systemsoperation and planning. The distribution system counterpart, however, has not receivedsufficient attentions in terms of developing near-real-time algorithms for improving itsreliability, efficiency, and resiliency. The aim of this dissertation is to introduce a newphilosophy in distribution system analysis and optimization. The available data fromsmart meters can be used to build an accurate model of DS and computationally efficientmethods can be developed to perform near-real-time optimization.A basic tool in most of DS analysis and optimization is the power flow solution. Theconventional methods, such as the Newton-Raphson method, suffer from the non-convexnonlinear structure of the formulation. In order to develop algorithms that are reliable,fast, and robust, a linear formulation is desirable. This is the motivation for the author todevelop a linear power flow model. Besides, the conventional power flow algorithms do nottake into account the voltage dependence of loads. This is particularly important for DSanalysis since the network is in the close vicinity of loads. To address this issue, the loadsare first identified at the meters level. The load composition for each node is then derivedat every time interval and its voltage dependence characteristic is determined based onthis composition. This idea motivated the author to create a load monitoring frameworkfor accurate load modeling.1The existing optimization techniques for DS need to be adopted to the emerging tech-nologies and measurement systems. The author was inspired by this fact to develop opti-mization algorithms that are fast and robust to meet the requirements of near-real-timeapplications. Within the course of developing optimization routines, many challenges havebeen faced in terms of mathematical formulations. For instance, having a set of currentflow solutions instead of voltages can significantly enhance the solution of the formulatedoptimization problem. As another example, imposing the radiality of the network topol-ogy in an optimization problem is not a trivial task and previous attempts are shown tofail to do so. These issues, which emerged during the process of research, motivated theauthor to seek further contributions.The emerging practices in the power industry, such as controlling the voltage andVAR resources, and the reportedly published promising outcomes have motivated theauthor to utilize the developed tools here in other applications such as Volt-VAR Opti-mization (VVO). The present techniques have some shortcomings, such as inaccurate loadmodels and sophisticated optimization algorithms, that need to be addressed. The authorwas motivated to find adequate solutions for the mentioned issues based on the developedtools within the present research.In the following, a literature review is provided on different topics covered in thisdissertation.1.2 Literature Review1.2.1 Load DisaggregationOne of the most important features of distribution systems is the load behavior, which isalso the main source of change in the state of the system during normal conditions. Athorough understanding of the load nature would make a great impact on the performanceof any control action in a system. In a particular case, in order to successfully managethe voltage and reactive power in the system, the load voltage dependencies are the mainplayers since the voltage level affects the active and reactive power consumptions. In or-der to reach a high level of accuracy in the load voltage dependence modeling, the loadcomposition at each node is needed. Since the measurements are usually provided in acomposite form, decomposition techniques are required to recognize the load components.The concept of load disaggregation was first used in the early 80’s for load forecastingpurposes [1]. Many studies have been conducted for load disaggregation using differenttechniques and/or different load features. To characterize the load signature, the plotsof active power P versus reactive power Q have been used in [2], which was commercial-ized under a nonintrusive appliance load monitoring system (NIALMS) [3]. An improved2version of this method was proposed in [4] taking into account the harmonic contents ofthe signals, which addresses one of the inadequacies of NIALMS. A voltage-current (V -I)trajectory was used as the load signature in [5] which has shown good performance indistinguishing between loads with similar patterns. A deeper vision of the V -I trajectoryfor load identification was introduced in [6], which improved the precision of predictionsby 10% as compared to [5]. The features used in [6] include P -Q, total harmonic dis-tortion, and V-I wave-shape. Four learning algorithms were also used in [6], includingArtificial Neural Network (ANN), ANN coupled with an evolutionary algorithm, SupportVector Machine, and Adaptive Boost. It was shown that the Adaptive Boost algorithmbased on V-I wave-shape produces the highest accuracy in the predictions. The Wavelettransform was utilized in [7] to decompose the measured signals into time and frequencydomains. Using a sampling frequency of 10 kHz, current waveforms were decomposed into6 wavelet components. This wavelet decomposition, similar to Fourier decomposition oftime series, can be used to classify loads.A thorough survey on the available methods for load disaggregation was done in [8].Based on this study, the data available from the smart meters have a high potential tobe useful for load disaggregation purposes. The utilization of smart meters data for loaddisaggregation has been recently reported, for example, in [9] and [10]. The smart meterdata was utilized in [9] and an Explicit-Duration Hidden-Markov Model with differentialobservations was applied to the data for detecting and estimating individual home appli-ance loads. An event window-based mechanism was proposed in [10] which uses the powerwaveforms and clusters the signatures based on some features in a selected time window.A real-time load monitoring approach was introduced in [11] which is based on ANN anduses P , Q, power factor, peak and RMS values of V and I.The Support Vector Machine algorithm as well as 5-Nearest Neighbors method wereemployed in [12] to classify loads using P , Q, and power factor data. An algorithm wasproposed in [13] to recognize major appliances in a household, namely water heater, base-board heater, washing machine, and refrigerator, which only uses P as the feature toclassify loads. When two loads have similar power consumption, the algorithm in [13] mayfail. A pattern recognition approach was introduced in [14] to decompose the energy useof major appliances in a household using P measurements. This method needs an on-sitetraining for the first time for about a week and it will work as long as no new applianceis added. A few possible load configurations can be found using the algorithm proposedin [15], which utilizes P and Q signals and applies a finite-state machine based on fuzzytransitions to classify loads. In an experiment, the method in [15] found two possibleconfigurations, one of which was the correct one.There are research groups worldwide concentrating on developing load disaggregationtechniques. The research group at Jin Wen University of Science and Technology in Taiwan3has recently published the results of their work based on multiple algorithms, such as ANN[16], Particle Swarm Optimization [17], and Wavelet Transform [18]. The results of theirstudy have shown noticeable recognition rates. The problem of load disaggregation waswell described in [19] and [20]. Different load features such as I, P , Q, harmonics, instan-taneous admittance waveform, instantaneous power waveform, eigenvalues, and switchingtransient waveform were used. Two disaggregation algorithms were also proposed basedon optimization or pattern recognition (based on ANN). A committee decision mechanismwas then adopted to render the best solution out of a candidate pool containing candidatesobtained by different features. This framework has shown good performance in terms ofaccuracy of the estimations.1.2.2 Power Flow Formulation for Distribution SystemsAn essential part of DS automation is to provide computationally efficient calculations foronline real-time automated actions. Power Flow (PF) solutions are a fundamental part ofthese calculations. There has been a large effort throughout the years towards enhancingpower flow calculations in terms of computational speed and convergence characteristics.At the transmission system level, fast and reliable methods such as Gauss-Seidel [21],Newton-Raphson [22], and fast decoupled PF [23] have been used for many years. Distri-bution systems, however, are different from transmission systems in a number of aspects,such as X/R ratio, the use of underground cables, radial structures, and the unbalanceof the phase currents. Due to these differences, PF algorithms developed for transmissionsystems often fail [24] or lose efficiency when applied to distribution systems. A numberof PF solution methods have been developed to account for the specific nature of DS, in-cluding the backward/forward sweep method [25].The backward/forward sweep method is an iterative algorithm which is relatively time-demanding, especially for real-sized DS. An improved version of this method was presentedin [26] for radial DS in which the linear proportional principle is adopted to find the ratio ofthe real and imaginary components of the specified voltages with respect to the calculatedvoltage at the substation bus during the forward sweep. Using this method in a typicaltest system, the execution time has been reduced by 35.7% compared to the conventionalbackward/forward sweep. The PF problem has also been formulated as biquadratic equa-tions in [27], which is still based on iterative computations of backward/forward sweeps.A direct method for PF solutions was proposed in [28] in which the loads are treated firstas current source injections, and simple matrix calculations and iterative computationsare used to find the nodal voltages. A recursive formulation was proposed in [29] whichincludes three nonlinear equations for each branch, called DistFlow equations, that can besolved using Newton’s method. By defining new variables in PF formulation, the problemwas converted into a convex optimization problem in [30], or specifically, a conic quadratic4problem; such problems can be solved by interior point methods. Current injection mod-els have also been proposed, for example, in [31]-[33]. The superior performance of thecurrent injection method against the backward/forward sweep method was shown in [34].In all the mentioned references on PF calculations, the loads are modeled as constantP-Q. However, as the system gets closer to the loads voltage levels, the voltage depen-dence of the actual loads becomes more important in the representation. For example, inthe framework of Voltage VAR Optimization (VVO), the load voltage dependence playsan inevitable role and the performance of the VVO algorithms is highly dependent on theaccuracy of the load models [35]. In addition, the introduction of Distributed Genera-tion (DG) in DS may change the unidirectional characteristics of the power flow, on whichsome of the previous methods have been built. Also, there might exist weakly-meshed DSin addition to the radial ones, which creates an additional limiting factor for some of theaforementioned approaches.Sensitivity AnalysisThe idea of sensitivity analysis in power systems has been widely used to avoid recalcu-lation of the power flow solution. In transmission systems, the parameters used in theseanalyses are the Power Transfer Distribution Factor (PTDF) and the Line Outage Dis-tribution Factor (LODF). PTDFs are defined as the changes in the line power flows dueto a change in power injection at a particular bus. LODFs are defined as the changes inthe line power flows due to the disconnection of a particular line [36]. The calculation ofthese sensitivity factors has gained more interest recently due to the need for fast on-linereadjustments in modern power systems.There are generally two approaches to calculate sensitivity factors for power systems.Considering a more realistic model, power flow equations form a nonlinear system of equa-tions. In order to find sensitivities, one need to find the Jacobian matrix at a particularsolution of the network, which yields sensitivity factors that are only valid for small varia-tions around the operating point [37], [38]. The second approach is to find an approximatelinear model that describes the system and find the sensitivity factors for the approximatelinear system. The decoupled power flow equations, for instance, were used in [39] tofind the sensitivity of reactive power flows to transmission line and power transformeroutages. A linear approximation of the power flow equations is the so-called DC powerflow, in which all voltage magnitudes are assumed to be one per-unit, line resistances areignored, and voltage angles are assumed to be small enough so that their sines are ap-proximately equal to the angles. In a DC power flow model, only active power flows canbe approximated and reactive power flows are not considered. Based on the linear DCpower flow equations, the PTDFs and LODFs can be derived for a transmission system [36].A generalization of LODF in the case of multiple simultaneous line outages in a network5was done in [40]. A direct way of calculating the sensitivity factors was discussed in [41].Besides analytical studies, there have been some efforts recently to derive the sensitivityfactors based on real-time data provided by phasor measurement units. For example, in[42], sensitivity factors were derived without the need for a power flow model.Despite broad discussions on derivations and applications of sensitivity factors at thetransmission level, there is rather limited work on the distribution system counterpart.Distribution systems have high R/X ratios, radial configurations, a mixture of cables andoverhead lines, and unbalanced loads. Due to these unique features, the power flow algo-rithms and sensitivity analysis derived for transmission systems are not always valid fordistribution systems. Using the adjoint network method, which is based on the applica-tion of Tellegen’s theorem to power systems, the authors of [43] derived the sensitivitiesof power losses and voltage magnitudes with respect to power injection at any node in thesystem. However, this method is only valid for radial distribution systems. Also, it doesnot consider the voltage dependence of the loads, which is an important consideration indistribution systems. The sensitivity factors derived in [43] cannot be used to calculatethe Branch Outage Distribution Factor (BODF).Similar to the case of transmission systems, there are many applications for sensitivityfactors in distribution systems. The problem of DG placement usually requires a knowl-edge of the impact of power injection at each node on certain quantities of the network,e.g., power losses or voltage profile. Capacitor placement, as an example of reactive powersources, also relies highly on the sensitivity of nodal voltages and system losses to the re-active power injection at each particular node. Network reconfiguration, for purposes suchas loss reduction and voltage profile improvement, has been a continuously-studied prob-lem in the area. The BODFs are substantially useful factors to improve the computationperformance of algorithms dealing with the reconfiguration problem. In many reconfig-uration methods, the algorithm starts with a meshed network, and then opens branchesone after another to reach a radial configuration, e.g., [44]. In some other approaches, thealgorithm starts with a radial topology and seeks a branch exchange that gives a betterobjective value, e.g., [29]. In such reconfiguration procedures, BODFs have great potentialapplications.1.2.3 Distribution System OptimizationThere are many control variables that affect the performance of a DS, such as adjustmentof Under-Load Tap-Changing (ULTC) transformers and switchable capacitor banks, DG,interruptible loads, and network configuration. Some widely used measures for evaluatingthe performance of DS are voltage profile, load balance among feeders and phases, networklosses, service restoration time, and security and reliability. Many studies have been con-ducted to find the optimum values for the control variables to achieve higher performance6in DS. A brief literature review is given in the following on each specific topic.ReconfigurationIt is more convenient to operate DS in radial configurations. Radiality, compared to otherpossibilities, leads to simpler control strategies, especially in protection. Nonetheless,there are normally-open switches between feeders in different locations that can createa loop. Similarly, there are normally-close switches that can disconnect one section of afeeder. The first group is called Tie Switch (TS) and the second group is called Section-alizing Switch (SS) [45]. By closing one TS and opening an appropriate SS, a new radialconfiguration can be achieved which may alter DS performance. This simple idea, namedbranch exchange algorithm, has previously been applied to find best configurations thatprovide minimum losses, e.g., [29], [46], and [47]. As a more detailed mathematical model,a binary variable can be designated to each branch admittance which will then appearin nonlinear AC power flow equations. This variable is zero if the branch is open andone otherwise. These binary variables along with the nonlinear power flow equations andobjective function form an optimization problem which belongs to the family of Mixed-Integer Nonlinear Programming (MINLP) problems. The combinatorial MINLP problemsare well-known to be computationally expensive to solve and, in large systems, almostintractable. Due to the lack of efficient mathematical approaches for solving this class ofoptimization problems within a reasonable time, heuristic algorithms have been widelyused in the literature to find fast, while suboptimal, solutions.There are many heuristics that have been applied to the problem of DS reconfigurationwith different characteristics. Genetic Algorithms have been utilized with improved char-acteristics in [48]-[50]. Other intelligent search algorithms have also been applied such asSimulated Annealing [51], Ant Colony [52], Harmony Search Algorithm [53], EvolutionaryAlgorithms [54] and [55], Artificial Neural Networks [56], and Particle Swarm Optimization[57]. Although some of these search algorithms have shown good performance, they mightnot be able to provide optimal solutions within the constraints of online applications. Inaddition, these methods may provide suboptimal solutions and their optimality is usuallynot guaranteed.Besides heuristic methods, there are also direct mathematical approaches available inthe literature dealing with the DS reconfiguration problem. The Benders decompositionhas been used in [58] to solve the MINLP problem by dividing it into two subproblems.Even though this approach has shown commendable performance, the master problem isformulated as an MINLP, for which the existing commercial solvers have no guarantee toprovide the optimum solution within a reasonable time.The authors of [59] have proposed two reformulations of the reconfiguration problem,a Mixed-Integer Conic Programming (MICP) version and a Mixed-Integer Linear Pro-7gramming (MILP) version. The computational efficiencies of both the MICP and the MILPformulations have been compared. It is shown that in order to reach the optimal or near-optimal solution, the MICP is more time-demanding than the MILP. However, in the caseof large networks, the Central Processing Unit (CPU) time is relatively long for both ofthese methods. In addition, there are extra variables and constraints in the formulationsthat make the size of the problem unnecessarily large. It has to be noted though that thisis an excellent study in the context of DS reconfiguration.The authors of [60] have proposed another reformulation of reconfiguration problemwhich utilizes the DistFlow method (proposed by Baran et al. in [29]) for power flow cal-culations. The problem is converted into three different types of well-known optimizationproblems, namely mixed-integer quadratic programming, Mixed-Integer Quadratically-Constrained Programming (MIQCP), and second-order cone programming. The computa-tion time for solving the first two problems has been shown to be relatively short; however,these formulations assume a fixed value for nodal voltages, which is not applicable whenthe system voltage profile has to be considered. The time required to solve the thirdformulation is drastically high and, therefore, not practical.The authors of [61] proposed a new reformulation of the network reconfiguration prob-lem using auxiliary variables and mixed-integer programming techniques. In their paper,the total power injected from the substation to the feeders is considered as the objectiveaiming at loss reduction. However, assuming the real-world loads, the actual demand ofvoltage-sensitive loads changes as one reconfigures the network (and implicitly alters thenodal voltages). Due to this fact, the considered objective function cannot be used forreal-world distribution systems for minimizing the power losses.Volt-VAR OptimizationIn the past couple of decades, a rapid movement towards DS automation has been startedby many power utilities worldwide to increase efficiency, reliability, resiliency, and powerquality. As the main part of DS automation, Distribution Management System (DMS) toolsare being introduced as counterparts to the Energy Management System (EMS) tools thatexist for the optimal operation of transmission systems. The DMS concept offers manyappealing functionalities, including state estimation, fault location, service restoration,and VVO. The VVO functionality is basically defined as the optimal control of the systemequipment (e.g., capacitor banks, voltage regulators, ULTCs, etc.) in order to satisfy ob-jectives such as minimizing the losses, improving voltage profile, etc.Many utilities have already deployed a number of VVO applications and promisingresults are being reported. Good examples are the Northern State Power Company [62],B.C. Hydro [63]-[64], Taiwan Power Company [65], and American Electric Power [66].Besides loss reduction and voltage profile improvement, VVO offers one more useful capa-8bility: the so-called Conservation Voltage Reduction (CVR). The CVR essentially refers toreducing the voltage magnitudes to the minimum allowable in order to reduce the demand,usually for peak shaving or energy saving purposes. It stems from the fact that loads arevoltage-dependent with a positive slope, i.e. reducing the voltage yields a reduction inboth active and reactive demands. An immediate reaction to this concept would be a con-cern towards the voltage quality at the customer terminals. The deployment of AMI hasprovided enormous visibility for almost all the nodes in the DS, facilitating the applicationof CVR with a high confidence in the quality of power delivery [64]. Practical studies haveshown that, upon proper implementation of CVR, a demand reduction of 1% to 6% in thetotal demand can be achieved [64].Several methodologies have been proposed for VVO in DS. An algorithm based on theoriented discrete coordinate descent method was proposed in [67], which gives priority tofinding a feasible but suboptimal solution in a short time rather than finding the globaloptimum, which may impose a high computational burden. A supervisory control schemewas developed in [68] that takes advantage of the measurements available at the substationto optimally adjust the voltage regulators and capacitor banks at the substation as well asover the feeders. The VVO problem was decomposed into two subproblems, at substationand feeder levels, in [69], which were then solved using a simplified dynamic program-ming and a fuzzy logic control algorithm. In [70], the VVO problem was solved using anadaptive neuro-fuzzy inference system. A multiobjective Genetic Algorithm, improved byfuzzy logic, was applied to the VVO problem in [71], in which the size of the combinatorialsearch space is reduced using expert knowledge.A Genetic Algorithm method was proposed in [72] for CVR using capacitor placement.While there are limited academic publications on the subject of CVR methodologies, theindustry is progressing gradually in implementing CVR. The key factor that governs theefficiency of VVO and CVR implementations is the voltage dependence of loads, as wasacknowledged in the report published by the US Department of Energy [73]. Preliminarytests at B.C. Hydro have shown a strong dependency of the active and reactive powerdemands on the voltage magnitude [63], which encourages the concept of voltage reduc-tion to save energy. However, some loads may need to stay “ON” for a longer period oftime to carry out the same task, compared to a normal voltage level. These loads areusually equipped with a thermostat, e.g., air conditioners, space heaters, etc. Althoughthese loads may stay “ON” for longer period as a result of CVR, the total consumed powercan be reduced at the peak and shifted to the next period of time after the peak. Theestimation of the loads’ voltage dependence has a significant impact on the performanceof the VVO algorithms.91.3 Research Objectives and Anticipated ImpactsThe overall view of the present research work is shown in Fig. 1.1. The work starts byattempting to find an accurate load model for each node in DS. To this end, a load de-composition technique is required to identify the main components of the load connectedto each node at every time interval. These two related objectives are sought in Chapter2. Based on the derived load models, a linear formulation is obtained in Chapter 3. ThisLinear Power Flow (LPF) formulation is expected to facilitate the near-real-time optimiza-tion routines required by the modern DMS. The LPF also yields very useful sensitivityfactors for DS. These sensitivity factors, also known as distribution factors, are good rep-resentatives of the effects of any change in a particular part of the network on the restof the network. Besides, a unique feature of the LPF is also taken advantage of to derivethe Linear Current Flow (LCF) equations. The LCF equations are particularly useful forthe optimization problems involving a change in the network topology. It is expected thatthe LCF formulation enhances the solution of distribution system reconfiguration problemsubstantially.Having the proper tools developed, i.e. the LPF and the LCF equations plus an accu-rate voltage-dependent load model, the second phase of the research begins. In the secondphase, the tools developed are applied to several DS optimization problems such as networkreconfiguration for loss reduction, load balancing, and voltage profile improvement, VVO,etc. In order to overcome some challenges in network reconfiguration problem, severaltechniques are developed. A DS is usually required to maintain a radial configuration.In order to impose this as a constraint in the optimization problem, a mathematical for-mulation is required. As it is shown in Chapter 4, the previously proposed formulationsare not sufficient for this purpose. A useful yet overlooked feature of a radial DS is itsrepresentation as a planar graph. This fact is utilized here to build a rigid mathematicalformulation for representing the radiality constraint in the optimization problems. In ad-dition, the solution speed might be an issue when applying the reconfiguration procedureat near-real-time. To meet this need, a very fast heuristic technique is also developedbased on the minimum spanning tree concepts in graph theory. This method can find anear-optimal solution for a network with over 10,000 nodes in less than a second.The final scope of the present study is to apply the developed tools to an importantproblem in modern DS, i.e. VVO. The urgent need for decreasing the power consumptionand network losses has pushed the utilities to come up with energy conservation tech-niques. One of the most recent advances in energy saving methods is the so-called CVR.The basis of CVR has been built on the fact that lowering the voltage can reduce the powerconsumption of load. There is more room for further investigations although the concepthas been proven by the industry. It is intended here to create a framework for a properly10Load Decomposition Load ModelingLinear Power FlowLinear Current FlowSensitivity AnalysisDistribution System OptimizationOptimal ReconfigurationVolt-VAR OptimizationFigure 1.1: Overall view of the proposed research plan.supervised VVO application to make sure that the actual outcomes will be similar to thoseplanned based on simulations.The main objectives of the present study are summarized in the following.Objective 1 Create an algorithm for load decomposition at smart meters level. This algorithmshould require minimal computation and storage requirements since it will eventu-ally be implemented on a smart meter device with limited processing and storagecapabilities.Objective 2 Build an accurate load model based on the load composition obtained by the previousstep. This load model should reflect the load voltage dependence characteristics atevery node in the system at a given moment.Objective 3 Formulate a system of linear equations to find the solution of the power flow problemin distribution systems. To this end, the voltage-dependent load model developedshould be adopted to ensure accurate power flow solutions, as compared to theconventional constant-power models.Objective 4 Derive a set of linear current flow equations based on the network synthesis obtainedby the LPF. The LCF equations are useful when the branch currents are needed tobe directly calculated, regardless of nodal voltages.Objective 5 Find the sensitivity factors for distribution systems based on the network synthesisbuilt by the LPF. These sensitivity factors are supposed to provide a fast and directway to recalculate all the system static variables, i.e. voltages and currents, upon achange in a particular entity.11Objective 6 Apply the developed analytical tools to DS optimization. In particular, optimalnetwork reconfiguration for various objectives is the main focus. Also, establish aframework for volt-VAR optimization in modern DS using the developed load modelsand the LPF.12Chapter 2Load Modeling2.1 Load Decomposition at Smart MetersThe importance of load decomposition techniques was previously discussed in Chapter1.2.1. Most of the mentioned research work on load disaggregation in Chapter 1.2.1 re-quires relatively high computational efforts in the recognition process. To this end, themeasured data needs to be sent to a central processor on which the optimization algorithmscan be run. It is desirable, however, to be able to implement the disaggregation algorithmat the smart meter level and just transfer the results using the limited communicationcapabilities of the current smart meter infrastructure. This requires a low-computationalgorithm that can be implemented using microprocessors and low-capacity memory stor-age located in the smart meter or in a box next to it. The aim here is to develop amethodology that satisfies those requirements.The main goal of this chapter is to characterize the load signature, i.e. its currentand/or voltage waveforms, in a way that is easily distinguishable. To achieve this goal,the concept of eigenloads is introduced. The eigenloads are, similar to the eigenvaluesconcept, the basic building blocks that can be used to reproduce all load signatures. Thisidea has been previously used in the well-cited paper by Turk and Pentland [74] for facerecognition. Some electric devices draw a transient current which varies over a period oftime. For example, a motor load draws a current equal to several multiples of its steadystate value during the starting period. Therefore, a time-frequency analysis is requiredto capture the full signature of a particular load. To address this, the S Transform (ST)is used here. The ST, also known as Stockwell transform, was introduced in [75] and itsadvantages over the continuous wavelet transform and short-time Fourier transform werediscussed therein. Using the ST, each current measurement is transformed into a photorepresenting the particular load, which opens the path for applying well-developed facerecognition mechanisms to the load identification problem without the need for new so-13phisticated optimization algorithms.The rest of the chapter is organized as follows. In Section 2.1.1, the load decom-position methodology is described. The application of this approach to real laboratorymeasurements is presented in Section 2.1.3. Finally, the voltage-dependent load modelingis described in Section 2.2.2.1.1 MethodologyThe first step in load decomposition is to collect measurements, i.e. the current and voltagewaveforms. The instantaneous power is chosen as the load signature. The waveform of thecurrent by itself may not give enough features to distinguish among loads since it does notaccount for the power factor (phase shift between voltage and current waveforms). Theinstantaneous power s(t) is the multiplication of current and voltage signals and, therefore,contains informations about power factor (see Section 2.1.4).Time-Frequency RepresentationDuring the starting period, each load draws a current with different frequency contents fordifferent durations. This means that a fixed-window transformation, such as the FourierTransform, is not useful to accurately represent load identities. To take this into ac-count, the author suggests the use of the family of time-frequency transformations thatprovide the frequency contents of a signal over time. There are some transformations thatgive a time-frequency representation of a signal. One of them is the Short-Time FourierTransform (STFT), which has been used frequently in many different disciplines. The dis-advantage of the STFT is that it only has a fixed resolution, which is determined by the“window function” used. For a better resolution in frequency, a wider window should bechosen whereas for a finer resolution in time a narrower window is required. However,when a good resolution is sought in both time and frequency, the STFT would not be theideal option.To overcome the mentioned shortcoming of the STFT, a multi-resolution transform isrequired. One of the commonly used techniques in this context is the Continuous WaveletTransform (CWT). The CWT of a function s(t) is given byW (t, d) = 1√d∫ +∞−∞s(τ)w∗(τ − td )dτ (2.1)where w(t) is called the “mother wavelet” and the asterisk stands for operation of complexconjugate; d and t are the dilation and translation factors, respectively. By varying tand d, the window can be shifted in time and frequency. There are many families ofmother wavelets proposed in the literature. The one that suits the application of load14decomposition needs to have good frequency resolution for low frequencies and a goodresolution in time for high frequencies. This is particularly important when the signal isthe instantaneous power because the fundamental frequency is located around 120 Hz (or100 Hz in a 50 Hz system) and a DC component is also present. Refer to Section 2.1.4 forexplanation. The ST [75] is a perfect choice for this purpose. This transformation has afrequency-dependent window size, which localizes the scalable Gaussian window dilations(d) and translations (τ). The ST of s(t) is defined as:S(t, f) =∫ +∞−∞s(τ) |f |√2pie−(t−τ)2f22 e−i2pifτ dτ (2.2)in which f is the frequency and S(t, f) is the ST of s(t). The ST has advantages over itscompetitors such as the STFT and the Wigner distribution function. It is shown in [75]that its competitors may fail to detect the high frequency contents of a signal while theST perfectly detects them. A fast version of the ST was proposed in [76] and is easy toimplement on a microprocessor.The ST of a time series can be visualized in a 2-D graph with time on the horizontalaxis, frequency on the vertical axis, and the magnitude of S(t, f) coded into colors. This2-D graph can be interpreted as an image. With this interpretation, the transient currentcorresponding to each load is mapped into a “photo” of that load. At this point, theproblem of load recognition is transformed into the familiar problem of face recognition.Principal Component AnalysisAssume a set of variables, N , which are possibly correlated. The Principal ComponentAnalysis (PCA) is a technique to develop an orthogonal basis with minimum number ofelements to represent N . Since some/all of the variables in N are correlated, it may bepossible to represent these variables using a smaller number of uncorrelated variables.Consider a set of data of n single-column vectors Li ∈ Rm. The average feature ψ of thisset can be calculated as:ψ = 1nn∑i=1Li (2.3)The difference between each vector and the average isΦi = Li − ψ (2.4)The PCA attempts to find n orthogonal vectors Vj and their corresponding eigenvalues λjthat can uniquely span the space of all Li. The eigenvectors defining this space can be15calculated from the following covariance matrix:C = 1nn∑i=1ΦiΦTi = AAT (2.5)where A = [Φ1 Φ2 . . . Φn] ∈ Rm×n. The eigenvectors of C, labeled as Vj , are the set oforthogonal basis for the data space. Some of the eigenvectors may be eliminated basedon their associated eigenvalues. That is, eigenvectors associated with small eigenvaluescan be eliminated since they contain negligible features (variances). Let all the remainingeigenvectors be placed in a matrix E = [V1 V2 . . . Vr], with r standing for the numberof remaining eigenvectors. Any Li can now be defined using a linear combination of thereduced eigenvectors:wi = ΦTi E (2.6)in which wi ∈ Rr is the vector of weights. To understand the concept, an example is givenin the following.Assume the following data is given:L1 =[12], L2 =[23.7], L3 =[35.5]Then the average isψ =[23.73]Subtracting the average from all Li yieldsΦ1 =[−1−1.73], Φ2 =[0−0.03], Φ3 =[11.77]The covariance matrix can then be calculated asC =[Φ1 Φ2 Φ3]Φ1Φ2Φ3 =[2 3.53.5 6.13]The eigenvalues and eigenvectors associated with C areV1 = κ1[−0.86830.4961], V2 = κ2[0.49610.8683], λ1 = 0.0004, λ2 = 8.1263where κ1 ∈ R and κ2 ∈ R are scaling factors and are equal to 1 to obtain a normalized(unity) Euclidean norm. The contribution of each eigenvector in each Li can be found by16−1 0 1 2 3 4 5123456L1L2L3 V2V1Figure 2.1: Illustration of PCA. V1 and V2 are scaled for proper representation.(2.6). Using these weights the vectors Li can be reconstructed as follows:L1 = 0.0084V1 − 2.0011V2 + ψL2 = 0.0000V1 − 0.0333V2 + ψL3 = 0.0081V1 + 2.0300V2 + ψThis is an exact representation of Li. However, it can be seen from the eigenvalues(λ1, λ2) that the first eigenvector does not contain significant information and, thus, canbe eliminated. By doing so, only one vector, i.e. V2, would be sufficient to reconstruct agood approximation of the original vectors. The followings are the reconstructed vectorsby only using V2:Lˆ1 =[1.0071.996], Lˆ2 =[1.9863.708], Lˆ3 =[3.0075.496]As can be seen, the error in the reconstructed vectors Lˆi is negligible.The above example has a conceivable interpretation too. The vectors Li are deliber-ately chosen close to the line y = 2x in a plane. Regardless of the number of points, onlyone basis is sufficient to represent them, since x and y are almost linearly dependent inthis case. This fact is illustrated in Fig. 2.1. Note that the length and position of the twoeigenvectors can be chosen arbitrarily, as long as the directions are maintained.EigenloadsAssume that measurement data is available for n different load signatures, each with thesame sampling frequency (fs) and time interval (tmax). This gives m = fs tmax samplesper measurement. The maximum frequency that can be detected in the signal is less than17the sampling rate. Applying the ST to the ith measurement results in a 2-D array Si(t, f)containing the magnitudes of the transformation with time and frequency coordinates. Itis possible to choose the resolution of the frequency axis to be equal to the resolution ofthe time axis. This way, the 2-D array Si(t, f) becomes an m by m square matrix.Each matrix Si is then rearranged into a single-column vector Li ∈ Rm2 . The PCA canthen be applied to the set of Li. Since A = [Φ1 Φ2 . . . Φn] ∈ Rm2×n, then the covariancematrix C is large (note that C ∈ Rm2×m2) and is computationally expensive to conductan eigenvalue analysis for. The method described in [77] avoids this problem (also usedin [74]). Instead of finding the eigenvalues and eigenvectors of C, one can find those forΛ = AT A, which is a n by n matrix. Conducting an eigenvalue analysis on a n by nmatrix is much easier than an m2 by m2 one. Assume the eigenvectors of Λ to be Uj . Theeigenvectors of C, denoted by Vj , can then be calculated as:Vj = AUj (2.7)The set of eigenvectors Vj are the orthogonal basis of the load signature space. For thisreason, we call them the eigenloads. Some of the eigenloads may be eliminated based ontheir associated eigenvalues.The eigenloads can be extracted from a library of load signatures. In this framework,load signatures form the feature space, or “load space”. The eigenloads may not resemblean individual load signature and have no physical meaning. They contain the most sig-nificant features of loads that can be used to reconstruct all the load signatures.Theoretically, the number of eigenloads to reconstruct all the signatures in the trainingset may be equal to the size of the training set. However, the number of eigenloads can bereduced to approximately reconstruct the signatures with fewer eigenloads. This reducesthe computation and storage requirements of the algorithm, but may introduce a smallerror in the recognition part. Nonetheless, this error is negligible, as is shown later. Thefollowing steps need to be taken to form a library of eigenloads:1. Collect the load signatures and form the training set T .2. Apply the ST to all the members in T .3. Calculate the eigenloads using the transformed load signatures in T .4. Keep the eigenloads with large eigenvalues (more significant) and form the library.Recognizing a New MeasurementOnce the library of eigenloads is formed, it can be used for near-real-time recognition ofinputs. The recognition algorithm is summarized in the flowchart given in Fig. 2.2a. The18StartGet the wave-forms of V (t)and I(t)Calculate X(t, f)by S transformFind Φx by (2.4)Project Φx ontoE using (2.6)Recognize theload using (2.8)End(a) RecognitionStartDetermine on/offaction using I(t)Recognizethe deviceswitched on/offUpdate the listof online devicesEnd(b) DetectionFigure 2.2: Flowchart of near-real-time load decomposition algorithm.current and voltage waveforms can be obtained from a smart meter with some samplingrate, e.g., 10 kHz. As soon as a load is turned on/off, the detection algorithm shown inFig. 2.2b captures the switching action, determines the load type using the recognitionalgorithm in Fig. 2.2a, and updates the list of ON loads. If a load was previously ON andis now turned off, the algorithm is able to detect this action based on the decrease in themagnitude of the total current. Turning off a device reduces the total current magnitude.When a device is turned on, its current waveform is measured and the triggering signal isset to zero by subtracting the present window from the previous one. A similar procedureapplies when a load is turned off.Current waveforms are monitored and when the magnitude changes beyond a certainthreshold, it is assumed that a load is turned ON/OFF. When a device is turned on, itscurrent waveform is measured and the triggering signal is set to zero by subtracting eachwindow from the previous window.Once a new signal is received, its ST is calculated (X in Fig. 2.2). Then the averagefeature ψ given by (2.3) is subtracted from X as in (2.4), denoted by Φx. The projectionof Φx into E gives the weights wx. The library of eigenloads contains the vector of weightswi for all the loads in the training set. The recognition problem here is to compare the19new vector wx with the vectors in the library and determine the most similar signature.This can be done by, for example, using normalized vectors and then compare the anglesbetween vector wx and all other vectors wi, i.e.:θi = cos−1( ~wi . ~wx‖ ~wi‖ ‖ ~wx‖)(2.8)The smallest θi is the most similar load signature to the new measurement.2.1.2 Load ClassificationMost of loads have distinct signatures, which makes the recognition easy and accurate.However, there are loads that are of different type, but may have similar signatures. Forinstance, conventional electric heating appliances (without power electronic interfaces)may have similar current waveforms. Good examples are irons, heaters, ovens, electrickettles, electric rice cookers, electric hair curlers, and even incandescent lights. Thesetypes of loads are essentially resistive loads and have negligible reactive power. They alsoproduce negligible harmonics and, therefore, have similar signatures. Although one maynot be able to detect the exact type of the load in such cases, it is still quite useful to beable to recognize the load as a member of a known family, e.g., resistive loads in this case.The above discussion motivates the author to give a more general direction to the loaddecomposition process. Instead of looking for specific loads, one may look for broaderfamilies of loads. To this end, a load classification for a typical residential customer isgiven in Table 2.1. Besides its meaningful physical interpretation, this categorizationis essentially driven by the eigenloads. When treating the load signatures as vectorsdefined by the weights obtained from (2.6), then the angles between these vectors givenby (2.8) represent the similarities between the different signatures. By carefully choosingan appropriate threshold, the loads with similar transient characteristics can be clusteredinto a distinct group. Mathematically, loads i and j belong to the same load cluster if thefollowing criterion is met:|θi − θj | ≤ θth (2.9)A proper value for θth can be obtained by try-and-error. In this case, experience showedthat θth = 10◦ is a viable choice to avoid overlapping between the different subclasses.Based on this load classification, instead of categorizing a measurement as a specific load,the load identification algorithm categorizes the measurement under a broader family ofloads. This is also a relaxation on the load identification algorithm and, thus, greater ratesof success in the recognitions are achieved.The proposed load decomposition algorithm is capable of identifying the loads undereach subclass given in Table 2.1. Further distinguishing between loads falling under the20Table 2.1: Load Classification for a Typical Residential CustomerClass No. SubclassMotor1 Small IM2 Medium IM3 Power electronic-driven IMResistive 4 Direct interface5 Power electronic interfaceFluorescent lamp 6 Magnetic ballast7 Electronic ballastTV 8 CRT9 LCD-LEDElectronic Cooker 10 Microwave oven11 Induction cookerElectronics 12 Chargers13 OthersClass No. Features ExamplesMotor1 PF lag Fan2 Inrush current, PF lag Vacuum cleaner, AC3 Non-sinusoidal Wash. machine, hand mixerResistive 4 Sinusoidal, unity PF Incand. light, heater5 Harmonics Oven, ironFluorescent lamp 6 Multi-stage start, PF lag -7 Multi-stage start, PF lead,harmonics -TV 8 Harmonics -9 Harmonics, PF lead -Electronic Cooker 10 Multi-stage, harmonics -11 Multi-stage, harmonics -Electronics 12 Multi-stage, harmonics Laptop/phone charger13 Harmonics PC, DVD-playersame subclass is not usually useful from the power utility point of view. For instance, ameat grinder and a vacuum cleaner may have similar signatures which makes them fallunder the second subclass, i.e. medium-size Induction Motor (IM). From a load modelingperspective, however, it suffices to know they belong to the class of IM rather than theclass of other appliances.2.1.3 Laboratory MeasurementsIn this section, the application of the proposed framework is shown through real measure-ments. A number of electric devices were analyzed in the lab and their current and voltagewaveforms were recorded for a period of a few seconds. The sampling frequency is 10 kHz.With this sampling rate, harmonics up to 1.0 kHz can be detected with high accuracy.A current probe and a voltage probe (differential probe) for measuring the correspondingquantities were attached to an oscilloscope. Furthermore, the measurements collected atthe Carnegie Mellon University [78], which are freely available to the public under thename BLUED, were also used. In total, 46 measurements were available, from which 2321were chosen for building the eigenloads and the rest were left for testing.Three examples of measurements and their ST are shown in Figs. 2.3-2.5. The fig-ures illustrating the ST results are actually the contour visualization of square matriceswith normalized elements. Some loads, such as the vacuum cleaner (induction motor)shown in Fig. 2.3, exhibit a fairly constant frequency content over time, while the magni-tude of the current has a decaying component (inrush current). The instantaneous powerof induction motors has a large DC component and a 120 Hz component decaying overtime. Some loads have multiple stages from the starting moment to the steady state.The current waveform from a MAC laptop charger shown in Fig. 2.4a, for example, hasthis characteristic. It shows higher frequency components in the middle stage and theystart damping out after a few cycles. Some power electronic-driven loads draw a currentrich in harmonics of different frequencies. The current waveform of a fluorescent lampwith electronic ballast, shown in Fig. 2.5a, has a low frequency envelope on top of thefundamental frequency. This feature produces a variety of frequency components in theinstantaneous power, which are reflected in the ST shown in Fig. 2.5b. The reason behindthe existence of the even harmonics in the ST is briefly investigated in Section 2.1.4.The 23 loads used to construct the eigenloads include incandescent lamps, fluorescentlamps (with magnetic and electronic ballasts), hair drier, refrigerator, fan, a variety ofinduction motors, blender, meat grinder, vacuum cleaner, hand mixer, phone chargers,laptop chargers, electric shaver, personal computer, TV (LCD, LED, CTR), DVD player,iron, heater, induction cooker, hot plate, and microwaves. These loads cover most of thetypical loads for a residential customer. Further analysis is to be done for commercial andindustrial loads. The idea, however, is applicable to all type of consumers.Evaluation of the Load Decomposition MechanismBased on the eigenvalue criterion, 10 eigenloads out of 23 were discarded, i.e. only 13eigenloads are sufficient to form a basis for the load signatures space. This reduces thememory storage and processing requirements substantially. This fact is further discussedin Section 2.1.4. In order to test the proposed load decomposition approach, measurementsfor devices that were not present in the training set were carried out. The results of sixtests are discussed here for illustration purposes. An LCD load manufactured by Samsungwas used in the training set, while an LCD manufactured by ViewSonic with differentratings was used to test the algorithm. The results for this case are shown in Fig. 2.6a.A hand mixer equipped with multiple stages, which was not used in the training set, wasused at its lowest and highest power settings and the recognition algorithm was applied.The results are shown in Figs. 2.6b and 2.6c. It should be noted that in its first powersetting, the hand mixer is driven by power electronics and, hence, the current drawn fromthe network is not sinusoidal (Subclass 3). In its second setting, however, the device draws220 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9−40−30−20−10010203040Time (s)Current (A)(a) Instantaneous current (I(t))Time (s)Frequency (Hz)0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10120240360480600720840(b) S Transform (S(t, f))Figure 2.3: Current waveform of a vacuum cleaner and the ST of its instantaneouspower.a fairly sinusoidal current and, thus, resembles an induction motor (Subclass 2).A heater equipped with a small fan, also not present in the training set, was tested andthe results are shown in Fig. 2.7a (Subclass 4). A compact fluorescent lamp was testedwhich has slightly different waveforms than the T5 fluorescent lamps used in the training230 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9−0.3−0.2−0.100.10.2Time (s)Current (A)(a) Instantaneous current (s(t))Time (s)Frequency (Hz)0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10120240360480600720840(b) S Transform (S(t, f))Figure 2.4: Current waveform of a MAC laptop charger and the ST of its instanta-neous power.set. The recognition algorithm, as shown in Fig. 2.7b, was still able to appropriatelydetermine the closest load type (Subclass 7). A personal computer was also tested andthe results are shown in Fig. 2.7c, in which the algorithm categorizes this load underSubclass 13. The proposed algorithm also demonstrated a high rate of success in other240 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1−1.5−1−0.500.511.5Time (s)Current (A)(a) Instantaneous current (I(t))Time (s)Frequency (Hz)0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10120240360480600720840(b) S Transform (S(t, f))Figure 2.5: Current waveform of a T5 fluorescent lamp with electronic ballast andthe ST of its instantaneous power.laboratory tests, not presented here.There is a trade-off between the number of eigenloads used and the accuracy of therecognition algorithm. Increasing the number of discarded eigenloads can save on mem-ory storage and computation requirements, while potentially affecting the accuracy of the251 2 3 4 5 6 7 8 9 10 11 12 13050100150Load Class No.θ(deg.)(a) LCD (32 W)1 2 3 4 5 6 7 8 9 10 11 12 13050100150Load Class No.θ(deg.)(b) Electric hand mixer-Stage 1 (300 W)1 2 3 4 5 6 7 8 9 10 11 12 13050100150Load Class No.θ(deg.)(c) Electric hand mixer-Stage 5 (300 W)Figure 2.6: Results of recognition tests on devices not present in the training set (θis defined by (2.8)).recognition process. For example, let us reconstruct the vacuum cleaner’s ST using a re-duced number of eigenloads. Figure 2.8 shows the changes in the ST when the number ofeigenloads used (r) is continuously reduced. As can be seen, down to r = 10, the image isstill close to the original one (r = 23). Beyond this point, however, the differences becomenoticeable. A good measure to compare the reconstructed image and the original image isthe angle difference given by (2.8). This is done by setting the weights corresponding to261 2 3 4 5 6 7 8 9 10 11 12 13050100150Load Class No.θ(deg.)(a) Heater (1 kW)1 2 3 4 5 6 7 8 9 10 11 12 13050100150Load Class No.θ(deg.)(b) Compact Fluorescent Lamp (15 W)1 2 3 4 5 6 7 8 9 10 11 12 13050100150Load Class No.θ(deg.)(c) Personal Computer (500 W)Figure 2.7: Results of recognition tests on devices not present in the training set (θis defined by (2.8)).the discarded eigenloads to zero in ωi. This procedure was carried out for all the 23 loadsand the results are shown in Fig. 2.9. As can be seen, by discarding up to 10 eigenloads,the difference between the original and reconstructed images remain within 10◦. There-fore, 13 eigenloads are found to be sufficient to define the load signature space.The number of eigenloads used can affect the success rate in the recognition process. Inorder to show this, eignloads were discarded one after another and the recognition proce-27 r = 23 r = 1 8 r = 10 r = 8 r = 6 r = 2 Figure 2.8: The S Transform of the vacuum cleaner reconstructed by a reducednumber of eigenloads (r is the number of used eigenloads).0 2 4 6 8 10 12 14 16 18 20 220102030405060708090Number of Discarded Eigenloadsθ(◦)Figure 2.9: Differences between the original and reconstructed images with reducednumber of eigenloads.dure was applied to all the test cases. The rate of success versus the number of eigenloadsdiscarded is illustrated in Fig. 2.10. As can be seen, discarding up to 11 eigenloads did notaffect the rate of success in the test cases. However, the rate of success is slightly droppedwhen more eigenloads are eliminated. By discarding 18 or more eigenloads, the rate ofsuccess is significantly deteriorated. Therefore, discarding 10 eigenloads keeps a balancebetween the accuracy and the computation and memory requirements of the proposedalgorithm.280 2 4 6 8 10 12 14 16 18 20 22020406080100Number of Discarded EigenloadsSuccessRate(%)Figure 2.10: Rate of success in the recognition process versus the number of dis-carded eigenloads.Robustness Against Noisy MeasurementsIn order to show the robustness of the proposed approach, the effect of noise on its per-formance is analyzed here. A random Gaussian noise with a zero-mean and a varianceequal to 30% of the RMS value of current was added to the original current waveforms.For example, Fig. 2.11 shows the original current waveform and the noisy version for aLCD unit. The noisy measurements were used to perform the recognition process. Theresults are shown in Fig. 2.12a. The same process was done for the hand mixer withthe results shown in Fig. 2.12b. These results show that the accuracy of the proposedapproach is robust under noisy measurements. The main reason is that the S transformfilters the higher frequencies since the upper bound on the frequency components used inthe training and recognition process is set to the fifteenth harmonic (900Hz).2.1.4 DiscussionAnalytical Derivation of Instantaneous PowerThe current waveform is a good choice for the load signature since each type of load, tosome extent, has a distinct current waveform. However, as discussed earlier, the phase shiftbetween the voltage and current waveforms (i.e., the power factor) would be lost if onlythe current waveform were used. For this reason, instantaneous power (s), which is thepoint-wise multiplication of voltage and current signals, is used here. Assume a sinusoidalwaveform of the voltage and a non-sinusoidal waveform of the current, represented by its29TimeCurrentFigure 2.11: Current waveforms of a LCD; bottom: original; top: noisy by a randomGaussian noise with a zero-mean and 30% of RMS value of current as itsvariance.1 2 3 4 5 6 7 8 9 10 11 12 13050100150θ(deg.)0% 30%(a) LCD (32 W)1 2 3 4 5 6 7 8 9 10 11 12 13050100150Load Class No.θ(deg.)0% 30%(b) Electric hand mixer-stage 5 (300 W)Figure 2.12: Effect of noise (as in Fig. 2.11) on recognition accuracy.truncated Fourier series:V (t) = V1 sin(ω t) (2.10a)30I(t) = I0 +N∑k=1Ik sin(k ω t+ φk) (2.10b)With these assumptions, s(t) can be calculated as:s(t) = I0V1 sin(ωt) +N∑k=1V1 Ik2[cos((k − 1)ωt+ φk)− cos((k + 1)ωt+ φk)](2.11)Consider, for example, that the current has only the first and the third harmonics (e.g.,induction motors). In this case, s(t) would have a DC component as well as the secondand the fourth harmonics of the fundamental frequency. This simple analysis explainsthe presence of the even harmonics in the S transforms. In some cases, e.g., some powerelectronic-driven loads, the current drawn by the device may have a DC component orsome even harmonics. Under such circumstances, s(t) has a frequency content of bothodd and even multiples of the grid frequency.ApplicationsIn terms of applications, load composition at the customer level is useful in the followingaspects.1- Tariff PoliciesMulti-rate tariff for energy consumption has been used as a tool for demand-side man-agement and/or energy saving purposes. This system of energy tariff is usually based ontime-of-use or level of energy consumption. In the first case, the $/kWh is different indefined time intervals, usually more expensive during the peak time and cheaper in lightload periods. In the second case, there are thresholds on the kWh/month and when eachone of those thresholds is exceeded, the $/kWh increases. In order to design an efficienttariff system in both cases, an understanding of the customer behavior is useful. If thecomponents of load used by the customer are known to the utility, it is then possible toestimate how much of this load is potentially responsive to a price signal or other similarpolicies. For instance, if the peak load occurs around 7 pm, and the measurements showthat many households use the washing/drying machines at this time, then there is a pos-sibility to shift part of the load by intensifying the customers. It is then a good measureto estimate the effectiveness of a tariff policy before applying it. More discussion on thistopic can be found in, e.g., the EPRI report on customer behavior [79].312- Customer BenefitsBy carefully studying the customer’s electricity consumption, it is possible to provideuseful feedback on possible ways to conserve energy. For example, by monitoring theconsumption of each appliance, inefficient appliances can be identified. A cost-benefitanalysis can be provided on replacing the inefficient appliance with a better one. Also,suggestions can be given on using particular devices during light-load periods to save onelectricity bill.3- Short-Term Load ForecastingThe average load is sensitive to several factors including ambient temperature, sun light,grid voltage and frequency, etc. There are some specific components of load that are moresensitive to the mentioned factors than others. For example, if temperature increasesabove the expected values, the air conditioners need to work more than expected. If theportion of load corresponding to air conditioning and the sensitivity of this type of loadto temperature are known, then a good estimation of load deviation from the forecast canbe derived. The load sensitivity to temperature is referred to as temperature elasticity in[80], in which the effects of weather condition on load is thoroughly discussed. Moreover,the load composition estimation can be used in building-level demand response in thedistribution system market [81].4- Near-Real-Time Load ModelingThe load voltage dependence has an important impact on the load power consumption.It is shown in the B.C. Hydro system that decreasing the substation voltage by 1%, theactive and reactive demand decrease by 1.5% and 3.4%, respectively [82]. Different typesof loads have been studied to assess their voltage dependence [83]. In distribution systemlevel, near-real-time control actions may be applied by the DMS. These actions includevolt-VAR control and optimization (VVO) and voltage-dependent power flow analysis [84],which require a good estimation of load response to voltage variations. A good voltage-dependent load model is beneficial for such cases. In transmission system level, near-real-time dynamic analysis such as transient and voltage stability assessments require a goodestimation of load model at each substation. The aggregated load composition at thesubstation can be used to build an accurate model of load static and dynamic behaviorseen from the transmission system [85].Implementation AspectsAdditional values can be added to the present advanced metering infrastructure by im-plementing the proposed load decomposition technique at the smart meters level. Load32composition can be sent along with the electricity consumption and voltage level to thecontrol center, where the received data is stored and analyzed. In order to implement theproposed algorithm at the smart meters, the memory storage and computation require-ments of the algorithm should not exceed the capabilities of available hardware.The largest memory requirement of the algorithm is pertained to the library of eigen-loads. It was shown in Section 2.1.3 that 13 eigenloads are sufficient to have satisfactoryresults. Also, the resolution of the images representing each eigenload can be reducedto save storage requirements, without affecting the accuracy of the results. Based onexperience, a resolution of 500 × 500 is sufficient. Using this resolution and assuminga double-precision floating-point format, each eigenload requires 500 × 500 × 8 bytes ofmemory storage. For 13 eigenloads, this adds up to 26 Mb. This is well within the capa-bility of external Random Access Memory (RAM) available for Field-Programmable GateArray (FPGA). An example is Stratix V from Altera that works with DDR3 memories at1,866 Mbps.The highest computation-intense part of the algorithm is the calculation of the STransform (ST) of a measurement. The discrete ST proposed in [76] has a computationalcomplexity of O(N logN), where N is the number of data points in the signal. This isequivalent to the complexity of the Fast Fourier Transform [76], which has been previouslyimplemented in many hardware platforms. The structure of the ST is in such a way thateach point in the time-frequency plane can be calculated independently. Therefore, effi-cient parallel implementation of the ST can be realized on common microprocessors suchas FPGA. A good example of real implementation of the ST on FPGA is [86].LimitationsSince the proposed method works based on the transient features of loads, it does not workproperly when two loads of different types are switched on simultaneously. This occursbecause the algorithm assumes the captured transient belongs to one single type of loadand, therefore, its time-frequency representation stands for a particular load. Nonethe-less, this is a rare event that two loads of different types are turned on simultaneously ina household.Another limitation of the proposed algorithm is that it is unable to further distinguishbetween loads with similar signatures but different applications. For example, most of theresistive loads have similar signatures, while they may have different functions. For thetargeted applications mentioned in Section 2.1.4, the resolution in the load compositionprovided by Table 2.1 is sufficient.The proposed method requires measurements with relatively high sampling rates. Al-though the present data available from the smart meters has usually low sampling rates,the internal process in the smart meter itself occurs at much higher sampling rates. There33is a great potential for implementing the proposed method at the smart meter and trans-mitting the load composition along with power consumption.2.2 Voltage-Dependent Load Modeling2.2.1 Conventional Load ModelsPower system loads show different reactions to grid voltage variations. A common wayof describing the dependence of active and reactive power consumption on the voltagemagnitude is the exponential model [87], as follows:P (V )P0=( VV0)α(2.12)Q(V )Q0=( VV0)β(2.13)where P and Q are the load’s active and reactive power consumption; V is the terminalvoltage magnitude and the zero subscript means the nominal value; α and β are the activeand reactive power exponents, respectively, which can be extracted from measurements.Some typical values for these exponents are given in [88] and [89].In addition to the exponential load model, the polynomial model (ZIP model) has alsobeen widely used in power system studies. This model consists of three major parts: aconstant-impedance (Z), a constant-current (I) and a constant-power (P). Mathematically,this model describes the load variation with voltage as [87]:P (V )P0= FZ( VV0)2+ FI( VV0)+ FP (2.14)Q(V )Q0= F ′Z( VV0)2+ F ′I( VV0)+ F ′P (2.15)where F and F ′ are fitting parameters; the subscripts Z, I and P stand for constant-impedance, constant-current and constant-power contributions, respectively. Note thatthere are only two independent parameters in (2.14) and (2.15) because FZ +FI +FP = 1and F ′Z + F ′I + F ′P = 1.2.2.2 Proposed Load ModelThe constant P-Q load model, as well as the voltage-dependent load models of (2.12)-(2.15) introduce nonlinearity in the solution of power flow equations. The load modelproposed here is an alternative to the voltage-dependent load models of (2.12)-(2.15) that34allows for a linear formulation of the power flow equations (this modeling approach willbe referred to as LPF load modeling from now on). The proposed model is described as:P (V )P0= CZ( VV0)2+ CI( VV0)(2.16)Q(V )Q0= C ′Z( VV0)2+ C ′I( VV0)(2.17)in which constants C and C ′ are calculated by a curve-fitting procedure. Note that there isonly one independent parameter in (2.16) and (2.17) because CZ+CI = 1 and C ′Z+C ′I = 1.From a mathematical synthesis point of view, exponents α and β in the model given by(2.12)-(2.13) can be calculated with a fitting procedure, while in the model described by(2.14)-(2.15) the polynomial coefficients are to be determined. The proposed synthesis of(2.16)-(2.17), when compared to the other two, introduces computational advantages informulation of power flow equations. This fact is discussed in Chapter 3.The proposed model approximates the load voltage dependence using an impedance Z,representing the quadratic voltage dependence, and a current source I, representing thelinear voltage dependence. The coefficients CZ , CI , C ′Z , C ′I can be found for a very goodfit to the measured data. The fitting process can be formulated in terms of a simple convexquadratic optimization problem. The fitting objective for the active and reactive poweris to minimize the difference between the fitted approximation and the measured data fora finite number of voltage points within a range of operating voltages (e.g., ±10%). Forexample, for P in p.u. values (writing in p.u. allows for dropping V0 and P0),MinimizeNv∑i=1[CZV 2i + CIVi − P (Vi)]2 (2.18a)subject toCZ + CI = 1 (2.18b)and similarly for reactive power Q. In (2.18), Nv is the number of points selected withinthe voltage range;(Vi, P (Vi))is the ith pair of measured (voltage, active power); CZ andCI are the unknowns. Note that (2.18b) reduces the number of basic variables to one.The solution of this problem is straightforward, as follows.In order to find the optimum solution of (2.18), the affine constraint (2.18b) is used todrop one dependent variable, say CI . Assume that the voltage and active power measure-ments are in per-unit with respect to their nominal values, i.e. V0 and P0. The reducedproblem can be written as35Minimize f(CZ) =Nv∑i=1[CZ(V 2i − Vi)− P (Vi) + Vi]2 (2.19)Taking the derivative of f with respect to CZ and putting it equal to zero (Karush-Kuhn-Tucker conditions for optimality), one has:CZ =Nv∑i=1(V 2i − Vi)(P (Vi)− Vi)Nv∑i=1(V 2i − Vi)2(2.20)Based on the fact that the objective function is convex quadratic, this is the optimalsolution. The solution for CI can be retrieved using (2.18b).The data for the synthesis in (2.18) can be obtained directly from load measurementsor previous synthesis by (2.12)-(2.15). The ZIP model, exponential model, and the LPFload model obtained by fitting a curve to experimental measurements for a 3-phase in-duction motor (measured in the lab) are depicted in Fig.2.13. Figure 2.13 also showsthe parameters for the fitted curves. As can be seen, the difference between the LPF andthe exponential model is negligible. Also, the ZIP model provides slightly tighter ap-proximation compared to the LPF and exponential models. For other types of loads thatshow stronger dependence on voltage, even tighter fits have been achieved using the LPFmethod. The parameters for some commonly-used types of loads are calculated throughmeasurements and are given in Table 2.2. It should be noted that since the ZIP modelhas two independent variables (since FZ + FI + FP = 1) for the fitting while the pro-posed LPF load model and the exponential model have only one, more accurate fittingsmay be achieved by the ZIP model. However, the differences between the resulting curvesare minor. Depending on the load type, it sometimes happens that negative values areobtained for F in (2.14)-(2.15), or C in (2.16)-(2.17). These negative values do not affectthe mathematical solutions, although they do not have a physical meaning.Table 2.2: Load Parameters in (2.16) and (2.17) Obtained by MeasurementsLoadTypeCompactFluorescentLampIncandescentLampElectricStovePersonalComputerCZ -0.394 0.578 0.967 -0.984C ′Z -1.857 -5.012 0.960 -1.351360.9 0.92 0.94 0.96 0.98 1 1.02 1.04 1.06 1.08 1.10.940.960.9811.021.041.06α = 0.5598FZ = 0.0255,FI = 0.5090,FP = 0.4655CZ = −0.4365,CI = 1.4365V (p.u.)P(p.u.)Measurements Exponential ZIP LPF(a) Active Power0.9 0.92 0.94 0.96 0.98 1 1.02 1.04 1.06 1.08 1.10.70.80.911.11.21.31.4β = 2.7866F ′Z = 5.0936,F ′I = −7.4387,F ′P = 3.3451C′Z = 1.7738,C′I = −0.7738V (p.u.)Q(p.u.)Measurements Exponential ZIP LPF(b) Reactive PowerFigure 2.13: Comparison of load models for a 3-phase induction motor (460 V, 3-phase, 1.4 HP, 1725 RPM).37Chapter 3Linear Power Flow Analysis3.1 The Linear Power Flow EquationsDistribution systems are essentially unbalanced three-phase systems. However, in mostof cases, the major control equipment used to adjust the voltages (transformers), reactivepower injections (capacitors), and network topology (switches) are usually three-phaseunits. In other words, all three phases operate simultaneously. Taking this fact intoaccount, it is justifiable to assume a balanced system for the control and automationpurposes at feeder levels. In studies directly addressing the phase unbalance problem, onthe other hand, a three-phase unbalanced load flow algorithm is required.Conventionally, the PF equations have been formulated based on a constant P-Q loadmodel. This makes the equations nonlinear and iterative methods are required to find thesolution. It is possible, however, with the load model proposed in Chapter 2, to reformulatethe PF equations into a linear system of equations. A similar attempt was made in [33] byreplacing the loads with current source injections. In this method (implicit Z-bus method),however, an iterative procedure is still required to update the injected currents at everyiteration. The updating mechanism for current injections at a generic node at iteration kisI¯(k) =( S¯(k)V¯ (k))∗(3.1)where values with a bar on top are complex numbers. The LPF load model suggestedin Chapter 2 introduces the modification that S¯ is not constant but voltage-dependent.Considering in (3.1) that P and Q are voltage-dependent and separating the real andimaginary parts of S¯ and V¯ , one can writeI¯(k) = P(k)(V (k))− jQ(k)(V (k))V (k)re − jV (k)im(3.2)38where Vre and Vim are the real and imaginary parts of voltage, respectively. Substitutingthe values of P and Q from (2.16)-(2.17) into (3.2), and temporarily dropping the iterationindex yields:I¯ = P0CZVre +Q0C′ZVimV 20︸ ︷︷ ︸Impedance+ P0CIVre +Q0C′IVimV0V︸ ︷︷ ︸Current+ j[ P0CZVim −Q0C ′ZVreV 20︸ ︷︷ ︸Impedance+ P0CIVim −Q0C′IVreV0V︸ ︷︷ ︸Current](3.3)In distribution systems, taking the voltage angle of the feeding substation as reference(zero value), the imaginary part of voltages Vim is smaller than the real part Vre by severalorders of magnitude, as is also assumed in [27]. This feature allows for eliminating theimaginary part of voltages in the current parts of (3.3). There are two common nonlinearterms in the current parts of (3.3), which are reproduced here:f1 =VreV =Vre√V 2re + V 2im(3.4a)f2 =VimV =Vim√V 2re + V 2im(3.4b)For the normal range of variables, the value of f1 is always close to 1, whereas the valueof f2 is close to zero. Since the voltage angles are relatively small in distribution systems(see Section 3.1.1 for explanation), the following approximations are then in order:f1 ≈ 1, f2 ≈ 0 (3.5)With these approximations, (3.3) can be simplified to its real and imaginary parts asIre = <{I¯} ≈Q0C ′ZV 20Vim +P0CZV 20Vre +P0CIV0(3.6a)Iim = ={I¯} ≈P0CZV 20Vim −Q0C ′ZV 20Vre −Q0C ′IV0(3.6b)As it is shown next, these approximations lead to a Norton equivalent of load synthesis.This load synthesis is then used to formulate the linear current flow equations. From nowon, the formulation obtained using (3.6) is referred to as the LPF.393.1.1 Voltage Variation in a Distribution FeederThe intention of this subsection is to show the variation of voltage magnitude and anglein a distribution feeder. Consider the simple 2-node distribution feeder shown in Fig. 3.1.Typical parameters for a 25 kV overhead line are R = 0.25Ω/km and X = 0.39Ω/km. Thelength of a feeder is limited by the maximum allowed voltage drop at the end of the feederfor the peak load (normally 6% drop is the maximum allowed according to the CanadianStandard Association, see Table 3.2). Assuming a load of 5 MW with a power factor of99%, the maximum feeder length would be 24 km. The voltage magnitude and angle forthis simple case can be obtained as follows:V1 =√(V2 +RP +X Q)2 + (X P −RQ)2 (3.7a)δ = tan−1( X P −RQV2 +RP +X Q)(3.7b)Using the typical parameters mentioned above, a phasor diagram is drawn in Fig. 3.2.This is one of the worst cases that could happen in a distribution feeder that the totalload is connected to the end of a long feeder. In this case, the maximum voltage angle is3.4◦ and the maximum voltage drop is 6%. In transmission systems, the X/R is relativelylarge, compared to distribution systems, and thus large voltage angles appear. This simpleanalysis shows the justification of assuming small voltage angles in distribution systems.I Rj XP + j QV1 6 δ V2 6 0Figure 3.1: A simple 2-node distribution feeder.~V1~V2~IR~IjX~IFigure 3.2: Phasor diagram for voltages in the 2-node distribution feeder.403.1.2 A Norton Equivalent of Load SynthesisMany software packages for system optimization are available that work only with realvariables. To use these packages, it is required to keep the real and imaginary parts ofthe complex voltages and currents separate. Writing the current drawn by a constant-impedance load as IL = YLV and separating real and imaginary parts of V and YL yieldsIL = YLV = (GL + jBL)(Vre + jVim) = (GLVre −BLVim) + j(BLVre +GLVim) (3.8)It can be seen from (3.8) that the voltage-dependent part of the drawn current given by(3.6) can be synthesized by an equivalent admittance with the following parameters:GL =P0CZV 20, BL = −Q0C ′ZV 20(3.9)while the rest of the terms in (3.6) can be represented as a constant-current source. Eventu-ally, the proposed load model for the LPF analysis is a parallel combination of an impedanceand a current source. The LPF load model can then be represented by circuit elements asshown in Fig. 3.3. The real and imaginary parts of the constant-current sources are givenby the last terms in (3.6), respectively, i.e.:Ip =P0CIV0, Iq = −Q0C ′IV0(3.10)3.1.3 The Linear Power Flow FormulationApplying Kirchhoff’s Current Law to a generic distribution system, considering the sub-station(s) as voltage source(s) and representing the loads using the LPF load model yields,YLi ILiYkYLjILjVi VjFigure 3.3: A generic part of distribution system: Norton equivalent.41in complex numbers form [YAA YABYBA YBB][VAVB]=[IAIB](3.11)in which Y , partitioned into four sub-matrices, is the system admittance matrix whichincludes the equivalent admittances of the loads given by (3.9) in the diagonal elements;VA is the vector of known voltages and IA is the unknown vector of corresponding nodalcurrent injections at the substation(s). IB is the constant-current part of the load givenby (3.10). VB is the vector of unknown voltages, which can be computed as:[VB] = [YBB]−1[IB]− [YBB]−1[YBA][VA] (3.12)Separating the real and imaginary parts of (3.11), the proposed LPF equations inrectangular coordinates can be represented by the following partitioned matrix equation:[G1,1:n −B1,1:nB1,1:n G1,1:n][G2,1:n −B2,1:nB2,1:n G2,1:n]...[Gn,1:n −Bn,1:nBn,1:n Gn,1:n]Vre,1Vre,2...Vre,nVim,1Vim,2...Vim,n=Ip,1Iq,1Ip,2Iq,2...Ip,nIq,n(3.13)in which G and B are the real and imaginary parts of the modified admittance matrix Y ,respectively. The same matrix reduction technique used in (3.11) and (3.12) can be usedhere to eliminate the known variables. By rearranging the rows in (3.13), a more compactform can be obtained as [G −BB G][VreVim]=[IpIq](3.14)The system of linear equations in (3.14) yields the solution for nodal voltages.3.1.4 Distributed Generation Modeling within the LPFThe number of distributed energy resources installed in distribution systems is constantlyincreasing. In most utilities, the grid code does not require distributed generations (DG)to act as P-V nodes. For small-scale DG, e.g., rooftop solar panels, there is no voltagesupport in terms of reactive power to the grid. These factors make it possible to modelDG as constant-power negative loads. The aim here is to show how to model the constant-power loads using the voltage-dependent load model introduced in Chapter 2.42A constant-power load can be defined as a load which its power consumption doesnot depend on the terminal voltage. Using a least-square method, the values reported inTable 3.1 are obtained for different voltage intervals. As can be seen, the error of thisapproximation grows when the voltage range is expanded. Nonetheless, it is declared thatthe application of the LPF method is intended for the system normal operation, for whichthe voltages are usually kept within ±10. Using the given parameters for this standardrange, exact power flow solutions based on the constant-power load model are comparedto the LPF solutions for several test systems. In all the cases, the solution error is lessthan 0.1%. More details in error analysis of the LPF solution is given in the next sections.3.1.5 Multi-Region Linear Power FlowThe LPF performs very well within the normal range of voltage variations, namely ±10%.However, if the voltages fall outside this range, the error of the LPF solution may increase.The author understands that this is an unusual case for the voltages to be outside the±10%range. According to the Canadian Standard Association, the standard voltage ranges fordistribution systems (below 50 kV) are defined for normal and extreme conditions in [90].Table 3.2 shows these voltage levels. According to this standard, 95% of the time theutility will not deviate from the normal operating range and that 99.9% of the time theutility will not deviate from the extreme operating range. Based on this standard, theapplication of the LPF method is well justified.For the sake of completeness, a Multi-Region Linear Power Flow (MR-LPF) is proposedhere to cover for the cases with abnormal voltages. In the MR-LPF, voltage levels aredivided into multiple regions. The first region is defined as D1 = {0.95 ≤ V ≤ 1.05}. Inthe first iteration, the solution is found by load equivalents calculated in D1. If there is anynodal voltages outside D1, the parameters representing the load for that node are updatedaccording to Fig. 3.4, depending on which region they fall into. In this figure, the loadsare modeled as constant-power and then approximated by the proposed voltage-dependentload model. Therefore, the MR-LPF is a two-iteration process and each iteration requiresTable 3.1: Error of Approximating a Constant-Power Load Model by (2.16)Voltage Range (p.u.) CZ CI Maximum Error (%)[0.95, 1.05] -0.998 1.999 0.16[0.90, 1.10] -0.997 2.000 0.68[0.80, 1.20] -0.994 2.008 2.86[0.70, 1.30] -0.993 2.022 6.71[0.60, 1.40] -0.990 2.041 12.50430.65 0.75 0.85 0.95 1.05 1.150.9940.9960.9981.0001.002V (p.u.)P(p.u.)D0[−0.8251.817]ξ+ = 0.04ξ− = 0.17D1[−12]ξ+ = 0.10ξ− = 0.16D2[−1.2322.221]ξ+ = 0.10ξ− = 0.23D3[−1.5592.499]ξ+ = 0.14ξ− = 0.27D4[−2.0362.856]ξ+ = 0.16ξ− = 0.38Figure 3.4: Voltage regions, divided by vertical dashed-lines, and the correspond-ing approximate constant-power loads (parabolas). CZ and CI are shown insingle-column matrices for each region. ξ+ and ξ− are the relative approxi-mation errors (%) above and below 1 p.u., respectively.solving a sparse system of linear equations.The MR-LPF can also be used in optimization routines, whenever needed, in the sameway as explained above. An optimization problem is first solved, assuming all the voltagesfall in D1. If the solution of the first optimization problem shows any voltages in regionsother than D1, then the corresponding parameters are updated and the problem is solvedfor the second time.3.2 Simulation Results-LPFIn this section, the LPF formulation is applied to sample test systems with 14, 70, and 135nodes, described in [46], [91], and [92], respectively. A larger system consisting of 3249nodes was also tested by combining the 135-node and 70-node systems, 16 times each.These test systems are described in Table 3.3. For simplicity, the voltage dependence ofall the loads is assumed to be identical. Also, both active and reactive power are consideredto have the same voltage dependence.A comparison between the results obtained using the implicit Z-bus method of [33],Table 3.2: Standard Voltage Ranges for Distribution Systems [90]Vmin(%) Vmax(%)Normal Operating Range -8 +4Extreme Operating Range -12 +644Table 3.3: Dimensions of the Employed Test SystemsTest Case Nodes Branches Feeders Load(MVA)14-node 14 13 3 28.70 + i17.3070-node 70 69 4 4.47 + i3.06135-node 135 134 8 18.31 + i7.933249-node 3249 3248 192 364.50 + i175.86assuming a voltage-dependent load model, and the proposed method was carried out. Theperformance of the proposed method was evaluated based on the relative error betweenthe voltage magnitudes calculated by the LPF method and the implicit Z-bus method as:η = 1nn∑i=1|V ′i − Vi|V ′i(3.15)where V ′ is obtained by the exact method and V is obtained by the LPF method. Thehistogram of errors for the 70-node, 135-node, and the 3249-node test systems are shownin Fig. 3.5. For the 70-node system, the maximum relative error is about 0.006% whichhappens when CZ = 0. For other values of CZ the error is always smaller. For othersystems, the relative error does not exceed 0.35% for the considered values of CZ . Figure3.6 shows the variation of η as a function of CZ for different test systems (C ′Z is assumedequal to CZ for the test). As can be seen in Fig. 3.6, the errors are very small (lessthat 0.11%) over the range of values considered for CZ . In the special case of CZ =1, i.e. constant-impedance (no current source in the model), the error is zero. It isworthwhile mentioning that, depending on the value of CZ , the Z-bus method needs atleast four iterations to reach a tolerance of 1e-4 p.u. in the magnitude of all nodal voltages.Increasing the value of CZ , the number of required iterations increases. This result isconsistent with the results found in [93]. In terms of iterations, the LPF model requiresonly one solution and it is, therefore, at least four times faster. A comparison betweencomputation times and number of iterations is given in Table 3.4.Figure 3.7 shows the voltage magnitudes and angles for the 135-node system, usingCZ = C ′Z = 0 which causes the maximum error. The figure compares the results obtainedusing the Z-bus method and the LPF method. The errors in voltage magnitudes does notexceed 0.3%, while there are errors of up to 6% in voltage angles. Note that in most ofdistribution system analysis, only the voltage magnitudes are of importance and voltageangles are rather neglected.The MR-LPF and the LPF are also compared for the abnormal voltage conditions.In order to create abnormal voltages, the loads are scaled up with a factor λ to stress thesystem. The 873-node and 10476-node test systems are used here for the simulations [94].450.0000 0.0008 0.0016 0.0024 0.0031 0.0039 0.0047 0.00550204060η (%)Density(%)CZ = 0 CZ = 0.5 CZ = 1.5(a) The 70-node test system0.00 0.05 0.09 0.14 0.19 0.24 0.28 0.33020406080η (%)Density(%)CZ = 0 CZ = 0.5 CZ = 1.5(b) The 135-node test system0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35020406080η (%)Density(%)CZ = 0 CZ = 0.5 CZ = 1.5(c) The 3249-node test systemFigure 3.5: Histograms of relative errors for the LPF solution for various CZ .46−0.5 0 0.5 1 1.5 2 2.500.020.040.060.080.10.12CZAverageError(%)14-node 70-node 135-node 3249-nodeFigure 3.6: Average relative errors for the LPF solutions as a function of CZ .Table 3.4: CPU Times for 20 Runs of the LPF and the Implicit Z-Bus MethodTest System Method Iterations CPU Time (s)3249-node Implicit Z-Bus 4 1.142LPF 1 0.24470-node Implicit Z-Bus 5 0.0511LPF 1 0.0075The 873-node system consists of 7 feeders, 900 branches (27 normally-open switches),and a total load of 124.9 + j 74.4 MVA. The 10476-node system consists of 84 feeders,10736 branches (260 normally-open switches), and a total load of 1490.7 + j 886.7 MVA.Loads are modeled as constant-power elements. In the initial configuration and loading,the minimum nodal voltage for the 873-node and 10476-node systems are 0.90 and 0.96p.u., respectively. The Newton-Raphson power flow is used as a reference to measure theaccuracy of the proposed methods. In the initial configuration of the 873-node system, themaximum relative error in nodal voltages, indicated by µ, is 0.018% and 0.076% for theMR-LPF and LPF, respectively. To further stress the systems, all loads are scaled up by afactor λ and the corresponding values for µ are calculated. The results of these simulationsare shown in Figs. 3.8a and 3.8b. These figures show that µ monotonically increases forthe LPF as the system loading is increased. However, the MR-LPF shows different variationsin the value of µ as λ increases. This occurs due to switching between the regions shownin Fig. 3.4. It should be noted that the maximum error for the MR-LPF in the heaviestloading condition is 0.4% and 0.7% for the 873-node and 10476-node systems, respectively.4720 40 60 80 100 1200.940.960.981NodeVoltagemagnitude(p.u.)Implicit Z-bus LPF(a) Voltage magnitudes20 40 60 80 100 120−3−2−10NodeVoltageangle(deg.)Implicit Z-bus LPF(b) Voltage anglesFigure 3.7: Voltage magnitudes and angles for the 135-node test system obtainedby the LPF method and the implicit Z-bus method.3.3 The Linear Current Flow EquationsThe linear power flow (LPF) equations derived in Section 3.1 are based on the Nortonequivalent of the load synthesis, as shown in Fig. 3.3. The equivalent current source and481 1.2 1.4 1.6 1.8 2 2.2 2.4 2.610−310−210−1100λµ(%)0.650.70.750.80.850.9V minµMR-LPFµLPFVmin(a) The 873-Node Test System1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.510−310−210−1100101λµ(%)0.60.650.70.750.80.850.90.95V minµMR-LPFµLPFVmin(b) The 10476-Node Test SystemFigure 3.8: Maximum relative error in nodal voltages (µ) obtained using the MR-LPF and LPF versus loading factor (λ) in the 873-node and 10476-node sys-tems. The minimum nodal voltage (Vmin) is also shown on the right axis.impedance, in complex numbers, are calculated as:YLi = GL,i + j BL,i (3.16a)ILi = Ip,i + j Iq,i (3.16b)Using nodal analysis and the proposed load model in Chapter 2, a linear formulation wasobtained in Section 3.1 that calculates the nodal voltages. In order to find the branchcurrents, extra computations are required after finding the nodal voltages.In a different approach, let us consider the branch currents, instead of the nodal volt-49ages, as the independent variables. Also, replace the Norton equivalent of the loads shownin Fig. 3.3 by their The´venin equivalent, as shown in Fig. 3.9. Next, assume an arbitrarycurrent direction for each branch, as is done in Fig. 3.9. Considering the loop formedthrough the ground by the load connected to Node i, Branch k, and the load connectedto Node j, we can apply Kirchhoff’s Voltage Law to obtain:− Ei + ZLi(−Ik +∑l∈N iniIl −∑l∈NoutiIl)− ZkIk− ZLj(Ik +∑h∈N injIh −∑h∈NoutjIh) + Ej = 0 (3.17)In (3.17), N ini /Nouti are the sets of branches entering/leaving Node i excluding Branch k.Rearranging this equation, one can write for Branch k,(ZLi + ZLj + Zk)Ik − ZLi∑l∈N iniIl + ZLi∑l∈NoutiIl+ ZLj∑h∈N injIh − ZLj∑h∈NoutjIh) = Ej − Ei (3.18)The linear current flow (LCF) equations can now be written as:Z I = E (3.19)in which Z ∈ Cm×m is the network impedance matrix (m is the total number of branches);I ∈ Cm×1 is the vector of branch currents; E ∈ Cm×1 is the vector of equivalent voltagesource for each branch. The kth element of E is the difference between the The´veninZLi+− EiIk ZkZLj+−EjIl IhNode i Node jFigure 3.9: A generic part of distribution system: The´venin equivalent.50voltage sources located at both ends of Branch k, as in (3.18). The non-zero elements onthe kth row of Z, assuming Branch k is connected between Nodes i and j, can be definedas:zk,l =ZLi + ZLj + Zk l = k−ZLi l ∈ N iniZLi l ∈ NoutiZLj l ∈ N inj−ZLj l ∈ Noutj(3.20)The definition in (3.20) is readily implemented in a computerized program to form Z fora generic distribution system. A noticeable feature of the Z matrix is that, unlike theadmittance matrix, it is not necessarily symmetrical.Many optimization software packages only allow for operations on real numbers. There-fore, it is worthwhile to derive the current flow equations (3.19) in rectangular (Cartesian)coordinates. Assume Z = R + jX, I = Ire + jI im, and E = Ere + Eim. It is trivial toverify the following variant of (3.19):[R −XX R][IreI im]=[EreEim](3.21)which has 2m variables and 2m equations. Using the LCF equations of (3.19) or (3.21),the branch currents can be directly calculated, independently from the nodal voltages.The accuracy of the LCF formulation is from the same nature of the LPF formulation.Therefore, no simulation results are reported for the LCF formulation here. However, theapplication of the LCF equations in the formulation of network reconfiguration problem isdiscussed in Chapter 4.3.4 Sensitivity Factors for Distribution SystemsIn the preceding sections, the LPF formulation was established. Based on this LPF for-mulation, Current Transfer Distribution Factors (CTDFs) and Branch Outage DistributionFactors (BODFs) are derived here for distribution systems. The proposed BODFs are able tofind a power flow solution in cases of a disconnected network, for which the normal powerflow algorithms fail to find a solution. In other words, opening any branch in a radialnetwork leads to two disconnected networks. In such instances, the system admittancematrix is singular and one needs to recalculate the admittance matrix to account for thedisconnected nodes and branches. Using the proposed BODFs, this problem will not beencountered.513.4.1 Current Transfer Distribution FactorsUsing nodal analysis, the LPF equations can be written, in complex numbers format, asY¯ V = I (3.22)in which Y¯ is the modified admittance matrix with load admittances added to its corre-sponding diagonal elements. Suppose the voltages are obtained for an initial case using(3.22). Suppose also that a change occurs in the current injection at Node k. The aimhere is to find the changes in the branch flows without having to solve (3.22) again. Since(3.22) is linear, one can write∆V = Z¯∆I (3.23)where Z¯ is the inverse of Y¯ . Now, suppose it is desired to find the changes in branchcurrents. The current flowing through the branch connecting Node i to Node j is calculatedasfij = (vi − vj)yij (3.24)Therefore, the changes in fij due to the changes in the current injection at Node k can becalculated as∆fij = (∆vi −∆vj)yij (3.25)Substituting ∆vi and ∆vj from (3.23) into (3.25), one has∆fij = dij,k∆Ik (3.26)where dij,k = (z¯ik − z¯jk)yij is the CTDF. For large distribution systems with thousands ofnodes, calculating and storing Z¯ may be computationally inefficient when the CTDFs foronly a few branches are needed. Under such circumstances, since only a few elements ofZ¯ are required for the analysis, solving the following equation suffices:Y¯ zi = [0 . . . 0 1 0 . . . 0]t (3.27)where the ith element of the right-hand side vector is 1 and the rest are zeros. The solutionof this equation, i.e. zi, is the ith column of Z¯. For instance, to obtain Cij,k, (3.27) has tobe solved twice to obtain all the elements on the ith and jth columns of Z¯. Note that oneneeds to calculate Z¯ once and the result can be used later in all the sensitivity analysis.The CTDFs are calculated off-line and stored in a matrix D ∈ Cm×n, where m is thenumber of branches and n is the number of nodes. This matrix contains the sensitivities ofthe current flowing through each branch to a current injection at each node. The conceptcan be easily extended to multiple injections. For instance, if two current injections at52two different nodes are considered simultaneously, say Nodes k and l, then the changes inthe voltage at Node i is calculated as:∆vi = z¯ik ∆Ik + z¯il ∆Il (3.28)substituting (3.28) into (3.25), the changes in the current of Branch i-j is:∆fij =((z¯ik − z¯jk)∆Ik + (z¯il − z¯jl)∆Il)yij (3.29)This analysis shows that superposition can be applied to the CTDFs. In other words:∆fij = dij,k∆Ik + dij,l∆Il (3.30)3.4.2 Branch Outage Distribution FactorsAssume a branch is disconnected from a network. In order to find the new power flowsolution, (3.22) has to be solved with a modified Y¯ . In a radial distribution system,disconnection of a particular branch leads to the disconnection of all the subsequent nodesand branches that are connected to the substation through that particular branch. Inorder to find the power flow solution under such circumstances, a careful modification isrequired to form the new admittance matrix for the rest of the network still connected tothe substation. In a weakly-meshed network, a branch disconnection may not always leadto two isolated networks. Conventional power flow solution methods fail to solve a networkwith isolated sections and a reformation of the set of equations is required. For both radialand weakly-meshed networks, in order to avoid rebuilding Y¯ , sensitivity factors can beused to calculate the new solution with respect to the initial solution.In order to model the outage of Branch i-j, an appropriate fictitious current (Ix) canbe injected into Node i and be drawn from Node j, such that the resulting current throughBranch i-j (f ′ij) is only circulating between the two added current sources, i.e. Ix = f ′ij .Under these conditions, this branch is considered as “disconnected” from the rest of thenetwork, which can be verified by applying the Kirchhoff’s Current Law (KCL) to Nodes iand j. This concept is illustrated in Fig. 3.10. It is also important to note that by addingthe new current sources to Nodes i and j, the flows in the rest of the network are alsoaffected, which is equivalent to the disconnection of Branch i-j. The new flow in Branchi-j due to the added current sources can be calculated asf ′ij = (v′i − v′j)yij = fij + ∆fij (3.31)53Ix0f ′ij = Ix yij Ix0Node i Node jFigure 3.10: Modeling of a branch outage using fictitious nodal current injections.in which the variables with a prime stand for the new quantities in the modified network.Substituting the value of ∆fij from (3.30) resulting from the two current injections shownin Fig. 3.10, the new flow becomesf ′ij = fij + (dij,i − dij,j)Ix (3.32)Remember that the aim is to have f ′ij = Ix. Therefore, the appropriate fictitious currentinjections at Nodes i and j areIx =11− (dij,i − dij,j)fij (3.33)By injecting the calculated Ix at both ends of Branch i-j, this branch can be considered asisolated from the rest of the network. Finally, the change in the current flowing throughthe other branches l-r is calculated using (3.30) as∆flr = hij,lrfij (3.34)where hij,lr is the BODF given ashij,lr =dlr,i − dlr,j1− (dij,i − dij,j)(3.35)The new flows in the remaining branches after the outage of a particular branch arecalculated by adding the changes obtained in (3.34) to the initial branch currents. For anetwork with m branches, all the hij,lr form a m×m matrix, H ∈ Cm×m, with −1 in thediagonal elements.3.4.3 Active/Reactive Power InjectionIn distribution systems, determining the location for installing DG or capacitor bankshas been of interest for both academia and industry. Objectives such as minimizingthe network losses or improving its voltage profile are usually sought when installingDGs/capacitors. CTDFs derived in the previous section are not directly applicable when54active/reactive power injections are to be considered instead of current injections. Nev-ertheless, it is possible to derive equivalent current injections that closely mimic the ac-tive/reactive power injections. When a current Ik is injected to Node k, the resultingpower injection is calculated as:Sk = Vk I∗k = (V rek Irek + V imk I imk ) + j(V imk Irek − V rek I imk ) (3.36)in which the superscripts indicate real and imaginary. The equivalent active and reactivepower injections are:Pk = V rek Irek + V imk I imk (3.37a)Qk = V imk Irek − V rek I imk (3.37b)Assume now that an active power Pk is injected at Node k with zero reactive power Qk = 0.The equivalent current injections can be found from (3.37) as:Irek =V rek|Vk|2Pk, I imk =V imk|Vk|2Pk (3.38)The voltages in (3.38) are the final voltages after injecting Ik. By applying Ik the valuesof V rek and V imk will change from their initial values. Since the final values are not knownbeforehand, an estimation is used here. Assuming V rek ≈ 1 and V imk ≈ 0 gives a goodapproximation for the final values of the voltages at the nodes where a power injection isintroduced. The accuracy of these assumptions is illustrated in the next section.A similar analysis can be done for reactive power injection Qk with zero active powerPk = 0. From (3.37), the equivalent current injections can be calculated as:Irek =V imk|Vk|2Qk, I imk = −V rek|Vk|2Qk (3.39)The previous assumptions for the final values of voltages, i.e. V rek ≈ 1 and V imk ≈ 0, applyequally here.3.4.4 Performance IndicesVoltage Profile IndexOne important study in DG/capacitor placement is to find out which node is a bettercandidate to install a new source to improve the voltage profile. A good measure toevaluate the voltage profile of a network is to find the sum of deviations from the oneper-unit value. A Voltage Profile Index (VPI) can be defined based on this measure, as55follows:VPI =√√√√ 1nn∑i=1(1− |V ′i |)2 (3.40)where |V ′i | is the voltage magnitude of Node i (in per-unit) after injecting a current ata particular node in the system. Conventionally, for a current injection at each node inthe system, a new power flow solution is required to find the VPI. Using the proposedsensitivity analysis, however, only one power flow solution for the initial system suffices.In other words, the VPI due to a current injection at Node k, Irek + j I imk , can be calculatedasVPIk =√√√√ 1nn∑i=1(1− |Vi + z¯i,k (Irek + j I imk )|)2 (3.41)where Vi is the initial voltage at Node i.A similar analysis can be conducted to find the VPI when a branch outage occurs usingthe proposed BODFs. This is very useful when trying to find the best radial network froman initially meshed network that will result in the best voltage profile.Active Power LossesThe total active power losses in a distribution system when a new current injection isapplied can be calculated as:P ′L =∑i,jRij |f ′ij |2 (3.42)in which Rij and |f ′ij | are the resistance and the new current magnitude of Branch i-j,respectively. When the changes in branch currents are due to a change in power injectionat Node k, the proposed CTDFs can be readily used to find the new branch flows withoutrequiring a new power flow solution. Using (3.26), the power losses due to the new injectionat Node k, Irek + j I imk , can be calculated as:P ′L,k =∑i,jRij |fij + dij,k (Irek + j I imk )|2 (3.43)where fij is the initial flow in Branch i-j.A similar analysis can be conducted to find the losses when a branch outage occurs,using the proposed BODFs. This is very useful when attempting to find the best radialnetwork from an initially meshed network that generates the minimum losses.563.4.5 Simulation Results-Sensitivity FactorsMultiple studies are done in this section to show different applications of the proposedsensitivity factors. In the first study, the problem of DG/capacitor placement in a distri-bution system for loss reduction/voltage profile improvement is considered. In the secondstudy, the problem of network reconfiguration for loss reduction/voltage profile improve-ment is studied. For the purpose of illustration, the 33-node system of [95] is used here.The power flow data of this system is given in [95]. The topology of this system is re-produced here, for reference, in Fig. 3.11. In order to bring the simulation results closeto the constant-power load model, for the sake of comparison with traditional methods,and without loss of generality, the parameters defining the load voltage dependence arechosen as CZ = C ′Z = −1 and CI = C ′I = 2, as discussed in Section 3.1.5. Using theseparameters, the voltages obtained by the LPF and conventional Newton-Raphson powerflow solution, for the initial network, have a maximum error of 0.08%.DG/Capacitor PlacementPlacing a DG/capacitor in a radial network has two aspects that need to be determined:location and size. To optimally locate a DG/capacitor is usually more sophisticated thanits sizing since, in a general formulation, it involves discrete variables representing itslocation. In this section, the formulations proposed for VPI and power losses in Section3.4.4 are used to find the best location(s) for installing a new DG/capacitor.The 33-node system shown in Fig. 3.11 has a total load of 3.715 + j2.3 MVA withinitial losses of 202.7 kW and an initial VPI of 0.0596. Assume now that a 1 MW DG is tobe installed at one of the nodes in the network. The appropriate current injection at eachnode is calculated using (3.38). The VPI for each node is then calculated using (3.41) andis shown in Fig. 3.12a. The best nodes to place a DG in order to obtain the minimumdeviations of all the nodal voltages from 1 p.u., i.e., the minimum VPI, are Nodes 17 and18. The results obtained by the proposed method and by performing 32 Newton-Raphsonpower flow solutions give the same conclusion.It is also important to locate the best node for installing, for example, a 1 MVARcapacitor bank. The VPI obtained using the proposed method and by performing 32Newton-Raphson power flow solutions are given in Fig. 3.12b. The two candidate nodesare again Nodes 17 and 18. It is important to notice here that adding 1 MW active poweror 1 MVAR reactive power source have very similar effects on the system voltage profile.The small differences between the VPI values obtained using the proposed method and theNewton-Raphson power flow solutions are essentially due to the fact that after injectingpower to a node, the voltage of that particular node will change, compared to the assumedvalues. This leads to a difference between the final power injection and the planned value.5712 3456789101112 131415 16 17 181920212223 24 2526 27 28 2930313233123456789101112131415 16 171819202122 23 2425 26 27 28293031323334353637Figure 3.11: Initial configuration of the 33-node system [95]. The dotted-lines indi-cate open switches.Nevertheless, since this is a comparative study, the relative sensitivity of VPI to powerinjection at each node can still be observed from the results.In some cases, it is intended to reduce the losses by locating the best node for installinga DG/capacitor. In this example, the power losses are calculated using (3.43) and byconducting 32 Newton-Raphson power flow solutions. The results for 1 MW DG and1 MVAR capacitor bank placement at each node are shown in Figs. 3.13a and 3.13b,respectively. The best candidate node for installing either a DG or a capacitor is Node29, which reduces the losses to 127 kW in the case of a DG and 146 kW in the case of582 6 10 14 18 22 26 30 330246 Initial ValueNodeVPI(×10−2 )LPF NR-PF(a) 1 MW DG2 6 10 14 18 22 26 30 330246 Initial ValueNodeVPI(×10−2 )LPF NR-PF(b) 1 MVAR capacitorFigure 3.12: Voltage profile indices obtained using (3.41) (denoted as LPF) andNewton-Raphson power flow solutions (denoted as NR-PF) when a 1MW/MVAR DG/capacitor is installed at each node.a capacitor bank. It should be noted that installing extra active/reactive power supplyat Nodes 20, 21, and 22 increases the losses. This increase in the losses occurs since thegenerated power at the mentioned nodes needs to travel through a longer electric path toreach the other loads, causing higher power losses.Network ReconfigurationOne common approach for optimal network reconfiguration is to first close all the tieswitches, forming a meshed network. The switches are then opened one at a time based592 6 10 14 18 22 26 30 3350100150200250Initial ValueNodeP L(kW)LPF NR-PF(a) 1 MW DG2 6 10 14 18 22 26 30 3350100150200250Initial ValueNodeP L(kW)LPF NR-PF(b) 1 MVAR capacitorFigure 3.13: Total power losses obtained using (3.43) (denoted as LPF) and Newton-Raphson power flow solutions (denoted as NR-PF) when a 1 MW/MVARDG/capacitor is installed at each node.on the improvement they make in the objective function. Assume the 33-node test systemis to be reconfigured to minimize the losses. All the switches are first closed and the BODFsare calculated for this meshed system. The active power losses upon disconnection of eachbranch can then be calculated using (3.42), while the new flows are calculated using (3.34).In the first iteration, the changes in the total power losses due to the outage of eachbranch, i.e. P newL − P initialL , is calculated and reported in Table 3.5, second column. Notethat disconnection of Branch 1 is equivalent to disconnecting all the nodes in the network60and, therefore, is not feasible. Disconnection of Branches 9 and 10 reduces the losses in themeshed network. Therefore, Branch 9 is selected as the first candidate to be opened. Byremoving Branch 9 from the network, a new set of BODFs are calculated and the process isrepeated to find the next candidate. It is important to note that in the second iteration,disconnecting Branches 10 and 11 reduces the losses. However, checking this with Fig.3.11 reveals that disconnecting Branches 10 or 11 leads to an isolating part of network.The next candidate branch which its disconnection has the smallest effect in increasing thelosses is Branch 14. In the second iteration, opening Branches 10-13, 21, and 35 reduce thelosses while isolating part of the network. The next accepted candidate branch is Branch32 in this iteration. In the fourth iteration, opening any of Branches 8, 10-13, 15-17, 21,29-31, and 34-36 leads to an isolated network. The next candidate is, therefore, Branch7. In the fifth and final iteration, Branch 37 is the best candidate.According to the above analysis, the final configuration has the following open switches:9, 14, 32, 7, and 37. This is the global optimal solution for this network, which reducesthe total losses from 202.7 kW down to 139.6 kW, i.e. 31.1% reduction. If the BODFs arenot used, one needs to run 36 × 35 × 34 × 33 × 32 = 45, 239, 040 power flow routines, ifno smart screening method is used to discard disconnections that lead to islanding. Thisnumber of solutions would be needed to evaluate the losses upon disconnection of eachbranch at each iteration. Recalculation of the BODFs at each iteration for the new networkrequires only one LPF solution, which in total requires 5 solutions. It should be noted thatsince the connectivity of the network is checked at every iteration, the final configurationis guaranteed to be radial.With the proposed method, it is easy to detect a disconnected network just based onthe values of the losses at each iteration. Upon opening the first few switches, the lossesmay decrease. In such cases, the switch that generates the minimum losses is chosen.After a few switches are opened, there will be an iteration at which a disconnected graphcorresponds to the minimum-loss configuration (e.g., the second iteration in the 33-nodecase). From this iteration on, all the solutions that generate less losses correspond to adisconnected graph. In other words, any further valid switch-opening action would ratherincrease the losses. Therefore, the switch-opening action that produces the minimumpositive difference between the previous losses and the new losses is the next candidateswitch that forms a connected graph. For instance, after the second iteration in Table 3.5,the minimum positive value at each column is the candidate, shown in bold.61Table 3.5: Changes in the Total Power Losses (kW) Due to a Branch OutageBranch Iter.#1 Iter.#2 Iter.#3 Iter.#4 Iter.#52 445.5 458.0 536.5 560.9 -84.73 51.7 52.2 54.1 54.6 75.44 40.9 41.3 43.0 43.3 58.35 36.4 36.7 38.2 38.5 51.26 4.5 4.6 4.7 6.5 -10.07 0.5 0.5 0.6 0.8 -8 5.5 7.4 14.9 -24.1 -24.59 -0.04 - - - -10 -0.03 -3.6 -3.3 -3.3 -3.611 0.2 -6.7 -6.1 -6.0 -6.512 3.5 5.2 -11.0 -11.0 -11.613 2.1 3.0 -7.9 -7.8 -8.314 0.3 0.2 - - -15 8.1 8.1 7.9 -18.7 -19.616 5.5 5.5 5.2 -15.5 -16.417 3.4 3.4 3.2 -11.9 -12.818 72.7 75.1 80.3 82.6 -44.519 62.1 64.2 68.4 70.3 -44.220 52.5 54.3 57.5 59.1 -43.721 13.4 20.8 -20.9 -20.8 -21.722 102.3 102.3 107.7 121.2 123.123 88.8 88.8 93.6 105.3 106.524 40.8 40.7 43.3 49.0 48.525 17.2 17.4 18.1 23.4 26.626 14.2 14.3 14.9 19.5 22.327 11.5 11.6 12.1 16.0 18.428 9.2 9.3 9.7 13.1 15.129 59.5 59.7 73.5 -54.6 -51.530 9.6 9.6 12.6 -28.2 -26.831 2.9 2.9 4.2 -17.6 -16.832 0.3 0.3 0.3 - -33 6.8 7.4 11.4 11.6 -31.134 3.3 5.2 10.9 -21.4 -22.135 9.1 14.2 -17.9 -17.8 -18.736 1.2 1.2 1.0 -5.4 -5.937 12.0 12.0 13.0 15.2 14.662Chapter 4Optimal Network Reconfiguration4.1 Network Reconfiguration Using Mixed-IntegerProgramming4.1.1 LPF-based Mixed-Integer ProgrammingBased on the load model described in Chapter 2, the LPF formulation was derived inChapter 3. The matrix formulation in (3.14) is reproduced here in expanded form, asfollows: n∑k=1(Gm,kV rek −Bm,kV imk)= Ip,m (4.1)n∑k=1(Gm,kV imk +Bm,kV rek)= Iq,m (4.2)where n is the number of nodes; m and k are the nodes at the two ends of the branchconnecting Node m to Node k; V re and V im are the real and imaginary parts of voltages.The most important part of the network reconfiguration formulation is how to modelthe status of the switches, i.e. “on”/“off”. The first method is to multiply the binaryvariables by the branch flow limits, as done in [60]. The second method is to multiplythe binary variables by branch admittances, which subsequently affects the admittancematrix. In the first method, if a branch is offline, then its flow is forced to zero while inthe second method, its impedance is forced to vanish. The second method is followed heresince the first method is not applicable in a nodal analysis formulation. The followingequation describes the modified admittance matrix:63Ym,k =um,k ym,k m 6= k−n∑k=1k 6=mum,k ym,k + jQC,m − YL,m m = k (4.3)in which Y (= G+ jB) is the modified admittance matrix; ym,k is the admittance of thebranch connecting Node m to Node k; QC and YL are the shunt capacitance and loadadmittance at each node, respectively; um,k stands for the status of Branch m-k. Notethat if there is no branch between Nodes m and k, the corresponding um,k is zero.The current flowing through Branch m-k (Im,k) is calculated as:I2m,k = um,k[G2m,k +B2m,k][(V rem − V rek )2 + (V imm − V imk )2] (4.4)The total network losses are given by the following formulas:Ploss =∑m,km1800[101] 280.2 611.3 0.959 0.977 4473MIQCP* 258.3 566.2 0.964 0.979 188.4880-node**Initial* 2858.8 2690.2 0.913 0.974 -[60] 999.1 1223.9 0.982 0.990 3192MIQCP* 958.7 1172.0 0.982 0.990 1134* For MIQCP and initial case, loads are modeled as voltage-dependent with CZ = C ′Z = 0.5. In other references, the loads aremodeled as constant-power.** Relative optimality gap is 1% for this case.tion. As can be seen, usually the first few switching actions produce the largest reductionin the total losses. For example, by halving the limit on S, only about 0.9% increase in thelosses was observed for the 136-node system. The CPU time was also reduced comparedto the previous case, while the number of integer variables is essentially the same. Itshould be noted that the solution speed for a Mixed-Integer Programming (MIP) problemis mainly governed by the number of integer variables. In a real system, the number ofswitches is small; hence, the computation time may not be a challenge for the applicationof the proposed framework for large-scale systems.In general, especially with mechanical switches, it is desired to have the best configura-71Table 4.4: Optimal Configurations for Test SystemsSystem Open Lines14-node 6-8,7-9,5-1433-node 7-8,9-10,14-15,32-33,25-2970-node 14-15,9-38,15-67,49-50,39-59,38-43,9-15,21-27,28-29,62-65,40-4484-node 7-6,13-12,18-14,26-16,32-28,34-33,39-38,43-11,72-71,83-82,55-5,41-42,63-62119-node 23-24,26-27,35-36,41-42,44-45,51-65,53-54,61-62,74-75,77-78,86-113,95-100,101-102,89-110,114-115136-node6-7,10-32,57-61,78-125,20-130,137-138,59-145,139-154,141-154,155-156,154-204,211-212,138-217,125-219,141-220,222-223,144-145,43-46,63-64,130-131,214-215880-node89-90,136-137,146-147,164-165,195-196,287-288,293-294,290-312,318-336,414-415,417-418,457-458,499-500,601-602,621-622,634-635,636-637,642-643,653-704,736-824,848-849,345-500,463-824,460-440,580-800, 850-250,10-470tion with the minimum number of required switching actions. In order to compromise be-tween these two objectives, a multi-objective optimization problem may be defined whichminimizes both the losses and the number of switching actions simultaneously. This re-quires an estimation of the monetary cost of each switching operation and of the losses.Since this data is not available to the author, in the examples an arbitrary limit on thenumber of switching actions is used.4.1.2 LCF-based Mixed-Integer ProgrammingThe LCF equations derived in Chapter 3.3 are used here to formulate the network recon-figuration problem. From a mathematical point of view, the LCF formulation suits thereconfiguration problem better than the LPF variant. This is mainly to avoid the bilinearterms resulting from multiplication of binary and continuous variables. With LCF, thedisconnection of each line is simply modeled by forcing its current to zero, while with LPFthe corresponding line impedance has to be taken out from the impedance matrix.The objective of the problem can be load balancing among feeders, loss reduction,voltage profile improvement, etc. Here, the loss minimization problem is studied. Thetotal losses in the network can be calculated as:Ploss =m∑k=1Rk(|Irek |2 + |I imk |2) (4.20)72Table 4.5: Results Obtained By Assuming A Limit On theNumber of Switching Actions*14-node Smax 4 2 0Ploss (kW) 553.7 570.7 593.933-node Smax 8 4 0Ploss (kW) 122.8 126.3 167.370-node Smax 10 6 0Ploss (kW) 172.3 173.9 189.284-node Smax 16 8 4 0Ploss (kW) 424.2 426.5 437.7 470.2119-node Smax 24 10 4 0Ploss (kW) 765.9 798.1 853.8 1029.8136-node Smax 24 12 4 0Ploss (kW) 258.2 260.8 263.4 287.4880-node Smax 48 12 4 0Ploss (kW) 958.7 1098.3 1470.0 2858.8* These results are obtained assuming: CZ = C ′Z = 0.5.The constraints to the minimum-loss network reconfiguration are established next.The Current Flow EquationsThe LCF equations can be written in expanded form asm∑k=1(RkIrek −XkI imk)= Erek (4.21a)m∑k=1(XkIrek +RkI imk)= Eimk (4.21b)Note that (4.21a) and (4.21b) add 2 ×m linear equality constraints to the optimizationproblem. This system of equations has a single unique solution. In other words, the matrixrepresenting the system of linear equations has full rank. When a few of the branches aredisconnected to retrieve a radial structure, the corresponding equations to those brancheshave to be eliminated. Otherwise, the system of linear equations is overdetermined and,thus, infeasible. The process of branches disconnection is carried out as follows.Let us denote the status of branch k by uk. It is then desired to activate the kthequation in (4.21a) and (4.21b) if k = 1, and deactivate it otherwise. This particular typeof variables, i.e. uk, are called indicator variables [102]. An indicator variable, when iszero, forces some of the problem variables to assume a fixed value, and, otherwise, forcesthem to belong to a convex set. Commercial mixed-integer programming solvers such as73CPLEX are capable of efficiently handling this feature of the problem [103].An indicator variable can also be represented as a disjunction [104]. Disjunctive pro-gramming deals with logic-based constraints, i.e. constraints that are only active whenan indicator variable is true. One reformulation of a disjunction, e.g., (4.21), takes thefollowing form:(uk − 1)M ≤m∑k=1(RkIrek −XkI imk )− Erek ≤ (1− uk)M (4.22a)(uk − 1)M ≤m∑k=1(XkIrek +RkI imk )− Eimk ≤ (1− uk)M (4.22b)in which M is a large-enough positive scalar. This reformulation is also known as “Big-M” formulation. It is tricky to find a good value for M since a large value leads to loosebounds and a small value may lead to an infeasible problem. Commercial solvers such asCPLEX are capable of finding good values for M during the solution process to tightenthe bounds. This useful feature of the solver is utilized here to enforce the disjunctions.Branch AmpacitiesThe magnitude of the current flowing through a branch is bounded by the thermal limitof the conductors. This constraint can be represented as:|Irek |2 + |I imk |2 ≤ |Imaxk |2 uk (4.23)where Imax is the maximum current magnitude allowed in a branch. Note that uk = 0forces both Irek and I imk to be zero.Nodal Voltage LimitsThe voltage at each node has to be within certain operational limits, e.g., ±5%. In orderto derive the voltages from the branch currents, the The´venin equivalent of loads shownin Fig. 3.9 is used. The voltage at Node i can be calculated asV rei + jV imi = −Ei + ZLi(−Ik +∑l∈N iniIl −∑l∈NoutiIl) (4.24)This equation is in the form of complex numbers. The real and imaginary parts of nodalvoltages, which are linear combinations of currents, can be used to impose the voltagemagnitude limits:V mini ≤√|V rei |2 + |V imi |2 ≤ V maxi (4.25)74Recall the assumption of small voltage angles in DS. This assumption allows for eliminatingthe term |V imm |2 in (4.25), giving the following box constraint:V mini ≤ V rei ≤ V maxi (4.26)where V rei can be retrieved from (4.24) as a function of currents.Radiality ConstraintWhen radiality of the final configuration is imposed, the formulation proposed in [59] isused. This formulation is given in (4.15). The only difference is in the notation of um,k,which stands for the status of the branch connecting Nodes m and k. In this section, thebranch index is used instead and, therefore, um,k is replaced by uk, representing the statusof Branch k.4.1.3 A MILP Formulation of the Network Reconfiguration ProblemThere are two quadratic terms in the MIQCP formulation, i.e. (4.20) and (4.23). In thissection, effective relaxation techniques are applied to the quadratic terms to find a MILPformulation for the problem.Active Power LossesThe total active power losses are calculated in (4.20). Each element of the summation, i.e.Rk(|Irek |2 + |I imk |2), is a convex quadratic term. A Piecewise Linearization (PWL) techniqueis proposed in [105] to linearize these type of functions. This technique has been provenefficient in many applications, e.g., in linearizing the generators cost functions [106]. Thedetails of this linearization process are provided in Appendix A. When applying the PWLtechnique, each term in the objective function takes the form of a maximum over a set ofaffine functions. Define λk asλk = max1≤i≤s{aiIrek + biI imk + ci} (4.27)which is the PWL approximation of |Irek |2 + |I imk |2. The objective function then takes thefollowing form:minm∑k=1Rk max1≤i≤s{aiIrek + biI imk + ci} (4.28)which is a min-max optimization problem. This can be easily reformulated as a linearprogramming problem, as discussed in [107]. The following problem is equivalent to (4.28):75minm∑k=1Rkλk (4.29a)s.t. λk ≥ aiIrek + biI imk + ci, ∀i = 1, . . . , s. (4.29b)The number of linear pieces, i.e. s, can be chosen in a way to satisfy a particular accuracyin the approximation.Branch AmpacitiesThe current ampacity constraint in (4.23), for a fixed uk = 1, is a quadratic constraint ofthe form (B.1). The branch status uk is enforced later by adding (4.30) to the problem.In Appendix B, a relaxation method for these constraints is proposed. Using a hexagonrelaxation (i.e. n = 6), the linear constraints of (B.6) replace (4.23). In addition, thefollowing two constraints are required enforce the branch status:− ukImaxk ≤ Irek ≤ ukImaxk (4.30a)− ukImaxk ≤ I imk ≤ ukImaxk (4.30b)Theses equations would be inactive if uk = 1, since they represent a larger feasible regionthat contains the hexagon defined by (B.6).Simulation ResultsIn this section, the proposed LCF-based MIQCP and MILP formulations are applied toseveral distribution test systems. The parameters for load models are assumed to beidentical as CZ = CI = 0.5. The value of losses and optimum configurations obtained usingthe MIQCP and MILP methods are the same as the ones obtained by the LPF-based MIQCPmethod reported in Tables 4.3 and 4.4. The linearization technique used to represent theobjective function in the MILP formulation, i.e. (4.29), introduces a negligible error in thevalue of losses. It should be noted that the direction at which the losses increase/decreaseplays a more important role in the optimization problem than its exact value. Due to thisreason, both the MILP and MIQCP methods find the same optimal configurations.Optimal reconfiguration for loss minimization invariably improves the voltage profileof the network. This fact has been observed by the author in all the test systems, andthe observations in [44] also support this fact. Therefore, it may be possible to drop thevoltage limit constraints, i.e. (4.26), in most of the cases when minimizing the losses. Theproposed formulation provides such a structure that allows for completely dropping thevoltage variables and constraints, without affecting the results. The problem dimensionwould then reduce substantially. This fact further supports the suitability of the proposed76Table 4.6: Computation Times (s) for the LCF-based FormulationsTest System (# Nodes)33 70 119 136 415 880**LCF-based MIQCP 35.2 485 1009 1785 1800+ 1800+LCF-based MILP 0.15 0.91 2.8 5.6 116 398LPF-based MIQCP 3.2 5.7 39.4 188.4 953 1134[60] 12.8 11310 N/A 2.35* N/A 3192*[59] N/A N/A N/A 1800+ 1800+ 1800+[101] 19 N/A 4007 4473 14256 N/A* Suboptimal solutions obtained by intense relaxations [60].** Relative optimality gap is considered 1% for this case.LCF-based formulation for network reconfiguration.The computation times for the test systems are given in Table 4.6. All the simulationswere done in the GAMS platform on an Intel Core i7-2600, 3.4 GHz CPU, 8 GB RAM,and 64-bit operating system machine.It is important to notice the difference between the computation time for MILP andMIQCP formulations. Although the size of the MILP problem is larger than the MIQCPversion, the computational complexity of the linear programming problem solved at eachiteration is less than the equivalent conic quadratic programming problem. The mainchallenge in the MIQCP problem is to handle the conic constraints. In other words, if theconic constraints are relaxed, then all the constraints are linear and only the objective isquadratic. Linearizing the quadratic objective only makes a small difference, as comparedto relaxing the conic constraints, in the solution time. When the control variables are onlydiscrete variables, which is the case in network reconfiguration, using linearized objectivefunction is well-justifiable.4.2 Mathematical Representation of Radiality Constraint4.2.1 Planar GraphA graph is called planar if, with any rearrangement of its vertices, it can be drawn ona plane without having intersecting edges. A planar graph with n vertices and m edgesdivides the plane into f faces. Euler’s formula [108] suggests the following relation for aplanar graph:f = n−m+ 2 (4.31)77For example, consider the graph shown in Fig. 4.3. The faces are the regions on theplane separated by the graph’s edges. The outside infinite face (shown by ”E”) is alsocounted. In Fig. 4.3, the faces are identified by capital letters. There are two necessary,but not sufficient, conditions for a graph to be planar [108]:n ≥ 32f , n ≤ 3m− 6 (4.32)Apart from those necessary conditions, there is a theorem in [108] that provides neces-sary and sufficient conditions for planarity. Before referring to that theorem, two particulargraphs, known as Kuratowski’s graphs, need to be introduced. The graphs shown in Fig.4.4, known as K5 and K3,3, are Kuratowski’s two graphs. Another preliminary concept isthat of homeomorphic graphs. Two graphs are homeomorphic if one can be obtained fromthe other by adding new edges in series to the existing ones or by merging already-existingedges that are in series. As per [108], a necessary and sufficient condition for a graphto be planar is that it does not contain either of Kuratowski’s two graphs, or any graphhomeomorphic to either of them.According to the author’s experience, all distribution systems encountered satisfy theconditions for planarity. Overhead lines are mainly built along land corridors, and be-cause they are geographically distributed in a plane (which is the definition of a planargraph), the natural intuition is that a DS has a planar graph representation. There arealso formulated algorithms to check whether an arbitrary graph is planar, e.g., [109].4.2.2 Dual GraphThe dual graph G∗ of a planar graph G is defined as follows [108]:• For each face of G, there is one corresponding vertex in G∗.12345678912345678912 4567893123456789DBECAABCD EFigure 4.3: The 9-node network. The letters show the faces of the graph.78(a) K5 (b) K3,3Figure 4.4: The Kuratowski’s two graphs.12345678912345678912 4567893123456789DBECAABCD EFigure 4.5: Dual graph (dotted lines) of the 9-node network.• For each edge joining two neighbouring faces in G, there is a corresponding edgebetween the two vertices in G∗.• For any pendent edge (an edge with only one vertex connected to it) in G, there isone self-loop at the corresponding vertex in G∗.From the above definition, it immediately follows that if G has n vertices, m edges and ffaces, then G∗ has f vertices, m edges and n faces [108]. Figure 4.5 depicts the dual graphof the 9-node graph shown in Fig. 4.3. As can be seen in Fig. 4.5, there may be morethan one edge between two vertices in the dual graph which have to be distinguished.4.2.3 The Spanning Tree ConstraintThe radiality constraint in a DS is identical to the spanning tree constraint in graphtheory. The minimum spanning tree in a weighted undirected graph is the subgraphthat is a tree and the sum of its weights is the minimum possible. This problem is well-79addressed in the literature [110]. Also, an integer linear programming formulation for thisproblem specifically designed for planar graphs is proposed in [110]. This method is brieflyexplained in the following.An undirected graph is first converted into a directed graph. Define the followingvariables and sets in G:• xij : status of the edge connecting vertex i to vertex j.• wij : weight of the edge connecting vertex i to vertex j.• Ni: set of vertices directly connected to vertex i.and in G∗:• y(k,l),e: status of the edge(s) connecting vertex k to vertex l. The index e is used todistinguish between multiple edges connecting the same two vertices.• Mk: set of vertices directly connected to vertex k.• Sk,l: set of multiple edges between vertices k and l.The minimum spanning tree problem is formulated as follows [110]:Minimize∑i,jwijxij (4.33)Subject to:∑j∈Nixij = 1, 1 ≤ i ≤ n− 1 (4.34a)∑l∈Mk∑e∈Sk,ly(k,l),e = 1, 1 ≤ k ≤ f − 1 (4.34b)xij + xji + y(k,l),e + y(l,k),e = 1, For all edges in G (4.34c)Note that since (4.34c) has one equation for each edge in G, it represents m constraints.Also, if i = n or j = n (k = f or l = f) in G(G∗), xij(y(k,l),e) only exists in one directionterminating in vertex n(f) in G(G∗), i.e. the last vertex.The set of constraints in (4.34) enforce the spanning tree constraint. These constraintsare used later in the network reconfiguration problem to impose radiality.4.2.4 Formulation of Radiality ConstraintsThe formation of the radiality constraint for the 9-node system shown in Fig. 4.5 isdiscussed here in detail as a reference.The needed sets are formed as follows:80N1 = {2, 6} N2 = {1, 3, 5} N3 = {2, 4}N4 = {3, 5, 8} N5 = {2, 4, 9} N6 = {1, 7, 9}N7 = {6, 8} N8 = {4, 7, 9} N9 = {5, 6, 8}MA = {B,C,D,E} MB = {A,D,E}MC = {A,D,E} MD = {A,B,C,E}ME = {A,B,C,D}SA,E = {1, 2} SB,E = {1, 2} SC,E = {1, 2}In the following, examples are given on how to write the problem constraints. In theprimal graph (4.34a) isFor i = 1, x1,2 + x1,6 = 1.For i = 2, x2,1 + x2,3 + x2,5 = 1.For i = 3, x3,2 + x3,4 = 1.For i = 4, x4,3 + x4,5 + x4,8 = 1.For i = 5, x5,2 + x5,4 + x5,9 = 1.For i = 6, x6,1 + x6,9 + x6,7 = 1.For i = 7, x7,6 + x7,8 = 1.For i = 8, x8,4 + x8,7 + x8,9 = 1.In the dual graph (4.34b) isFor k = A, yA,B + yA,C + yA,D + yA,E,1 + yA,E,2 = 1For k = B, yB,A + yB,D + yB,E,1 + yB,E,2 = 1For k = C, yC,A + yC,D + yC,E,1 + yC,E,2 = 1For k = D, yD,A + yD,B + yD,C + yD,E = 1For (4.34c), each branch has one equation:For Branch 1-2, x1,2 + x2,1 + yA,E,1 + yE,A,1 = 1For Branch 2-3, x2,3 + x3,2 + yB,E,1 + yE,B,1 = 1For Branch 3-4, x3,4 + x4,3 + yB,E,2 + yE,B,2 = 1For Branch 4-5, x4,5 + x5,4 + yB,D + yD,B = 1For Branch 1-6, x1,6 + x6,1 + yA,E,2 + yE,A,2 = 1For Branch 6-7, x6,7 + x7,6 + yC,E,1 + yE,C,1 = 1For Branch 7-8, x7,8 + x8,7 + yC,E,2 + yE,C,2 = 1For Branch 8-9, x8,9 + x9,8 + yC,D + yD,C = 1For Branch 2-5, x2,5 + x5,2 + yA,B + yB,A = 1For Branch 6-9, x6,9 + x9,6 + yA,C + yC,A = 181Table 4.7: Values for βij According to (4.15) for the Network in Fig. 4.6Value Variables1 β61 β52 β23 β34 β84 β45 β95 β760 β21 β96 β89 β87For Branch 5-9, x5,9 + x9,5 + yA,D + yD,A = 1For Branch 4-8, x4,8 + x8,4 + yD,E + yE,D = 1If there is a pendent node in the network, i.e. a node that has only one branch connectedto it, that branch is definitely in the spanning tree since its disconnection renders thegraph disconnected.4.2.5 Inadequacy of Existing Methods in Representing the RadialityConstraintThere are several different methods proposed in the literature for representing the radialityconstraint. These methods, however, may not be adequate, as is shown through exampleshere. The first method of representing the radiality constraint, which is the simplest, isto require the number of branches to be exactly equal to the number of nodes minus one.In other words,∑ij∈Wuij = n− 1 (4.35)where uij is the binary variable standing for branch i-j status (0: “open”, 1: “close”); Wis the set of all branches. This criterion has been used to enforce the radiality constraintin, e.g., [17], [111], [48]. However, this constraint does not guarantee the connectivity ofthe resulting network, as is also acknowledged in [101], [112]. As a counterexample, lookat the topology shown in Fig. 4.6 for the network shown in Fig. 4.3. It is trivial to checkthat the network in Fig. 4.6 satisfies (4.35), but is not connected.The second approach used in the literature for representing radiality is to model thenetwork as a directed graph, e.g., [60], [59]. The clear statement of the constraints isgiven in (4.15). It is trivial to check that the network shown in Fig. 4.6 satisfies all theconstraints in (4.15). The values for β constructing this network are given in Table 4.7.Only non-zero values of βij are reported in Table 4.7.The reason that the mentioned constraints representing radiality still work, whenembedded in a network reconfiguration problem, is discovered by the authors. The dis-connected network leads to an infeasible power flow solution. In other words, the connec-tivity constraint is imposed by the power flow equations. This fact is also emphasized in[101] and the network flows are used to impose the connectivity of the network. However,during the process of solving a mixed-integer programming problem such as the one in8212345678912345678912 4567893123456789DBECAFigure 4.6: The network obtained by applying the conventional radiality con-straints.[59], [60], and the proposed formulations in Section 4.1, infeasible configurations may begenerated, which is due to insufficiency of the constraints representing the radiality. Theinfeasible subproblems in a branch-and-cut algorithm lead to a larger number of iterationswhich, in turn, increases the CPU time of the whole solution process. The most severecase is when decomposition algorithms are used to solve mixed-integer problems, e.g., [58].When Bender’s Decomposition is used [58], the master problem, which mainly deals withthe integer variables, may generate infeasible configurations that would not be known tobe infeasible until the subproblem is solved. Moreover, in most of the heuristic methods,many infeasible configurations are generated first and then discarded by checking the radi-ality constraint. The process of determining and discarding the infeasible solutions slowsdown the whole solution process, leading to an unnecessarily large CPU time.Another method for representing the radiality constraint is to employ the branch-to-node incidence matrix, e.g., [52], [113]. The elements of the incidence matrix are all -1, 0,or 1. It is known that if a graph represents a spanning tree, then the determinant of theincidence matrix must be -1 or 1 [114]. In other words, if the determinant is zero, thenthe obtained subgraph is Not a spanning tree. One deficiency of this method is that itcannot be explicitly stated in a mathematical formulation that can be used in conventionaloptimization models. Another problem with this method is that it needs calculations withhigh computational complexity. Therefore, its direct application in mathematical modelsof network reconfiguration problem is impractical.4.2.6 Simulation ResultsIn this section, the problem of network reconfiguration for loss reduction, as formulatedin Section 4.1.1, is solved for several test systems using the new method for representing83Table 4.8: CPU Times When Different Meth-ods are Used To Represent Radiality Con-straintTest Case Losses(kW) T1(s)* T2(s)**Initial Optimum14-node 514 468.3 0.14 0.1633-node 202.7 139.6 3.2 3.070-node 227.5 201.4 5.7 4.284-node 532 469.9 9.4 7.8119-node 1298.1 869.7 39.4 30.1135-node 320.4 280.2 188.2 132.5* T1: CPU time when (4.15) is used.** T2: CPU time when (4.34) is used.the radiality constraint. The mixed-integer programming problem is written in AIMMSenvironment [99] and GUROBI is used to solve it.Simulation results for the test systems obtained using the radiality constraint in (4.15)as well as the results obtained using the radiality constraint in (4.34) are given in Table 4.8.The power losses are calculated in the post-process using a constant-power load model,after the configuration is determined by the voltage-dependent load model formulation.This is done for the purpose of comparison with reported values in other literature, e.g.,[59]. The optimality of the results is guaranteed by the solver. A comparison of theoptimal losses and configurations with previous literature is done in Section 4.1.1 and is notreproduced here. It is important to notice that the amount of saving in CPU time increasesas the size of the problem increases. For instance, 17%, 24%, and 30% reduction in CPUtime are achieved for the 84-node, 119-node, and 135-node systems, respectively. This isdue to tighter relaxations provided by the subproblems generated during the branch-and-cut algorithm. When a disconnected network is generated by (4.15), the correspondingquadratic programming relaxation would be infeasible since the power flow equations areimpossible to satisfy in a disconnected network. This fact does not affect the quality of thefinal solution. However, it causes the whole solution process to take longer to converge,as can be clearly seen in the reported CPU times. The optimal configurations are thesame as the ones reported in Table 4.4. There is no difference between the configurationsobtained when (4.34) is used and the configurations obtained when (4.15) is used.844.3 Network Reconfiguration Using Minimum SpanningTreeThe tendency towards optimizing the utilization of the current infrastructure in powersystems has been significantly increased in the past decade. The topological reconfigura-tion problem was introduced in Section 4.1 in this chapter. Although the deterministicapproaches developed for network reconfiguration, such as [59], [60], provide many advan-tages in terms of flexibility of the formulation (e.g., for adding extra constraints, aimingat different objectives, knowledge about the optimality gap at every iteration, etc.), theydo not take advantage of a good initial solution to start the MIP solver. Many commercialMIP solvers allow for a “warm- start”, i.e. starting from a known solution, to speed upthe search for the global optimum, e.g., CPLEX [103] and GUROBI [100].In this section, a fast method is sought which provides a suboptimal solution for theminimum-loss reconfiguration problem with a small optimality gap. This can replace thefirst stage of the MIP solution process that most of the commercial solvers have as a built-inroutine. This routine searches for a feasible solution using some general-purpose heuristicmethods. Those heuristics, however, are relatively slow and often find a solution with alarge optimality gap. Providing a good initial solution to start with, the solution time canbe significantly reduced [103].In this section, the problem of minimum-loss DS reconfiguration is mapped into a Min-imum Spanning Tree (MST) problem based on some assumptions. The main assumptionhere, based on engineering knowledge, is that the meshed network is a good solution (ifnot the best) for loss minimization when the radiality constraint is relaxed. The sameassumption has been made in [112] with a slight difference. In [112], it is assumed thatthe meshed network generates the least possible losses, which may not be always true aswill be shown later on in this chapter. However, the difference between the losses in themeshed network and the best-possible network configuration is negligible. The secondassumption in the present method is that reconfiguring the network in order to reducethe losses will implicitly improve the voltage profile. This assumption has been previouslyvalidated in [44], and is also validated here on various test systems. The voltage profilecan also be modified by adjusting capacitor banks, voltage regulators, and transformerstap positions, outside of the reconfiguration process.Based on the assumptions mentioned in the previous paragraph, the problem of DSreconfiguration for loss reduction is defined as “finding a spanning tree that imitates, asclosely as possible, the same flow pattern as the meshed network”. The idea of startingfrom a meshed network and opening switches sequentially has been proposed in [44], calledDISTOP. There are two main differences between DISTOP and the proposed methodhere. The first difference is that after every switching action, a new power flow solution85is required in DISTOP, while only one power flow solution is sufficient for the proposedalgorithm. The other difference is how the algorithms check whether a configuration isradial. It is stated in [44] that the switches are opened one after another until the networkbecomes radial. However, it is not explicitly explained how to check the radiality. In manycases, disconnecting a switch that carries the minimum current leads to a disconnectednetwork. It is not trivial to check if the network remains connected after opening a par-ticular switch. The proposed algorithm here, on the other hand, is guaranteed to providea radial configuration.The proposed heuristic method here has the advantage of using a fast and robust al-gorithm that also finds a high-quality solution (with small optimality gap). The efficientMST algorithms developed in the literature such as Kruskal [115] and Prim [116] algo-rithms can be employed to solve the formulated problem. A refinement to the solution isdone by sensitivity analysis around the neighborhood of the candidate tie switches.4.3.1 Background: Minimum Spanning TreeA spanning tree is a subgraph of an undirected graph that contains all the vertices (i.e. itis connected) and has no cycles (has no loops). The MST in a weighted undirected graph isthe subgraph that is a spanning tree, and the sum of its weights is the minimum possible.This problem is well-addressed in the literature under the name “minimum spanning tree”(MST) and there are efficient algorithms such as Kruskal [115] and Prim [116] to solve theproblem. By negating the weights, the algorithms for MST will find the maximum spanningtree [117]. The time complexity (the amount of time taken by the algorithm to run as afunction of the network size) of a version of Prim’s algorithm which is suitable for sparsegraphs is O(e log v), which is similar to the time complexity of Kruskal’s algorithm.Prim’s algorithm is used in this paper to find the MST. There is basically no advantageof using Prim’s algorithm over Kruskal’s, and one may choose the other. It is importantto emphasize that these algorithms provide the optimal solution to the MST problem. Thestructure of the Prim’s algorithm for a weighted undirected graph G(V,E) with V verticesand E edges is summarized in the following steps.1. Choose any vertex r in V to be the root node. Set Vt = {r} and Et = ∅.2. Find an edge with the smallest weight such that one of its end points is in S andthe other is in V \Vt. Add this edge to Et and its new vertex to S.3. If V \Vt = ∅, then terminate. Otherwise, go back to Step 2.Prim’s algorithm is applied to the 14-node graph shown in Fig. 4.7a. The final spanningtree obtained using the MST algorithm is shown in Fig. 4.7b. The details of each step aregiven in Table 4.986Table 4.9: MST Search For The 14-Node Graph Using Prim’s AlgorithmIteration 1 2 3 4 5 6 7New Edge 1-6 1-2 1-11 6-7 2-3 7-10 11-13New Vertex 6 2 11 7 3 10 13Iteration 8 9 10 11 12 13New Edge 13-14 3-9 2-4 11-12 12-8 4-5New Vertex 14 9 4 12 8 54.3.2 Reconfiguration as A Minimum Spanning Tree ProblemIn this section, the MST algorithm is applied, as a heuristic method, to the DS reconfigura-tion problem. This problem aims at finding a radial topology which admits the minimumresistive losses. It is stated in [112] that the minimum-loss network is achieved when allthe switches are closed. This statement, however, does not always hold. In order to showthis, the problem is solved without imposing the radiality constraint. The results of thesesimulations obtained by the MIQCP method proposed in Section 4.1.1 are reported in Sec-tion 4.3.4. The following assumptions are derived based on the engineering knowledge ofDS and are used throughout this paper. When all the tie-switches are closed, forming aweakly-meshed network, then• the losses are close (if not exactly the same) to the losses of the best-possible con-figuration.• the load is distributed between the feeders and a moderate balance (but not thebest-possible) among all the feeders is achieved.• the voltage profile of the network is close to the best-possible voltage profile.The validity of these assumptions is demonstrated through examples in Section 4.3.4.Note that loss reduction, load balancing between feeders, and voltage profile improvementare three different objectives and may lead to different configurations. However, lossreduction implicitly improves the voltage profile, as was also previously acknowledged in[44], and maintain a moderate balance between feeders, as was also recognized in [118].It is shown here, mathematically, that the loss minimization is in alignment withvoltage profile improvement. The total losses in a DS are obtained as:Ploss =∑(i,j)∈SLGij[(V rei − V rej )2 + (V imi − V imj )2](4.36)8712345678910111213140.07630.11080.11210.07430.04860.02340.02900.05970.02840.01270.00810.02980.03950.00640.02010.0254(a) Weighted 16-node graph1234567891011121314(b) MST of the 16-node graphFigure 4.7: The 16-node graph and its MST.where G is the branch conductance; SL is the set of all the branches; V re and V im arethe real and imaginary pars of the nodal voltages. The voltage at the substation node isusually considered to be a known quantity. Normally, the following values are assumed forthe substation voltage: V res = 1 and V ims = 0. Recall the fact from Section 4.1.1 that thenodal voltage magnitudes are mainly governed by their real parts, i.e. V re. Minimizing(4.36) can be translated into minimizing the differences between the real parts of the nodalvoltages. Since one of the voltages, i.e. V res , is already fixed at 1 p.u., the other nodal88voltages are also pushed to stay close to 1 p.u. to minimize the objective in (4.36). Thisargument reveals the fact that minimizing losses implicitly improves the voltage profile ofa distribution network.Let us look at a meshed DS as a weighted undirected graph, where the negative of thecurrent magnitudes are taken as the weights. The problem of finding a tree for the meshednetwork that closely imitates it can be converted into a MST problem. The MST algorithmtries to keep the branches with the largest current magnitudes (note the negative sign) andopen the branches with the least currents to form a tree. For example, take the 14-nodegraph shown in Fig. 4.7a, which is the 14-node system of [46], with the current magnitudesshown for each branch. Figure 4.7b shows that the branches with the least currents, i.e.5-14, 6-8, and 7-9 are open in the MST solution.Although the above discussed method appears reasonable, network constraints, such asnodal voltage limits and feeder ampacities have not been explicitly considered. However,it is shown in Section 4.3.4 that the system operational constraints are implicitly satisfiedin the final results. The flowchart of the proposed algorithm is shown in Fig. 4.8. Recallthat the MST algorithm proposed here is meant to quickly find a good initial solution tobe used in a direct mathematical approach, such as the MIQCP problem in Sections 4.1.1and 4.1.2. Most of the commercial solvers, e.g., CPLEX, have embedded heuristics tofind an initial solution to start the branch-and-cut procedure. These heuristics need toprovide a feasible solution in a short time. The proposed algorithm is to replace thosegeneral-purpose heuristics to quickly find a good initial solution.4.3.3 Series Neighborhood SearchIn order to further improve the optimality of the solution obtained by the MST algorithm, aLocal Search (LS) is performed based on BODFs derived in Section 3.4. Each new tie switchobtained by the MST algorithm breaks a particular loop in the graph. The simulationresults, discussed in the next section, show that the tie switches in the actual optimalsolution are mostly in the series neighborhood of the tie switches obtained by the MSTalgorithm. Two branches are called “Series Neighbors” if their common node is onlyconnecting those two branches. Mathematically, the degree of the common node mustbe two. Figure 4.9 shows examples of series neighbors (e.g., branches 3-4 and 4-5) andnon-series neighbors (e.g., branches 5-6 and 6-7). Note that one branch may have multipleseries neighbors. For instance, branches 3-4, 4-5, and 5-6 are all in series. By closingan open switch and opening its series neighbor, the radiality of the solution is preserved,while the optimality of the solution may improve. Note that the branch exchange betweentwo non-series neighbors may not necessarily lead to a tree and, hence, is not sought here.In order to reevaluate the losses after each branch exchange, a new power flow solutionis required. This is, however, computationally expensive to run a power flow after each89StartRun LPF formeshed networkRun MST usingcurrent magnitudesStart the LS usingthe 1st tie switchFind series neighborsfor the tie switchCalculate LODFfor series neighborsFind the besttie switchAny tieswitch left?Go to the nexttie switchFinal configurationNoYesFigure 4.8: Flowchart of the MST+LS algorithm for network reconfiguration.123 4 5 678Figure 4.9: Visual examples of series and non-series neighbors.branch exchange. The BODFs are used here to avoid running a new power flow solution. LetSTS represent the set of tie switches obtained by the MST algorithm with NTS members.Assume that each tie switch k in STS has nTSk series neighbors. Thus,∑NTSk=1 nTSk powerflow solutions are required to search for an improved solution. Instead, for every tie switch,the BODFs are calculated for its series neighbors assuming this particular switch is closedand all other tie switches are open. For instance, if 3-4 in Fig. 4.9 is a tie switch in theMST solution, both 4-5 and 5-6 are candidates for a branch exchange. The BODF for theoutage of these two branches (i.e. 4-5 and 5-6), assuming all other tie switches are open,90are calculated. The new current flow in each branch is then computed. Using the newcurrents, the new losses can be calculated as in (3.43). The BODF calculation routineneeds to be called only NTS times, whereas in the normal case the power flow routineneeds to be called∑NTSk=1 nTSk times. Also, each BODF evaluation needs less computationthan running a normal power flow. The size of the BODF calculated for each tie switchdepends on the number of series neighbors for that tie switch, which is obviously quitesmaller than the total number of branches in the network. The flowchart of the proposedalgorithm is shown in Fig. 4.8.4.3.4 Simulation ResultsIn this section, the differences between losses for the meshed network and the best-knownconfiguration are shown through several test cases to demonstrate the assumptions madein Section 4.3.1. Then, the proposed MST and LS algorithms are applied to sample distri-bution test systems to show their effectiveness.Meshed Network Versus Best-Known ConfigurationThe MIQCP formulation of Chapter 4.1.1 is used here to obtain the optimal configurationby dropping the radiality constraint. Several test systems are used and the problem issolved using CPLEX. Common test systems used in the literature are the 14-node [46],33-node [29], 70-node [91], 84-node [54], 119-node [98], 135-node [92], 415-node [59], 873-node [94], and 10476-node [94] systems. Simulation results for the meshed network and theoptimal solutions are provided in Table 4.10. As can be seen in Table 4.10, the differencebetween the losses in the meshed network (P1) and the optimal solution (P2) is negligible,which partially confirms the assumption in [112]. These results support the assumptionsstated in Section 4.3.1.Radial Configurations and Active LossesThe algorithm depicted in Fig. 4.8 is applied here to the test systems. The results ofthe analysis are reported in Table 4.11. The optimum solution for each case is obtainedusing the MIQCP method. In Table 4.11, 1 is the relative optimality gap for MST solution,defined as1 =|Popt − PMST|Popt× 100 (4.37)And 2 and 3 are defined similarly for the LS and DISTOP solutions. As can be seen, theoptimality gaps for the MST solutions are mostly less than 3%. For the 10476-node system,the optimal solution is not known to the authors. Nonetheless, it is possible to providea lower bound on the optimal solution using the losses in the meshed network, which is917838.7 kW. This gives a gap of 4%, which guarantees that the solution provided by theMST routine is within an optimality gap of less than 4%. For the 10476-node system, theproposed algorithm provides 47.4% loss reduction, and the same value is also reported in[47] for loss reduction. However, this reduction in losses is achieved in [47] in 472s. Thecomputation time of the proposed method is discussed later in this section.The column denoted by 2 in Table 4.11 shows the relative gap for the local search (LS)results. In all the cases, LS has led to an improved solution as compared to the MST results.The rightmost column of Table 4.11 shows the optimality gap for the DISTOP solution.Since there is no explanation on how to check the connectivity of the network at eachiteration, it was not possible for the authors to fully implement the DISTOP algorithm tocompare CPU times. It was also not possible to check the connectivity constraint manuallyfor the 10476-node system due to its dimensionality. Despite more computations requiredby this algorithm, the optimality of its solution is, in most of the cases, worse than theproposed MST+LS method in this paper.The tie switches in the MST solution and after applying the LS are listed in Table4.12. The branches that are different from the optimal solution are shown in bold. It iseasy to check that the branches in bold are electrically close to the ones in the optimalsolution. For example, in the 33-node system, to get the optimal solution Branch 28 hasto replaced by Branch 37, which have Node 29 in common, but are not in series. For the70-node system, the following two replacements gives the optimal solution: 39 → 70 and71→ 79. These pairs of branches have one node in common, but are not in series. Similardifferences were observed for other cases.The aforementioned results motivate the idea of focusing on the branches close to theMST+LS solution to find an improved solution in the branch-and-cut procedure. The initialsolution is used as the incumbent solution which helps prune all subproblems for which thevalue of the objective function is no better than the incumbent. Moreover, in the branch-and-cut process, it is possible to issue priority orders for variables in the branching strategy.Variables with higher priorities will be branched on first. The knowledge obtained fromthe MST+LS algorithm about the variables that may provide an improved solution canbe invoked here to assign priority orders. An appropriate variable ordering can result inspeed up in the solution process by branching on more influential variables [103].Voltage ProfileIt is interesting to look into the voltage magnitudes in the initial configuration and theradial configuration found by the MST method. For illustrative purposes, the 119-nodeand 873-node systems are considered here. The voltage profiles for these test systemsbefore and after reconfiguration are shown in Fig. 4.10. The minimum and average valuesfor the voltage magnitudes are also shown in the figure. As can be seen, a significant92Table 4.10: Comparison of The Best Possible Configuration andMeshed NetworkTest Case P1(kW)1 P2(kW)2 ∆P (%)3 Open Lines14-node 427.814 427.814 0 -33-node 123.291 123.253 0.03 970-node 198.425 198.355 0.03 51, 75, 7884-node 462.682 459.824 0.62 13, 63, 83, 84119-node 819.359 818.342 0.12 23, 51, 75, 129, 130135-node 271.846 270.83 0.37 84, 106, 135, 138873-node 990.336 986.346 0.40 132, 241, 412, 451, 596,639, 882, 888, 890, 9001 Losses for meshed network.2 Losses for best-known configuration. 3 Relative difference.Table 4.11: Comparison of Network Losses Obtained by The Proposed Algo-rithms and Other References#Nodes Losses(kW) ρ(%)i 1ii 2iii 3ivMIQCP MST LS DISTOP14 468.3 468.3 468.3 468.3 8.9 0 0 033 139.5 140.7 139.9 139.5 31.0 0.8 0.2 070 201.4 208.7 203.6 203.9 10.5 3.6 1.1 1.284 469.9 471.7 470.08 471.4 11.3 0.38 0.04 0.32119 869.7 894.3 883.5 891.9 31.9 2.8 1.5 2.5135 280.2 289.4 286.4 295.9 10.6 3.3 2.2 5.6415 2352.4 2358.6 2355.3 2357.2 11.4 0.26 0.12 0.20873 991.3 1001.9 1000.3 991.8 70.1 0.92 0.76 0.0510476 N/A 8173.4 8158.1 - 47.5 N/A N/A N/Ai Loss reduction by MST+LS w.r.t the initial configuration.ii Relative optimality gap for the MST solution (%).iii Relative optimality gap for the LS solution (%).iv Relative optimality gap for the DISTOP solution (%).93Table 4.12: Radial Configurations Obtained by The LS Algorithm#Nodes Tie Switches14 7,8,1633 7,9,14,28,3270 14,30,39,45,51,66,71,75,76,77,7884 7,34,39,42,63,72,83,84,86,89,90,92119 23,26,34,39,42,52,58,70,73,75,95,109,122,129,130135 9,35,51,54,90,96,106,126,135,136,138,141,143,144,145,146,147,148,150,151,1554157,13,34,39,42,63,72,83,84,86,89,90,92,103,109,130,135,138,159,168,179,180,182,185,186,188,199,225,231,234,255,264,274,276,278,280,281,282,284,295,321,327,330,351,360,370,372,374,376,377,378,380,391,417,423,426,447,456,466,468,470,472,473,474,476,481,482,483,484,485,486,487,488873 84,130,144,159,190,282,288,306,330,409,412,434,451,494,596,616,629,631,637,698,818,843,885,888,890,896,900improvement in voltage profiles has been achieved. These observations support the ideathat the radial network obtained by the MST routine provides improved voltage profiles.4.3.5 Feeders LoadingAlthough it is not its explicit objective, the minimum-loss network reconfiguration nor-mally results in a moderate load balance among all the feeders [118]. However, there arecases in which a good load balance cannot be achieved due to the increment in the lossesthat would be imposed. Another case is when there is no connection between one clusterof feeders to the other cluster. In this case, every cluster can be optimized independentlybecause they do not have any interdependencies. For instance, for the 84-node system in[54], Feeders F1, F7, and F8 are in one cluster and the rest are in another cluster. Nobranch exchange can be performed between these two clusters to reduce losses or improvethe load balance. Therefore, maintaining a perfect load balance may not be achievable insome cases. For illustrative purposes, the feeder currents are shown in Fig. 4.11 for the84-node and the 873-node systems, when the network is reconfigured for loss reductionusing the MST+LS method. As can be seen, the MST+LS solution generates fairly closecurrents to the ones in the optimal solution.It is of importance to notice that the main objective of the reconfiguration here is tominimize the losses, not to find the best-possible load balance among the feeders. In orderto balance the load and minimize the losses simultaneously, a multi-objective optimization9410 20 30 40 50 60 70 80 90 100 1100.860.880.90.920.940.960.981NodesNodalVoltages(p.u.)InitialMST+LS(a) The 119-node system100 200 300 400 500 600 700 8000.90.920.940.960.981NodesNodalVoltages(p.u.)InitialMST+LS(b) The 873-node systemFigure 4.10: Voltage profiles of the 119-node and 873-node systems under initialand obtained configurations.problem needs to be solved. The solution provided by the MST+LS routine is also a goodinitial solution to start such multi-objective optimizations and reduce the computationalburden of the branch-and-cut procedure.4.3.6 Computational ComplexityThe time complexity of Prim’s algorithm to find the MST is O(e log v). The MST algorithmand the LS algorithm were implemented in the MATLAB platform on an Intel Core i7-2600 CPU @ 3400 MHz with 8GB of RAM. The average CPU times for each test system is95F1 F2 F3 F4 F5 F6 F7 F8 F9 F10 F110.10.20.30.40.50.6Current(p.u.)MST+LSOptimal(a) The 84-node systemF1 F2 F3 F4 F5 F6 F70.10.20.30.40.5Current(p.u.)MST+LSOptimal(b) The 873-node systemFigure 4.11: Feeders currents for the MST+LS and the optimal solutions.divided into the MST time and LS time, as shown in Table 4.13. The computation times ofthe proposed method are, by orders of magnitude, smaller than many previously proposedheuristics. Also, the scalability of the MST method is a significant advantage. There areonly a few references in the literature that use systems with over 10,000 nodes to demon-strate the applicability of their method on large-scale systems. The CPU times in Table4.13 shows that the proposed method is able to find the solution for a system with over10,000 nodes in about a second, where to achieve the same optimality, the algorithm in[47] requires 472s.The speed up in the DS reconfiguration problem is achieved by sacrificing some possibleoptimality of the solution. However, this is a fast and robust algorithm which providesa good initial solution to start the more sophisticated algorithms, such as the MIP algo-96Table 4.13: Comparison of CPU Times for Different MethodsSystem MST LS Reference 1 Reference 233-node 0.015 0.022 7.2 ([53]) 0.3 ([50])70-node 0.021 0.085 3 ([91]) 2.4 ([49])84-node 0.023 0.151 36.15 ([54]) 0.3 ([50])119-node 0.023 0.261 8.61 ([53]) 9.04 ([98])135-node 0.026 0.551 7.3 ([59]) 0.4 ([50])873-node 0.048 4.85 20.9 ([47]) 874 ([60])10476-node 1.360 37.11 472 ([47]) -rithms proposed in Sections 4.1.1 and 4.1.2, as well as [59] and [60], to find the optimalsolution. Based on our experience, initializing the mixed-integer programming solvers us-ing the solution found by the MST algorithm can reduce the total solution time up to 40%,depending on the system under study.97Chapter 5Volt-VAR OptimizationAt the distribution system level, the deployment of smart meters has provided enormousamount of data on load behavior. Since the basic functionality of most of the commercialsmart meters involves capturing the waveforms of current and voltage, and conductingcalculations based on those waveforms, there is a significant opportunity to leverage theavailable waveforms to recognize the load. Many algorithms have been developed for loaddisaggregation at the smart meter level, such as the “Eigenload” approach proposed inChapter 2. Once these algorithms are implemented, the type of connected appliances ateach customer terminal is available at any time instant. Knowing the type of load, one canestimate the voltage dependence of the aggregated load at each distribution transformerlevel. In an extensive study conducted by the Electric Power Research Institution (EPRI),the voltage dependence characteristics of major appliances are derived based on labora-tory measurements [83]. Having both the load types and the voltage dependence of eachtype, an accurate-enough estimation of the aggregated load voltage dependence can beperformed.In this chapter, a framework for implementing VVO in modern distribution systemsis proposed. This framework takes into account the voltage dependence of loads, utilizesthe LPF formulation, and formulates the operation planning problem as a MIQCP problemthat can be solved efficiently to its global optimum. The analysis is truncated at thedistribution transformer level, which is the primary network, and the secondary networkis modeled as an aggregated load with a known composition. Capacitor banks, voltageregulators, and ULTC transformers are considered as the control variables and the problemis formulated for a day-ahead period.The rest of this chapter is organized as follows. In Section 5.1, the load modelingprocedure is described. The VVO algorithm is explained in Section 5.2. Simulation resultsare provided in Section 5.3.985.1 Load Characterization For Volt-VAR OptimizationTraditionally, loads are roughly categorized as residential, commercial, and industrial.For each load type, some typical characteristics are then derived based on historical mea-surements. For today’s distribution system, however, this level of decomposition is notsufficient. The smart meters installed at the customer level are capable of capturing loadvariations and even recognizing the load types. The load composition at each meter canthen be aggregated at the distribution transformer level to create an estimate of the totalload composition at the secondary network fed by each distribution transformer. Withoutloss of generality, a residential area is considered here. The number and types of on-lineappliances vary during a 24-hour period. The pattern for weekdays are quite similar,while the weekends have different patterns than the weekdays. The load composition alsochanges according to weather conditions, season, geographical location, etc. This impliesthe need for load modeling and estimation schemes for every single feeder during a periodof time to be able to create a relatively accurate load model and forecast.There are many different appliances in a typical residential home. Nonetheless, for thepurpose of load modeling for applications such as VVO, some appliances may be ignoreddue to their minor impact on the total demand response. For instance, a phone charger,which is typically rated at 5 W, does not influence the demand response of a house with 2kW total demand. This fact also relieves the accuracy requirement on the disaggregationscheme implemented at smart meters. In addition, there are appliances with relativelyhigh power demand which are only used for a very short period of time. For instance, amicrowave oven, which may be rated at 1.5 kW to 3 kW, is typically used for only a fewminutes to warm up foods. These types of load are also of no importance for the VVOapplication, assuming an action time window of, say, 30 minutes. Taking into account thementioned points, the number of appliances that need to be detected and properly mod-eled is reduced to a manageable extent. A list of typical major appliances for a residentialcustomer are given in Table 5.1 along with a load tag to be used later in this chapter.Note that this list is generated according to a typical cold-climate household and it mayvary for different countries and/or climates.Residential customers can be further divided into two categories: 1) with residentsleaving the place for work during the day, and coming back during the afternoon/evening(House 1); 2) with residents staying at home most of the time (House 2). Besides, aroundthe residential areas, it is very common to find small businesses such as restaurants andconvenience stores. Also, street lights are connected to a designated distribution trans-former on the same feeder. This information is obtained by studying a typical residentialarea in Vancouver, BC, Canada, and sample data is generated for the purpose of simula-tion. For the four types of customers mentioned earlier: House 1, House 2, Restaurant,99Table 5.1: List of Typical Appliances for a Residential AreaType Water Heater Space Heater RefrigeratorTag No. 1 No. 2 No. 3Type Electronics Incand. Light Fluor. LightTag No. 4 No. 5 No. 6Type Dishwasher Washing Machine Clothes DryerTag No. 7 No. 8 No. 9Type Electric Range Vacuum Cleaner Ventilation-FanTag No. 10 No. 11 No. 12Type Induction Cooker HP Sodium Light Air ConditionerTag No. 13 No. 14 No. 15and Street Light, a typical load composition over a 24-hour period during a weekday isshown in Figs. 5.1a-5.1d. A combination of these loads may be connected to a distribu-tion transformer. For thermostatic loads, e.g. refrigerators, the ON/OFF time intervalsvary from customer to customer so that, at an aggregated level, there are always a fewrefrigerators on-line. In a more sophisticated version of this framework, every appliancemay be from different manufacturer, and may exhibit different voltage dependence char-acteristics. For instance, the fluorescent lamps may be driven by an electronic or magneticballast [119], which may exhibit different voltage-power characteristics. In this chapter, asimplified version of load types is assumed for illustration purposes.The voltage dependence of a variety of appliances were derived in [83] using extensivelaboratory measurements. For the sake of completeness, a sample of the data is reproducedhere for some typical devices. Figure 5.2 depicts the per-unit values of active and reactivepower consumptions of the devices for a range of 0.9 p.u. to 1.05 p.u. of the terminalvoltage. Once the total demand for each type of device at a distribution transformer isknown, the corresponding graph in Fig. 5.2 is multiplied by those values and a voltage-dependent model for that load is generated. This is done for all other types of devices andthe sum of all the scaled graphs will provide the voltage dependence characteristic of thetotal load connected at a particular distribution transformer at a specific time.The outcome of the above analysis is to have an estimation of the parameters in (2.16)and (2.17) for each node in the system to be incorporated in the power flow analysis.5.2 Volt-VAR Optimization AlgorithmThe problem of distribution system VVO is, in some ways, similar to the unit commitmentproblem in the transmission system. In the unit commitment problem, the status of each1000 2 4 6 8 10 12 14 16 18 20 2202468Time (h)P(kW)(a) House 10 2 4 6 8 10 12 14 16 18 20 22024Time (h)P(kW)(b) House 20 2 4 6 8 10 12 14 16 18 20 220102030Time (h)P(kW)(c) Restaurant0 2 4 6 8 10 12 14 16 18 20 2202040Time (h)P(kW)No. 1 No. 2 No. 3 No. 4 No. 5 No. 6 No. 7No. 8 No. 9 No.10 No.11 No.12 No.13 No.14(d) Street Light (100×HP Sodium light)Figure 5.1: Typical load profile for a residential area in a weekday.1010.9 0.92 0.94 0.96 0.98 1 1.02 1.040.80.911.1V (p.u.)P(p.u.)No. 1 No. 3No. 4 No. 5No. 8 No. 9(a) Active Power0.9 0.92 0.94 0.96 0.98 1 1.02 1.040.811.2V (p.u.)Q(p.u.)No. 1 No. 3No. 4 No. 5No. 8 No. 9(b) Reactive PowerFigure 5.2: Voltage dependence of typical appliances [83].unit is to be determined during a certain period of time, say for a day-ahead, to minimizethe total generation cost. Similarly, in the VVO problem, the statuses of capacitor banks,voltage regulators, and ULTCs are to be determined for a day-ahead period, assuming aforecast for the demand, in order to satisfy objectives such as minimizing network losses.This is particularly important because of the fact that the number of switching actions hasto be limited. Switch operations will reduce the lifetime of the equipment. Therefore, theproblem cannot be decoupled for each time interval, i.e. the solution at the current timestep is dependent on the previous and next time steps. This characteristic requires theproblem to be formulated for the whole period of the planning horizon which, in turn, dras-tically increases the problem dimensionality. To cope with this problem, simplificationsare sought here to reduce the dimensionality while respecting a certain level of accuracy.Once the status of all system equipment are determined, minor adjustments can bemade in near-real-time to compensate for possible errors in the predictions and/or to re-spond to a contingency. This mimics very well the concept of optimal power flow for102real-time adjustments, based on the unit commitment schedule. For instance, equipmentwith continuous adjustment capabilities such as distributed generators, distributed energystorages, distribution static compensators, and interruptible loads are perfect candidatesfor near-real-time adjustments, especially for reactive power support. In addition, the sta-tus of the mentioned equipment does not need to be determined ahead of time since theseare almost always connected to the gird. In this study, the problem of day-ahead schedul-ing of switchable equipment is studied and the problem of near-real-time adjustments isleft as future work of the research group.5.2.1 System ModelingDistribution systems (DS) are known for the unbalance of the phase currents, which stemsfrom the non-uniform distribution of single-phase loads at the household level. To runan accurate load flow, a three-phase algorithm is required. However, for some specificapplications, it may be accurate-enough to model the single-phase loads as balanced three-phase loads. Most of the capacitor banks, voltage regulators, and ULTCs are three-phaseunits, i.e. a switching action involves all three phases simultaneously. In addition, all theanalysis done in this chapter deals with the primary distribution network, which includesthe LV network between the substation and the secondary distribution transformers. Thesecharacteristics make it a valid approximation to model primary distribution system as abalanced network.The LPF formulation presented in (4.1) and (4.2) is adopted here. A capacitor bankis modeled here as an admittance element B¯c, which then appears in the correspondingdiagonal elements of B¯, i.e.:Bci = ui B¯ci (5.1)Here ui ∈ {0, 1} represents the status of the capacitor bank connected at Node i.A transformer with ULTC capability is modeled using the pi-equivalent circuit [120], asshown in Fig. 5.3. In this figure, a stands for the transformer ratio and yt is the transformeradmittance for a = 1. The two shunts on both sides compensate for the different tap ratiosto model the voltage buck/boost. For the sake of simplicity, the transformer impedanceis assumed to be inductive. In order to linearize the term a (a − 1) for 0.9 ≤ a ≤ 1.1,it is approximated by a − 0.9963 which places a maximum error of 0.08% in the voltagemagnitudes in the simulated systems. It should be noted that a is a discrete variable andcan only take values between 0.9 and 1.1 with steps of, e.g., 0.01, which yields 20 steps intotal. This leads to the following mathematical representation:a =20∑j=1bjxj ,20∑j=1xj = 1 (5.2)103a (a− 1) yta yt(1− a) ytVm VkFigure 5.3: The pi-equivalent circuit of a ULTC-equipped transformer.in which xj ∈ {0, 1} and bj ∈ {0.9, 0.91, . . . , 1.1}. The last equality constraint in (5.2) is aSOS which can be efficiently handled by commercial mixed-integer programming solvers.A multiplication of binary-continuous variables may be encountered by adding (5.1)or (5.2) to the modified admittance matrix elements in (4.1) and (4.2). This can be dealtwith similar to what was done in Chapter 4.1.1 using (4.7).The VVO problem may have different objectives, such as loss reduction, voltage profileimprovement, and CVR. These objectives as well as the corresponding system operationalconstraints are described in the following subsections.Active Loss ReductionThe active power losses are calculated as:Ploss =∑m