UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

EMTP-based load disaggregation at feeder levels Shojae, Shima 2015

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2015_november_shojae_shima.pdf [ 4.71MB ]
Metadata
JSON: 24-1.0166755.json
JSON-LD: 24-1.0166755-ld.json
RDF/XML (Pretty): 24-1.0166755-rdf.xml
RDF/JSON: 24-1.0166755-rdf.json
Turtle: 24-1.0166755-turtle.txt
N-Triples: 24-1.0166755-rdf-ntriples.txt
Original Record: 24-1.0166755-source.json
Full Text
24-1.0166755-fulltext.txt
Citation
24-1.0166755.ris

Full Text

EMTP-based Load Disaggregation at Feeder LevelsbyShima ShojaeB.Sc., Sharif University of Technology, 2011A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFMaster of Applied ScienceinTHE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES(Electrical and Computer Engineering)The University of British Columbia(Vancouver)September 2015c© Shima Shojae, 2015AbstractPower distribution networks play an important role in electricity grid. Distribution sys-tem components require becoming smarter and more automated for the sake of improv-ing their reliability and increasing their operational efficiency. Smart meters are one ofthe powerful devices that achieved this goal. However, their data are of minimal use —grid or load information obtained from smart meters are shallowly analysed. This thesistakes advantage of the shortcoming by accurately calculating the load information usingEMTP-based load disaggregation method. The approach is applicable to residential loadsat small scale and feeders at large scale.In this thesis we first give our theoretical method for load disaggregation inspired byEMTP computational program. Then with simulation and experimental results, we demon-strate that our work outperforms past solutions by the following advantages:1. EMTP-based load disaggregation is applicable at every point of interest, i.e., fromdistribution feeder down to the customer’s entry point.2. Unlike other methods, our method employs both transient and steady state loadproperties.3. Last but not least, our solution is capable of determining load’s electrical parame-ters.In this thesis we stress on three major eigen-loads: (1) motor, (2) resistive, and (3) purelyinductive. Then we report how much of the load is made of each eigen-load. We exam-pled our method on a number of PSCAD simulation cases and a few real appliance mea-surements. Our results prove load disaggregation shall assist power system engineers inevaluating the power flow on accurate load observations.iiPrefaceThe originality of this thesis is based on the author’s authentic research work. All theory,simulation, and experimental results have been produced in UBC Electrical EngineeringPower Systems Lab exercised within the last two years.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xivGlossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xixAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxi1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Power Distribution Networks . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1.1 An Example of Distribution: BC-Hydro Distribution Network . . . . 31.2 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41.2.1 Energy Management . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51.2.2 Smart Grid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61.2.3 Voltage VAR Optimization (VVO) . . . . . . . . . . . . . . . . . . . . 81.2.4 Linear Power Flow (LPF) . . . . . . . . . . . . . . . . . . . . . . . . . 92 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.1 Appliance Load Monitoring . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.1.1 Intrusive Load Monitoring (ILM) . . . . . . . . . . . . . . . . . . . . . 122.1.2 Non Intrusive Load Monitoring (NILM) . . . . . . . . . . . . . . . . . 132.2 Load Disaggregation Elements . . . . . . . . . . . . . . . . . . . . . . . . . . 14iv2.2.1 Data Acquisition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.2.2 Feature Extraction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.2.3 Load Classification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.3 Literature Reviews of Previous Studies . . . . . . . . . . . . . . . . . . . . . 202.4 Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 323.1 EMTP Software . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 333.2 PSCAD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 343.3 Discretization Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 353.3.1 Forward Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 363.3.2 Backward Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 363.3.3 Trapezoidal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 373.3.4 Simpson’s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 383.4 The Thesis Method Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . 443.4.1 Load’s Electrical Circuit Order Determination Knowing the Struc-ture of the Load . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 463.4.2 Load’s Electrical Circuit Order Determination Unknowing the Struc-ture of the Load . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 463.4.3 Examples of Parametric Equivalent RLC Circuits Trapezoidal Models 503.5 Network Synthesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 563.5.1 Network Synthesis Description . . . . . . . . . . . . . . . . . . . . . . 563.6 Foster Forms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 573.7 Cauer Forms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 583.8 RL-RC Network Synthesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 623.8.1 Realization of ZRC in Foster 1 Form: . . . . . . . . . . . . . . . . . . . 623.8.2 Realization of ZRC in Foster 2 Form: . . . . . . . . . . . . . . . . . . . 623.8.3 Realization of ZRC in Cauer 1 Form: . . . . . . . . . . . . . . . . . . . 633.8.4 Realization of ZRC in Cauer 2 Form: . . . . . . . . . . . . . . . . . . . 643.8.5 Realization of ZRL in Foster 1 Form: . . . . . . . . . . . . . . . . . . . 643.8.6 Realization of ZRL in Foster 2 Form: . . . . . . . . . . . . . . . . . . . 653.8.7 Realization of ZRL in Cauer 1 Form: . . . . . . . . . . . . . . . . . . . 653.8.8 Realization of ZRL in Cauer 2 Form: . . . . . . . . . . . . . . . . . . . 66v4 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 674.1 Discrete-Continues Transformation . . . . . . . . . . . . . . . . . . . . . . . . 694.1.1 Bilinear Transform Method . . . . . . . . . . . . . . . . . . . . . . . . 694.2 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 704.2.1 Simulation One: Parallel Resistive and Inductive Load . . . . . . . . 704.2.2 Simulation One: Paralleled RL Load Parameters Calculations . . . . 734.2.3 Simulation Two: Resistor in Series with a Parallel Inductance andCapacitance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 744.2.4 Simulation Three: Motor Load . . . . . . . . . . . . . . . . . . . . . . 834.2.5 Simulation Three: Single Motor Load Parameters Calculations . . . . 914.2.6 Simulation Four: Paralleled Resistive and a Motor . . . . . . . . . . . 924.2.7 Simulation One: Paralleled Motor and Resistance Load ParametersCalculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 994.2.8 Simulation Five: Single Motor in Parallel with One Resistance andOne Inductance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1004.2.9 Simulation Five: Single Motor in Parallel with One Resistance andOne Inductance Parameters Calculations . . . . . . . . . . . . . . . . 1044.2.10 Simulation Six: Two Series Motor in Parallel with a Resistive and anInductive Load . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1064.2.11 Simulation Seven: Three Series Motor in Parallel with a Resistiveand an Inductive Load . . . . . . . . . . . . . . . . . . . . . . . . . . . 1124.2.12 Simulation Eight: Two Parallel Motors Loads . . . . . . . . . . . . . . 1154.2.13 Simulation Eight: Two Parallel Motors Load Parameters Calculations 1234.2.14 Simulation Nine: Two Parallel Motors in Parallel with One ResistiveLoad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1254.2.15 Simulation Nine: Two Parallel Motors in Parallel with A ResistiveLoad Parameters Calculations . . . . . . . . . . . . . . . . . . . . . . 1294.2.16 Simulation Ten: Two Parallel Motors in Parallel with A Resistiveand An Inductive Load Parameters Calculations . . . . . . . . . . . . 1325 Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1355.1 Appliances Electrical Load Types and Classification . . . . . . . . . . . . . . 1355.1.1 Resistive Loads . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1355.1.2 Inductive Loads . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136vi5.1.3 Motor Loads . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1365.1.4 Non-Linear Loads . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1365.2 Measuring General Observations . . . . . . . . . . . . . . . . . . . . . . . . . 1375.2.1 Fan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1375.2.2 Fan Parameters Calculation . . . . . . . . . . . . . . . . . . . . . . . . 1385.2.3 Florescent Lamp . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1395.2.4 Florescent Parameters Calculation . . . . . . . . . . . . . . . . . . . . 1395.2.5 Fan in Parallel With Florescent Lamp Parameters Calculation . . . . 1405.2.6 Parallel Florescent and Fan Parameters Calculation . . . . . . . . . . 1405.2.7 Hair Dryer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1415.2.8 HTC Charger . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1425.2.9 Refrigerator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1435.3 Chapter Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1446 Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1456.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1456.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147viiList of TablesTable 3.1 Different loads order prediction. . . . . . . . . . . . . . . . . . . . . . . . . 55Table 3.2 Foster and Cauer impedance/admittance expansions . . . . . . . . . . . 61Table 4.1 Values of current at time t and t− ∆t, voltage at time t− ∆t for the firstsampling scenario in the case of a parallel RL load. . . . . . . . . . . . . . 72Table 4.2 Values of voltage at time t for the first sampling scenario in the case of aparallel RL load. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72Table 4.3 Values of current at time t and t−∆t, voltage at time t−∆t for the secondsampling scenario in the case of a parallel RL load. . . . . . . . . . . . . . 72Table 4.4 Values of voltage at time t for the second sampling scenario in the caseof a parallel RL load. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73Table 4.5 Values of current at time t and t− ∆t, voltage at time t− ∆t for the firstsampling scenario in the case of a resistance in series with a paralleledLC load. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76Table 4.6 Values of voltage at time t for the first sampling scenario in the case of aresistance in series with a paralleled LC load. . . . . . . . . . . . . . . . . 76Table 4.7 Values of current at time t and t−∆t, voltage at time t−∆t for the secondsampling scenario in the case of a resistance in series with a paralleledLC load. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76Table 4.8 Values of voltage at time t for the first sampling scenario in the case of aresistance in series with a paralleled LC load. . . . . . . . . . . . . . . . . 77Table 4.9 Values of current at time t, t− ∆t and t− 2∆t, voltage at time t− ∆t andt− 2∆t for the first sampling scenario in the case of a resistance in serieswith a paralleled LC load. . . . . . . . . . . . . . . . . . . . . . . . . . . . 77Table 4.10 Values of voltage at time t for the first sampling scenario in the case of aparallel RL load. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78viiiTable 4.11 Values of current at time t, t− ∆t and t− 2∆t, voltage at time t− ∆t andt − 2∆t for the second sampling scenario in the case of a resistance inseries with a paralleled LC load. . . . . . . . . . . . . . . . . . . . . . . . . 78Table 4.12 Values of voltage at time t for the second sampling scenario in the caseof a resistance in series with a paralleled LC load. . . . . . . . . . . . . . . 79Table 4.13 Values of current at time t, t − ∆t, t − 2∆t and t − 3∆t, voltage at timet− ∆t, t− 2∆t and t− 3∆t for the first sampling scenario in the case of aresistance in series with a paralleled LC load. . . . . . . . . . . . . . . . . 79Table 4.14 Values of voltage at time t for the first sampling scenario in the case of aresistance in series with a paralleled LC load. . . . . . . . . . . . . . . . . 80Table 4.15 Values of current at time t, t − ∆t, t − 2∆t and t − 3∆t, voltage at timet− ∆t, t− 2∆t and t− 3∆t for the second sampling scenario in the caseof a resistance in series with a paralleled LC load. . . . . . . . . . . . . . . 80Table 4.16 Values of voltage at time t for the second sampling scenario in the caseof a resistance in series with a paralleled LC load. . . . . . . . . . . . . . . 81Table 4.17 Values of current at time t and t− ∆t, voltage at time t− ∆t for the firstsampling scenario for a single motor load . . . . . . . . . . . . . . . . . . 84Table 4.18 Values of voltage at time t for the first sampling scenario for a singlemotor load. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84Table 4.19 Values of current at time t and t−∆t, voltage at time t−∆t for the secondsampling scenario for a single motor load . . . . . . . . . . . . . . . . . . 85Table 4.20 Values of voltage at time t for the second sampling scenario for a singlemotor load. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85Table 4.21 Values of current at time t and t− ∆t, voltage at time t− ∆t for the thirdsampling scenario for a single motor load . . . . . . . . . . . . . . . . . . 85Table 4.22 Values of voltage at time t for the third sampling scenario for a singlemotor load. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86Table 4.23 Values of current at time t, t− ∆t and t− 2∆t, voltage at time t− ∆t andt− 2∆t for the first sampling scenario for a single motor load. . . . . . . . 86Table 4.24 Values of voltage at time t for the first sampling scenario in the case of asingle motor load . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87Table 4.25 Values of current at time t, t− ∆t and t− 2∆t, voltage at time t− ∆t andt− 2∆t for the second sampling scenario for a single motor load. . . . . . 87ixTable 4.26 Values of voltage at time t for the second sampling scenario in the caseof a single motor load . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88Table 4.27 Values of current at time t, t− ∆t and t− 2∆t, voltage at time t− ∆t andt−−2∆t for the third sampling scenario for a single motor load. . . . . . 88Table 4.28 Values of voltage at time t for the third sampling scenario in the case ofa single motor load . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89Table 4.29 Values of current at time t and t− ∆t, voltage at time t− ∆t for the firstsampling scenario in the case of a paralleled resistance and a motor. . . . 93Table 4.30 Values of voltage at time t for the first sampling scenario in the case of aparalleled resistance and a motor. . . . . . . . . . . . . . . . . . . . . . . . 93Table 4.31 Values of current at time t and t−∆t, voltage at time t−∆t for the secondsampling scenario in the case of a paralleled resistance and a motor. . . . 93Table 4.32 Values of voltage at time t for the second sampling scenario in the caseof a paralleled resistance and a motor. . . . . . . . . . . . . . . . . . . . . 94Table 4.33 Values of current at time t and t− ∆t, voltage at time t− ∆t for the thirdsampling scenario in the case of a paralleled resistance and a motor. . . . 94Table 4.34 Values of voltage at time t for the third sampling scenario in the case ofa paralleled resistance and a motor. . . . . . . . . . . . . . . . . . . . . . . 94Table 4.35 Values of current at time t, t− ∆t and t− 2∆t, voltage at time t− ∆t andt− 2∆tfor the first sampling scenario in the case of a parallel resistanceand a motor resistance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95Table 4.36 Values of voltage at time t for the first sampling scenario in the case of aparalleled resistance and a motor . . . . . . . . . . . . . . . . . . . . . . . 95Table 4.37 Values of current at time t, t − ∆t and t − 2∆t, voltage at time t − ∆tand t − 2∆tfor the second sampling scenario in the case of a parallelresistance and a motor resistance. . . . . . . . . . . . . . . . . . . . . . . . 96Table 4.38 Values of voltage at time t for the second sampling scenario in the caseof a paralleled resistance and a motor. . . . . . . . . . . . . . . . . . . . . 96Table 4.39 Values of current at time t, t− ∆t and t− 2∆t, voltage at time t− ∆t andt− 2∆tfor the third sampling scenario in the case of a parallel resistanceand a motor resistance. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97Table 4.40 Values of voltage at time t for the third sampling scenario in the case ofa paralleled resistance and a motor. . . . . . . . . . . . . . . . . . . . . . . 97xTable 4.41 Values of current at time t, t − ∆t, t − 2∆t and t − 3∆t, voltage at timet− ∆t, t− 2∆t and t− 3∆t for the first sampling scenario in the case of amotor load in parallel with an inductance and a resistance. . . . . . . . . 101Table 4.42 Values of voltage at time t for the first sampling scenario in the case of amotor load in parallel with an inductance and a resistance. . . . . . . . . 101Table 4.43 Values of current at time t, t − ∆t, t − 2∆t and t − 3∆t, voltage at timet− ∆t, t− 2∆t and t− 3∆t for the second sampling scenario in the caseof a motor load in parallel with an inductance and a resistance. . . . . . . 102Table 4.44 Values of voltage at time t for the first sampling scenario in the case of amotor load in parallel with an inductance and a resistance. . . . . . . . . 102Table 4.45 Values of current at time t, t − ∆t, t − 2∆t and t − 3∆t, voltage at timet− ∆t, t− 2∆t and t− 3∆t for the first sampling scenario in the case oftwo series motor in parallel with a resistive and an inductive Load. . . . 107Table 4.46 Values of voltage at time t for the first sampling scenario in the case oftwo series motor in parallel with a resistive and an inductive Load. . . . 107Table 4.47 Values of current at time t, t − ∆t, t − 2∆t and t − 3∆t, voltage at timet− ∆t, t− 2∆t and t− 3∆t for the second sampling scenario in the caseof two series motor in parallel with a resistive and an inductive Load. . . 108Table 4.48 Values of voltage at time t for the second sampling scenario in the caseof two series motor in parallel with a resistive and an inductive Load. . . 108Table 4.49 Values of current at time t, t−∆t, t− 2∆t, t− 3∆t and t− 4∆t, voltage attime t−∆t, t− 2∆t, t− 3∆t and t− 4∆t for the first sampling scenario inthe case of two series motor in parallel with a resistive and an inductiveLoad. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109Table 4.50 Values of voltage at time t for the first sampling scenario in the case oftwo series motor in parallel with a resistive and an inductive Load. . . . 109Table 4.51 Values of current at time t, t− ∆t, t− 2∆t, t− 3∆t and t− 4∆t, voltageat time t− ∆t, t− 2∆t, t− 3∆t and t− 4∆t for the second sampling sce-nario in the case of two series motor in parallel with a resistive and aninductive Load. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110Table 4.52 Values of voltage at time t for the second sampling scenario in the caseof two series motor in parallel with a resistive and an inductive Load. . . 110xiTable 4.53 Values of current at time t, t− ∆t, t− 2∆t, t− 3∆t, t− 4∆t and t− 5∆t,voltage at time t− ∆t , t− 2∆t, t− 3∆t, t− 4∆t and t− 5∆t for the firstsampling scenario in the case of three series motors in parallel with aresistive and an inductive load . . . . . . . . . . . . . . . . . . . . . . . . . 113Table 4.54 Values of voltage at time t for the first sampling scenario in the case of athree series motors in parallel with a resistive and an inductive load . . . 113Table 4.55 Values of current at time t, t− ∆t, t− 2∆t, t− 3∆t, t− 4∆t and t− 5∆t,voltage at time t−∆t , t− 2∆t, t− 3∆t, t− 4∆t and t− 5∆t for the secondsampling scenario in the case of three series motors in parallel with aresistive and an inductive load . . . . . . . . . . . . . . . . . . . . . . . . . 114Table 4.56 Values of voltage at time t for the second sampling scenario in the caseof a three series motors in parallel with a resistive and an inductive load. 114Table 4.57 Values of current at time t, t − ∆t, t − 2∆t and t − 3∆t, voltage at timet− ∆t, t− 2∆t and t− 3∆t for the first sampling scenario in the case oftwo parallel motors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116Table 4.58 Values of voltage at time t for the first sampling scenario in the case oftwo parallel motors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116Table 4.59 Values of current at time t, t − ∆t, t − 2∆t and t − 3∆t, voltage at timet− ∆t, t− 2∆t and t− 3∆t for the second sampling scenario in the caseof two parallel motors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117Table 4.60 Values of voltage at time t for the second sampling scenario in the caseof two parallel motors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117Table 4.61 Values of current at time t, t−∆t, t− 2∆t and t− 3∆t and t− 4∆t, voltageat time t−∆t, t− 2∆t, t− 3∆t and t− 4∆t for the first sampling scenarioin the case of two parallel motors. . . . . . . . . . . . . . . . . . . . . . . . 118Table 4.62 Values of voltage at time t for the first sampling scenario in the case oftwo parallel motors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118Table 4.63 Values of current at time t, t− ∆t, t− 2∆t and t− 3∆t and t− 4∆t, volt-age at time t− ∆t, t− 2∆t, t− 3∆t and t− 4∆t for the second samplingscenario in the case of two parallel motors. . . . . . . . . . . . . . . . . . . 119Table 4.64 Values of voltage at time t for the second sampling scenario in the caseof two parallel motors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119xiiTable 4.65 Values of current at time t, t−∆t, t− 2∆t and t− 3∆t and t− 4∆t, voltageat time t−∆t, t− 2∆t, t− 3∆t and t− 4∆t for the third sampling scenarioin the case of two parallel motors. . . . . . . . . . . . . . . . . . . . . . . . 120Table 4.66 Values of voltage at time t for the third sampling scenario in the case oftwo parallel motors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120Table 4.67 Values of current at time t, t− ∆t, t− 2∆t, t− 3∆t and t− 4∆t, voltageat time t−∆t, t− 2∆t, t− 3∆t and t− 4∆t for the first sampling scenarioin the case of two parallel motors in parallel with a resistive load. . . . . 126Table 4.68 Values of voltage at time t for the second sampling scenario in the caseof two parallel motors in parallel with a resistive load. . . . . . . . . . . . 126Table 4.69 Values of current at time t, t−∆t, t− 2∆t, t− 3∆t and t− 4∆t, voltage attime t−∆t, t− 2∆t, t− 3∆t and t− 4∆t for the second sampling scenarioin the case of two parallel motors in parallel with a resistive load. . . . . 127Table 4.70 Values of voltage at time t for the second sampling scenario in the caseof two parallel motors in parallel with a resistive load. . . . . . . . . . . . 127Table 5.1 Experimental appliances list. . . . . . . . . . . . . . . . . . . . . . . . . . . 136xiiiList of FiguresFigure 1.1 Distribution network from substation to feeders from [2]. . . . . . . . . 3Figure 1.2 Typical distribution feeder schematic. . . . . . . . . . . . . . . . . . . . . 4Figure 1.3 World’s total primary energy supply from 1965 to 2035 by fuel [5]. . . . 5Figure 1.4 DMS (Distribution Management System) schematic [9]. . . . . . . . . . . 7Figure 1.5 ZIP-ZI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11Figure 2.1 Power consumption of a single family for 40 minutes [18]. . . . . . . . . 13Figure 2.2 Overview diagram of NILM [19]. . . . . . . . . . . . . . . . . . . . . . . . 14Figure 2.3 Load electrical features categorization [27]. . . . . . . . . . . . . . . . . . 16Figure 2.4 Load clustering using P and Q [16]. . . . . . . . . . . . . . . . . . . . . . 19Figure 2.5 Energy consumption during a day for a house [35]. . . . . . . . . . . . . 21Figure 2.6 REDD box instalment inside a house . . . . . . . . . . . . . . . . . . . . . 21Figure 2.7 Load three dimensional characteristics[36]. . . . . . . . . . . . . . . . . . 23Figure 2.8 Metered vs reconstructed consumption of a dish washer and refrigera-tor using probabilistic method [28]. . . . . . . . . . . . . . . . . . . . . . 24Figure 2.9 Graphical illustration of voltage-current trajectories from six differentappliances from REDD household 3 [26]. . . . . . . . . . . . . . . . . . . 25Figure 2.10 The STFT results on turn on current transients. (a) 160 horse power in-duction motor, (b) 123 horse power induction motor driven by variablevoltage driver. (c) A bank of loads supplied by a six-pulse thyristorrectifier for AC power [37]. . . . . . . . . . . . . . . . . . . . . . . . . . . 26Figure 2.11 The DWT results of turn on instantaneous power transients for a 160hp induction motor [37]. . . . . . . . . . . . . . . . . . . . . . . . . . . . 27Figure 2.12 The DWT results of turn on instantaneous power transients for a 123hp induction motor driven by variable voltage driver [37]. . . . . . . . . 28xivFigure 2.13 The DWT results of turn on instantaneous power transients for a bankof loads supplied by a six-pulse thyristor rectifier for AC power [37]. . . 29Figure 2.14 A computer and an incandescent light bulb in the ∆ P, ∆ Q and ∆ 3rdharmonic coordinate system [38]. . . . . . . . . . . . . . . . . . . . . . . . 30Figure 2.15 Frequency spectrum of a particular light switch being toggled. . . . . . . 31Figure 3.1 Different integration methods description. . . . . . . . . . . . . . . . . . 35Figure 3.2 Approximation of x(t) by Forward-Euler method [41]. . . . . . . . . . . . 36Figure 3.3 Approximation of x(t) by Backward Euler method [41]. . . . . . . . . . . 37Figure 3.4 Approximation of x(t) by Trapezoidal method [41]. . . . . . . . . . . . . 38Figure 3.5 Approximation of x(t) by Simpson method . . . . . . . . . . . . . . . . . 39Figure 3.6 Inductance equivalent EMTP model. . . . . . . . . . . . . . . . . . . . . . 40Figure 3.7 Voltage across the inductor [42]. . . . . . . . . . . . . . . . . . . . . . . . 41Figure 3.8 Capacitance equivalent EMTP model. . . . . . . . . . . . . . . . . . . . . 42Figure 3.9 Voltage across the capacitance [42]. . . . . . . . . . . . . . . . . . . . . . . 43Figure 3.10 Voltage across the resistance [42]. . . . . . . . . . . . . . . . . . . . . . . . 43Figure 3.11 Different length windows having the same initial point as one of “volt-age/current sampling points” selection method. . . . . . . . . . . . . . . 48Figure 3.12 Same length windows having different initial point as one of “volt-age/current sampling points” selection method. . . . . . . . . . . . . . . 48Figure 3.13 Proposed method: Process from smart-meter data capturing to loadidentification. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49Figure 3.14 Parallel RL case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50Figure 3.15 A parallel LC in series with a resistance case. . . . . . . . . . . . . . . . . 51Figure 3.16 Induction motor type. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53Figure 3.17 First Foster form-partial fraction expansion in the impedance form. . . . 57Figure 3.18 Second Foster form-partial fraction expansion in the admittance form. . 58Figure 3.19 Ladder network topology [52]. . . . . . . . . . . . . . . . . . . . . . . . . 59Figure 3.20 First Cauer form continued fraction expansion about infinity. . . . . . . 60Figure 3.21 Second Cauer form continued fraction expansion about zero. . . . . . . 60Figure 3.22 First Foster, 2nd Foster, first Cauer and 2nd Cauer forms method sum-mary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61Figure 3.23 Realization of ZRC in F1 form [53]. . . . . . . . . . . . . . . . . . . . . . . 62Figure 3.24 Realization of ZRC in F2 form [53]. . . . . . . . . . . . . . . . . . . . . . . 63xvFigure 3.25 Realization of ZRC in C1 form [53]. . . . . . . . . . . . . . . . . . . . . . . 63Figure 3.26 Realization of ZRC in C2 form [53]. . . . . . . . . . . . . . . . . . . . . . . 64Figure 3.27 Realization of ZRL in F1 form [53]. . . . . . . . . . . . . . . . . . . . . . . 64Figure 3.28 Realization of ZRL in F2 form [53]. . . . . . . . . . . . . . . . . . . . . . . 65Figure 3.29 Realization of ZRL in C1 form [53]. . . . . . . . . . . . . . . . . . . . . . . 65Figure 3.30 Realization of ZRL in C2 form [53]. . . . . . . . . . . . . . . . . . . . . . . 66Figure 4.1 Mapping of the s-plane to the z-plane with the bilinear transform method. 69Figure 4.2 Inductive load paralleled with a resistive load. . . . . . . . . . . . . . . . 71Figure 4.3 First order impedance coefficients for the first and second scenarios ofa paralleled RL load. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73Figure 4.4 Paralleled RL Load . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73Figure 4.5 Paralleled RL parameters calculation . . . . . . . . . . . . . . . . . . . . . 74Figure 4.6 Resistive load in series with a paralleled LC. . . . . . . . . . . . . . . . . 75Figure 4.7 First order impedance coefficients for the first and second scenarios ofa resistance in series with a paralleled LC load. . . . . . . . . . . . . . . . 82Figure 4.8 Second order impedance coefficients for the first and second scenariosof a resistance in series with a paralleled LC load. . . . . . . . . . . . . . 82Figure 4.9 Third order impedance coefficients for the first and second scenarios ofa resistance in series with a paralleled LC load. . . . . . . . . . . . . . . . 83Figure 4.10 Motor load case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83Figure 4.11 First order impedance coefficients for three different scenarios of a sin-gle motor load. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90Figure 4.12 Second order impedance coefficients for three different scenarios of asingle motor load. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90Figure 4.13 Single motor load. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91Figure 4.14 Motor load parameter calculation. . . . . . . . . . . . . . . . . . . . . . . 91Figure 4.15 Paralleled resistive and a motor. . . . . . . . . . . . . . . . . . . . . . . . 92Figure 4.16 First order impedance coefficients for three scenarios of a paralleledmotor and resistance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98Figure 4.17 Second order impedance coefficients for three scenarios of a paralleledmotor and resistance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98Figure 4.18 Paralleled motor and resistance. . . . . . . . . . . . . . . . . . . . . . . . 99Figure 4.19 Motor load in parallel with a resistive load parameter calculation. . . . . 99xviFigure 4.20 Single motor in parallel with one resistance and one inductance load. . . 100Figure 4.21 Third order impedance coefficients for three different scenarios of a mo-tor load in parallel with an inductance and a resistance. . . . . . . . . . 103Figure 4.22 Fourth order impedance coefficients for three different scenarios of amotor load in parallel with an inductance and a resistance. . . . . . . . 104Figure 4.23 Single motor in parallel with one resistance and one inductance . . . . . 104Figure 4.24 Motor in parallel with an inductive and a resistive load parameter cal-culation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105Figure 4.25 Two series motor in parallel with a resistive and an inductive load. . . . 106Figure 4.26 Three series motor in parallel with resistive and inductive load. . . . . . 112Figure 4.27 Two parallel motors loads. . . . . . . . . . . . . . . . . . . . . . . . . . . . 115Figure 4.28 Third order impedance coefficients for two different sampling scenariosof two parallel motor loads. . . . . . . . . . . . . . . . . . . . . . . . . . . 121Figure 4.29 Fourth order impedance coefficients for three different sampling sce-narios of two parallel motor loads. . . . . . . . . . . . . . . . . . . . . . . 122Figure 4.30 Two parallel motors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123Figure 4.31 Two parallel motors load parameter calculation. . . . . . . . . . . . . . . 124Figure 4.32 Two parallel motor in parallel with one resistive load. . . . . . . . . . . . 125Figure 4.33 Third order impedance coefficients for two different sampling scenariosof two parallel motors in parallel with one resistance load. . . . . . . . . 128Figure 4.34 Fourth order impedance coefficients for two different sampling scenar-ios of two parallel motors in parallel with one resistance load. . . . . . . 128Figure 4.35 Fifth order impedance coefficients for two different sampling scenariosof two parallel motors in parallel with one resistance load. . . . . . . . . 129Figure 4.36 Two parallel motors in parallel with a resistive load. . . . . . . . . . . . . 129Figure 4.37 Two parallel motors load parameter calculation. . . . . . . . . . . . . . . 131Figure 4.38 Two Parallel motors in parallel with a resistive and an inductive load. . 132Figure 4.39 Two parallel motors in parallel with a series RL load parameter calcu-lation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133Figure 4.40 General overview of or work. . . . . . . . . . . . . . . . . . . . . . . . . . 134Figure 5.1 Voltage measured at main outlet. . . . . . . . . . . . . . . . . . . . . . . . 137Figure 5.2 Current waveform in time for a fan. . . . . . . . . . . . . . . . . . . . . . 138Figure 5.3 Current waveform in time for a florescent lamp . . . . . . . . . . . . . . 139xviiFigure 5.4 Current waveform in time for a hair-dryer. . . . . . . . . . . . . . . . . . 142Figure 5.5 Current waveform in time for an HTC laptop charger. . . . . . . . . . . . 143Figure 5.6 Current waveform in time for a refrigerator. . . . . . . . . . . . . . . . . 143xviiiGlossaryThe following are abbreviations and acronyms used in this thesis, listed in alphabeticalorder:SCADA supervisory control and data acquisitionAGC Automatic Generation ControlDMS Distribution Management SystemVVO Voltage VAR OptimizationSG Smart GridDS Distribution SystemOLTC On Load Tap ChangerLPF Linear Power FlowALM Appliance Load MonitoringNILM Non Intrusive Load MonitoringNALM Non Intrusive Appliance MonitoringILM Intrusive Load MonitoringTHD Total Harmonic DistortionREDD Residential Energy Disaggregation DataSTFT Short Time Fourier TransformxixSVM Support Vector MachineANN Artificial Neural NetworkDWT Discrete Wavelet TransformLPF Linear Power FlowSFA Shifted Frequency AnalysisAMI Advanced Metering InfrastructurexxAcknowledgementsFirst, I would like to thank my supervisor, Professor Dr. Jose Marti. I had such an honourto meet Dr. Marti when for the first time, I visited UBC in 2013. I am super grateful to haveDr. Marti recruiting me, especially during hard times of this research. He never stoppedmotivating me yet encouraged to work harder. I remember whenever I was stuck in un-charted areas of the work, he kept telling that, there always has to be a solution. He wasright, there was a solution. Without his constructive feedback, technical support, patientguidance, and his big heart, I would have not been able to finish this work.Second, I would like to thank my graduate course instructors, in particular, Dr. HermannDommel, Dr. Nagpal Mukesh, Mr. Nathan Ozog, and Dr. Ebrahim Vaahedi. I learned alot in their classes. They opened up my eyes in the power system world, the challengesand the solutions. They provided me technical insight.Third, my thanks go to UBC power systems lab graduate students. We all were like afamily working together and pushing each other towards reaching our goals. In partic-ular, Dr. Hamed Ahmadi, I had a close collaboration with him during my research. Healways spot me with difficult questions, which triggered to think deeper to find genuineanswers.Fourth, my parents, the angles of my life. Whoever I am today and I will be in future,undoubtedly is because of them. My father was my role model from the beginning of myeducation. I learned his hard working while inherited his sharpness and patience. Due tohis treat, I have always been number one in school to enter Sharif. He taught me to workhard for our goals to achieve them. My mom, the symbol of scarification who supportedme emotionally during every stage of my life. I am indebted and proud of these two won-derful people for the whole rest of my life.Lastly and the most importantly, many many thanks to my handsome and kindheartedMehdi. He is the one who wants me have the best, makes me laugh hard and be thexxiluckiest girl. Mehdi always pushes for stretching the boundaries to get closer to my bigdreams. Without him, this thesis could have never come together. During the past sixmonths, he did a lot for me specificity patiently walking with me for writing my thesis.xxiiChapter 1IntroductionEvery power systems network is comprised of three main sections: (1) Generation, (2)Transmission, and (3) Distribution. Generation takes care of producing electrical powerfrom the potential water accumulated behind dams or burning the coal or natural gas.Transmission is the process of transferring the generated power to feeders — electricalpower is transmitted at high voltage through overhead transmission lines. Distribution isabout dividing the electrical power to lower voltages to supply residential and industrialcustomers. In this thesis target domain is distribution. Distribution networks, as the lastpart of a power system chain are reportedly less studied [1], thus require more attention.This chapter gives the background overview of power distribution network then the mainmotivation behind the thesis work. In the following, Section 1.1, an overview of the powersystem network is provided, then in Section 1.2, the problem statement and fundamentalneeds for this thesis are discussed.1.1 Power Distribution NetworksPower distribution compared to generation and transmission is less studied [1]. This iswhile distribution networks play a key role in distributing the energy where monitoringthe quality of power and voltage is significantly important. However adding monitorsin distribution networks carry some challenges. Traditional distribution networks, con-sist of feeder elements such as reclosers, switches, fault indicators, capacitor banks, andvoltage regulators that are incapable of sending critical information to the control center.Besides, the major obstacle is the lack of communication infrastructure across the distribu-1tion systems. There basically is no or little possible communication way between controlcentres and distribution network equipment. There certainly is a need for automationinvolvement in distribution similar to transmission network — in transmission, faults’propagation is managed by monitoring, localizing, and predicting the fault through au-tomation.To overcome the shortage, utilities launched deployment of smart meters. Smart meters,collect customers’ consumption information regularly. Smart meters emergence, sincelast decade, encourages most or all sections of a network to become automated and in-telligent. This thesis will cover more on smart meters in Section 1.2.2. However, morework needs to be done to have a smart, automated and efficient Distribution Manage-ment System (DMS). An automated DMS will fulfil the following requirements for futuregrids:1. Integration of distributed renewable energy resources in substitute of traditionalfossil fuels.2. Providing management for optimising demand, i.e., providing incentives for cus-tomers to use electrical energy in a more efficient manner.3. Improving Automatic Generation Control (AGC) system due to the peak demandrise, especially for emerging demand of electrical vehicles.4. Managing fast-response behaviour of power electronics components of the grid ver-sus slow-response nature of synchronous generators.Additionally, industry has developed the following smart tools for main parts of adistribution system. Some of these smart tools are:1. Voltage VAR Optimization (VVO) : VVO algorithm, optimises the voltage profilealong the distribution feeders and customers to avoid voltage flickers and fluctua-tions. VVO is not practical without being informed about description of the load.2. Stochastic load and generation forecasting in consideration of large renewable pen-etration.3. Optimal network reconfiguration: Providing an optimum switching manoeuvre forthe feeder. Disconnecting some parts of the feeder compromised to reduce the en-ergy loss and supply more critical loads. Load disaggregation reveals descriptionsof feeder connected loads beneficial to network reconfiguration.2Having less accurate real time measurement devices, not accessing to the efficient algo-rithms corresponding to the real time data management, and possessing no accurate loadmodels restrict the distribution system’s future improvement. To improve the operationalbehaviour of the distribution network, feeder planning, outage prediction and customers’reliability improvement are well advised. Figure 1.1 shows a basic schematic of a smartdistribution network. Figure 1.2 shows a typical radial distribution feeder layout.Figure 1.1: Distribution network from substation to feeders from [2].In following section, BC-Hydro distribution network as an example is discussed.1.1.1 An Example of Distribution: BC-Hydro Distribution NetworkA typical BC-Hydro distribution network consists of an infinite bus, which is connectedvia a primary distribution line to the distribution substation. The low voltage bus at thesubstation is generally regulated with one of the following elements:1. Automatic load tap changer (LTC) with a set point of about 124 V (7.44/12.89 kV or14.88/25.77 kV).2. 3-phase 300/400 Amp feeder position voltage regulator.3. Bus regulator with a set point of 122-123 Volt (120.0 Volt on a secondary basis is3Figure 1.2: Typical distribution feeder schematic.equivalent to Bc-Hydro nominal primary voltages of 7.2/12.47 kV and 14.4/25.0kV) [3].To this end, we introduced the structure of a typical distribution network. We pointed outthe difference between traditional distribution network and a smart distribution networkas well.Distribution system dynamic behaviour is all about the nature of its loads. Load disag-gregation function is to separate the aggregated load to its components. Knowing thecomponents of the total load, supports power engineers towards a more efficient powersystem’s operation and maintenance.1.2 MotivationThe basic goal of a power system cycle is to supply different types of the load. Thus,having an accurate information about customer loads plays a critical role in the distribu-tion system planning and operation. On the other hand, advent of distributed generatorssuch as solar panels, onshore/offshore wind farms and PVs forces conventional grid to4operate more dynamically, thus capturing load’s dynamic behaviour turns into be morecrucial matter. In the following Sections 1.2.1 and 1.2.4, the most important applicationneed of load disaggregation are discussed.1.2.1 Energy ManagementToday, energy management is challenging due to the fast growing nature of energy de-mands (loads). Electricity demand is increasing twice as fast as overall energy use. Onthe other hand, energy supplies are declining to two thirds by 2035 as shown in Fig-ure 1.3. This increasing need and decreasing energy production, coupled with the factthat coal fuels are 40% of energy supply contributor, makes electricity energy generationas the highest cause of CO2 emission. Hence, world needs to increase electricity energysupply preferably clean energy resources for the next 20 years [4]. Accessing adequateFigure 1.3: World’s total primary energy supply from 1965 to 2035 by fuel [5].clean energy supply resources without tackling demand side management is not practi-cal. Thereby, to cope with the growing energy demand while producing less greenhousegasses, it is important to understand the type of loads connected to the grid as the main5players’ demand type. Real-time information helps out in having a less lossy distributionnetwork and a more stable and reliable transmission system thus saving energy.1.2.2 Smart GridSmart grid refers to a modernized and smart healing network. Smart Grid (SG) answersmost of the upcoming challenges which Distribution System (DS) is facing since 20th cen-tury. SG adds up communication and automation features to grid’s elements [6]. Strongcommunication capability enables utilities, substations, and distribution customers mon-itor each other in real-time. For example, if a contingency occurs and one generator ortransmission line trips, distribution system is able to get updated in less than a secondreducing outage duration.Smart MeterSmart grid development, brought smart meters and in-home energy displays as the solu-tion to address the issue of accessing customers’ real-time demands information. Smartmeters are intelligent energy meters that obtain customer load’s information and send itback to control centre. Meter measures time domain voltage and current waveforms fora customer and sends the obtained data for evaluation to supplier [7]. Received real-timedata are analysed for the following purposes:1. Giving customers an insight of their energy consumption.2. It offers tariff models in electricity market.3. Providing customers with incentives about the time of usage.4. Enabling the utilities to monitor the grid’s aggregated load and energy in real time.5. Helping out distribution engineers in optimising VVO techniques.Load disaggregation exploits smart meter readings to provide customers with an intuitivequantitative understanding of the amount of energy consumed by them. The main aim ofload disaggregation is energy saving.Distribution Management SystemWith the addition of millions of new data resources in DS, it is imperative for utility op-erators to monitor and control a high volume of real-time data in an optimum manner.6DMS is the decision support tool, which collects data all over from smart measurementdevices. It takes the data back to the control centre. Operators assess the grid based onreceived data. The most important application of DMS is that, it improves operators anddistribution engineers’ visibility to the distribution system. DMS is valuable because ofits application in the following:1. Advanced control of the voltage and reducing network’s loss through the automaticdistribution optimized reconfiguration.2. Integration of distributed renewable resources.3. Integration of dynamic loads and having dynamic feeder protection.In other words, DMS is an effective mechanism of managing the AMI data and perform-ing complex analysis and calculation [8]. Figure 1.4 shows a large picture of a DMS net-work including all subsystems and devices incorporated together. Load disaggregationFigure 1.4: DMS (Distribution Management System) schematic [9].algorithm has to be embodied in a DMS, towards achieving a better performance from allsections of the DMS platform.71.2.3 VVOVVO DefinitionVVO is defined as maintaining the voltage level at a standard range (+/- 15% nominalvoltage according to CSA standard) through the distribution network. VVO has been ofgreat interest by utilities for many years. Load’s behaviour and consequently voltage/-VAR level are changing throughout the day. Therefore, distribution engineers need toapply different methods to control the voltage at the feeders. Deployment of voltage regu-lators and shunt capacitors are the two common methods of controlling the voltage withina limit. Voltage and reactive power have a proportional relation, meaning that increasingthe reactive power leads to voltage rise and reactive power drop will cause voltage dropas well. Because of this, inductive and capacitive loads play an important role in con-trolling the voltage profile of the system. Recently, with advent of Advanced MeteringInfrastructure (AMI), engineers try to leverage real time data achieved from smart metersto control the voltage along the feeders more realistically. Voltage regulators, capacitorbanks and On Load Tap Changer (OLTC) are common equipment employed, to controlthe voltage and the reactive power. OLTC increases and decreases the voltage level ac-cording to the load changes. For example, if load increase across the feeder, tap changerswill increase the voltage to account for the excessive voltage drop on the feeder. Voltageregulators and capacitor banks perform the same.Load Disaggregation Role in VVO DefinitionLoads play an important role in VVO techniques. Motor loads form more than 80% ofthe total load. They have inductive characteristics. They absorb Q from the feeder whichleads to a voltage drop. Also, for a capacitive load case, they generate reactive powerwhich may cause a voltage rise in the system. Thus, loads nature has a direct impacton the system voltage and power levels. Being able to detect the accurate type of loadhelps VVO to control the system more accurately.supervisory control and data acquisition (SCADA) enables optimising the voltage andreactive power, according to the load changes.All in all, we need load voltage and power information at the feeder level to adjust thevoltage and reactive power along a feeder [10].81.2.4 Linear Power Flow (LPF)Power flow calculation is one of the strongest tools to analyse distribution network. Powerflow study solves a group of non-linear equations to calculate the bus voltage magnitudesand phase angles. Advance distribution system automation should be fast enough tomeet the real time demand response. Linearising the power flow algorithm gives thisopportunity to analysis the distribution system faster and more efficiently. Linear PowerFlow (LPF) is introduced for the first time in [11]. There are different load models that areused in solving a power flow problem as described in following section.Load ModelsLoads modelling has not been taken into serious consideration for the last 2 decades inpower systems area [12]. Based on BC-Hydro experimental results, there is 0.6 percentKWh (active power P) reduction per 1 percent load’s voltage reduction, and a 3 percentKVAR (Reactive power Q) reduction per 1 percent voltage reduction which brightens thesignificance of Loads models [13]. Loads are classified to two types of load modellingmethodologies:1. Static: So far, static load models are integrated in a power flow solution. Static loadmodels are not dependant on time, i.e. static models express P&Q relationship tovoltage at the same time instant while dynamic approaches describe P and Q asa function of voltage at each time instant. Power system study tools mostly usepolynomial and exponential load models under the static group.2. Dynamic: Dynamic loads characteristics change in time. So, they should be mod-elled in real time.For static models, Equation 1.1 formulates the voltage dependency of load characteristicsby an exponential model.P(V) = P0(VV0)aQ(V) = Q0(VV0)b(1.1)Polynomial models (ZIP model) also are widely used for load’s behaviour formulation.There are three major components in (ZIP model), (1) constant power, (2) constant current,9(3) constant impedance. They are formulated in Equation 1.2.P(V)P0= Fz(VV0)2 + FI(VV0) + FpQ(V)Q0= (F′z)(VV0)2 + F′I(VV0) + F′p(1.2)Where, FZ,FI ,FP are constant impedance, constant current and constant power respec-tively. Both polynomial and exponential models add non-linearity to the power flowequations. A new static load model is proposed in [11], which demonstrates a linearpresentation of load’s behaviours. Equation 1.3 describes the new load model (ZI model).It assumes a load, as the combination of constant impedance ( Z) and constant current (I),which leads to a linear power flow formulation.P(V)P0= Cz(VV0)2 + CI(VV0)Q(V)Q0= (C′z)(VV0)2 + C′I(VV0)(1.3)Where (CI + CZ = 1)and(C′I + C′Z) = 1.CI ,CZ,C′1,C′Z are derived from measurement data. A curve fitting will determine the val-ues of C and C’. It is a matter of solving a fitting optimization problem as formulated inEquation 1.4.minimizeN∑i=1(Cz(Vi)2 + CIVi − P(Vi))2 (1.4)subject to Cz + CI = 1. Figure 1.5 shows a comparison between three different explainedload modelling forms, applied on a three phase induction motor with the ratings of: 460V, 3-phase, 1.4 hp, 1725 r/min .Load Disaggregation Application in LPFLoads are behaved as voltage dependent instead of constant P and Q as described inEquatuion 1.3. Voltage dependency characteristics of the loads are required for solvinga linear power flow (LPF). Different type of loads have different voltage dependency in-dexes. Therefore, we should know the type of a load to estimate its voltage dependency10Figure 1.5: Comparison of the exponential, ZIP and proposed LPF load mod-els. [11].Data in this figure are from [14].behaviour. Load disaggregation result is fed into the LPF problem. EPRI ( Electric PowerResearch Institute) conducted a study on how to derive voltage dependency of differenttype of loads from laboratory measurements [15]. Above mentioned factors motivated usto devise an innovative tool to recognize grid connected types of loads. The main focus ofthis thesis is load disaggregation.11Chapter 2Literature ReviewThis chapter provides previous work made in research and industry. In Section 2.1, basedon surveyed literature, appliance load monitoring as the first topic discussion for loaddisaggregation is discussed. In Section 2.2, based on a collection of papers, load disag-gregation main elements are listed and the contributions in each domain are reviewed.Finally, in Section 2.3, popular studies in the area of load disaggregation are discussed.2.1 Appliance Load MonitoringAppliance Load Monitoring (ALM) is essential for energy management. It allows to ob-tain appliance specific energy consumption. It is approachable in the following two meth-ods:1. Intrusive Load Monitoring (ILM) which is known as distributed sensing method.2. Non Intrusive Load Monitoring (NILM) which is known as single point sensingmethod [16].In the following, these two methods are reviewed.2.1.1 ILMILM requires deploying smart power outlets on every appliance due to extra hardwarecost and installation complexity. However, it is accurate in measuring appliance’s specificenergy consumption.122.1.2 NILMNILM, also called, load disaggregation, is an area in computational sustainability thattries to discern what electrical loads (i.e. appliance) are running within a physical area.Such areas can include communities, industrial sites, office towers, buildings, homes, oreven within an appliance [17].NILM receives the aggregated load data from home entry electrical panel with no needto install sensors for each appliance inside the home. George W. Hart was the first onecoming up with this idea about load disaggregation [18].NILM HistoryNILM is initiated almost two decades ago by Hart [18]. He suggested to install Non In-trusive Appliance Monitoring (NALM) devices connected to the total load circuit. Basedon switching on/off events and specific signatures for each individual load, he was ableto determine the nature and number of loads connected. Moreover, he was able to checkout each load’s active and reactive power consumption and its time of the day variation.Figure 2.1 shows the power consumption versus the time of a single family for 40 minutesFigure 2.1: Power consumption of a single family for 40 minutes [18].period. There are four step changes in the active power consumption. This states existenceof four different appliances. However, with knowing each individual load’s consumption13pattern, we can say when each load is switched on or off and for how long. Figure 4.2Figure 2.2: Overview diagram of NILM [19].presents NILM as a block diagram which is fed in with two main inputs, measurementsand appliance data information. NILM connects with the total load using standard meteruser socket. It measures information such as voltage, current, P, Q and Total HarmonicDistortion (THD) for total load of a house entry. Appliance data block contains specificappliance’s information such as P signatures, Q signatures and also type of the applianceobtained by the manufacturer [19]. The output of NILM process is the daily status of theanalysed appliance.2.2 Load Disaggregation ElementsThis section elaborates three main elements of a load disaggregation process: (1) data ac-quisition, (2) feature extraction and (3) learning/classification.Researchers working in load disaggregation domain, argue about different techniquesand state of the art algorithms which are applicable in any of the three major modules.Most researchers consider the active and reactive power patterns as the appliance distinctfeatures, to distinguish between different types of home appliances. The common prop-erty among most of the load disaggregation studies, is to report the types of loads lookingat a single customer’s point of entry. Their algorithm will state which appliances turn onor off at each moment.142.2.1 Data AcquisitionData Acquisition receives loads electrical waveforms measured at smart meters. Differentenergy meters have different sampling rates. We need to choose the appropriate samplingfrequency based on the type of signal chosen for load’s feature. For example, if steadystate real power and reactive power are the metrics, low sampling rate (i.e., 120 Hz) isadequate. But If harmonics are studied as loads signatures, we need to sample the mea-sured data at a higher frequency rate. For instance, to capture the harmonics up to 6, theminimum sampling frequency should be 360 Hz. The other important point to consider isNyquist-Shannon criteria. Nyquist theorem states that sampling frequency has to be morethan twice of the highest harmonic frequency in a signal. Since there is no higher than11th harmonics available, maximum sampling rate is between 1.5-02 KHz [20]. Higherfrequency meters are useful in capturing more accurate information of a signal to reachhigh accuracy of load detection. They are able to capture “Microscopic” features suchas harmonics and instantaneous signal values. In compare with Microscopic signatures,power consumption changes are “Macroscopic” features, which can be detected with lowfrequency meters.2.2.2 Feature ExtractionAfter obtaining load data, we have to choose a distinguishable feature of the load. Rawdata received from meter needs to be processed towards achieving a load feature metric,to distinguish between different types of loads. Loads’ signatures are classified into steadystate and transient categories.1. Steady state features: RMS values of a waveform like voltage and current and anevent change in a metric like power consumption are examples of steady-state fea-tures. Sampling time of equal or larger than 1 second in adequate to extract steadystate features. For example power consumption changes usually happens wheneveran appliance turns on or off. This happens in more than one second time duration.Capturing this type of slow changing load properties does not need high volumeof memory or fast paced measurement devices. This is known to be the advantageof using steady state load features as far as they are able to differentiate betweendifferent types of loads. The other dominance property is that they do not requirerepeatable patterns unlike transients signals [18] and [21].152. Transients features: An appliance’s starting current, starting power, shape and du-ration of high order harmonics and also some other transient load’s characteristicsthat happen in less than one milliseconds are examples of transient features. Tran-sient features recognize the non-linear type of loads. Interestingly, transient signalscontain some unique characteristics of the load since, some characteristics, happenin high frequencies [22] [23] [24] [25] [26]. Figure 2.3 demonstrates different cate-gories of load features.Figure 2.3: Load electrical features categorization [27].Selecting an appropriate method either time domain approaches or frequency domain,depends on the level of disaggregation accuracy required and also the meter type. Forexample, if harmonics are chosen as appliances’ distinctive features, we employ frequencydomain analysis (fft, dft).2.2.3 Load ClassificationLearning process is grouped into two categories:1. Supervised Learning Approaches (Classification)2. Unsupervised Learning Approaches (Clustering)16Supervised LearningSupervised classification algorithms, assign a label for each input data. Given a trainingdata in the form of xi which represents ith example data, and yi which is the ith class label,the algorithm finds a model in which A(xi) = yi. The extracted load features are analysedwith an optimization method, or with a pattern recognition (machine learning) technique.Researches in this domain, usually pattern recognition technique to train a seen library ofdata and detect the unseen load based on trained data information. Learning systemallows training a classifier based on available load measurements signals. Any patternrecognition composes of three major steps:1. Selecting a category of load features which are distinctive enough to classify be-tween different types of loads.2. Synthesizing a training set of feature’s data.3. Training a classier with the observed data towards identifying the new type of loadafter training.Supervised learning uses KNN, Byes theorem, SVM, Neural Networks and Optimizationmethods as the classification tools to make the distinctive clusters of different load types.In following, each of above mentioned classifiers are briefly introduced.KNN KNN works based on the smallest distance between the input data and trainedone. Accordingly the class label assigned to smallest k, is given to that input.Byes Theorem Byes theorem assigns the class which has the highest probability is as-signed to the observed data [21] [28]. We can assign class c to an unknown examplewithin feature sets: X = (x1, · · · ,xn) such that: c = maxc P(C = c | x1, · · · ,xN).The P probability can be formulated as : P(C = c | x1, · · · ,xN) = P(C=c).P(x1,··· ,xN |C=c)P(x1,··· ,xN)SVM SVM is mainly used for binary classification. SVM method separates an ‘n’ di-mensional space into two classes. However, since data set is not always linearly separa-ble, sometimes kernel-induced feature space is introduced to cast the space into higherdimensional one [23] [25] [26].17Neural Network Neural network works based on three layers of trained neurons: (inputlayer, hidden layer and output layer). Observed data set is fed into the input layer. Dif-ferent weights are assigned to the input data due to a function which could be designedbased on problem’s objective. Input data it passed through the hidden layer. Finally, Out-put layer is a value which shows the class label of the corresponding input. To illustratemore on this one,in order to identify an event associated with operation of an new appli-ance, we match data set which includes all available event based appliances’ signatures,with new load’s event based value. There are many literature arguing different availablerecognition techniques [22] [26] [29] [30].Optimization methods Optimization methods consider load disaggregation problem as anoptimization problem. It works well when we need to identify an unknown single loadwhich is not present in the data set. It compares the new measured feature vector likep(t), to the feature vectors which exist in the data set library. The objective is to reduce thematching error. The optimization problem is formulated in Equation 2.1class = argmini|| y∧i − yi || (2.1)In which yi is the new measured feature vector and y′i is the feature vector which exists inthe load library. The limitation appears when, when we are working with an unknown ag-gregated load data. The reason arises form the fact that, an optimization problem can hasone objective function at each time instant. To identify the components of an aggregatedload, the method should consider all possible combinations of the appliances present inthe data set, which are probable to synthesis the unknown signal. The combination, whichleads to the least matching error will be selected as the components of the total load. Inte-ger programming and generic algorithms are the examples of an optimization problem’ssolvers [31]. However, in the presence of a complex load feature vector, it is getting com-putationally expensive and time consuming to approach the load disaggregation as anoptimization problem.Un-Supervised LearningRecently, researchers are developing unsupervised learning techniques to resolve someof supervised learning algorithm’s drawbacks. Supervised learning system is expensiveand not efficient due to the fact that preparing a complete data set covering all the possible18types of loads, needs new meters and sometimes sensors’ installation. Unlike supervisedmethods, unsupervised algorithms, tackles the load disaggregation, in a smarter way withno need to the data set pool. K-means and clustering approaches, are two examples ofan unsupervised learning method [32] [33] [34]. For example in Figure 2.4 P and Q areselected as the main load features to constitute the major load clusters. The goal behindFigure 2.4: Load clustering using P and Q [16].any unsupervised learning approach, is to achieve the minimum distance between anunknown load and the existing clusters in order to recognize, which main category ofloads, the new load belongs to.Most of unsupervised load disaggregation methods, need to full fill following criteria dueto their application with smart meters and due to having a more simplified computationeffort.1. Power consumption data should be captured in real time.2. No training dataset is necessary.3. Number of total appliances usually does not exceed 20-30.4. These methods cover ON/OFF appliances, multi state appliances, continuous con-suming appliances and permanent consuming appliances.5. Methods are scalable in a sense of robustness.192.3 Literature Reviews of Previous StudiesHart adopted active power as the main signature of loads. NILM’s task is to decomposethe P value in Equation 2.2, into its single load’s power components [16]. Total value ofactive power mainly is formulated as the sum of all the individual load’s consumption asfollowing:[P(t) = p1(t) + p2(t) + ...+ pn(t)] (2.2)Where pi is the power consumption of individual appliance and n is the total number ofappliances. There are literatures which employs unsupervised learning methods to cap-ture a load’s type. Dominik Egarter, takes advantage of HMM (Hidden Markov Model)for the appliance modelling and FHMM (Factorial Hidden Markov Model) for the house-hold consumption modelling. All the hidden states of appliances, which are generatedwith HMM, are transferred to FHMM. FHMM aggregates the appliance’s power con-sumption. He proposes popular Monte Carlo (Partial Filtering) estimation method, toestimate the disaggregated appliance’s states. Partial Filtering calculates the weightedparticles to obtain PDF values for each state. Particles are propagated over time to ob-tain the new particles and new the weights [32]. For the validation, this method is, testedon a household synthetic power draw numbers. It was capable of recognizing up to 18different appliances. He evaluated his approach on Residential Energy DisaggregationData (REDD) data set.REDD contains detailed power usage information from several homes [35]. Figure 2.7demonstrates how to install a REDD box inside a house.We conclude that unsupervised learning methods, are advantageous because they donot need any training library. They require a general knowledge about the appliancestructures such as, the number of operations ( whether on/off and multi state appliances)and also, power demand measured values of each appliance.In [33] load disaggregation is viewed as a single channel source separation problem.In this study, tensor (Multi-way array) factorization algorithm is applied for the electricalsource modelling. Multi-way array considers the power consumption of each load ateach smart home as a tensor. The Goal is to achieve all appliance’s significant informationthrough the available measurements data resulted from tensors factorization. Authors inthis paper intend to develop a model based on a non-negative sparse coding, in particulara non-negative matrix factorization algorithm, which enforces the sparsity conditions to20Figure 2.5: An example of energy consumption during a day for one of the house [35].Figure 2.6: REDD box instalment inside a housethe load models. This approach learns a sparse model for each load, which is later used forthe purpose of disaggregation. This method assumes an aggregated power consumptionsignal x(t), corresponding to the sum of all of the appliance’s usage over time period of T:x(t) = f (x1(t), ...,xk(t)) (2.3)21In Equation 2.3, f is assumed as a linear combination of constitution load. In this case wecan formulate f as :x(t) =∑ xi(t) (2.4)To solve a single source separation problem, all the load’s models are used to extract thesignificant characteristics of each source xi. Load’s models can be learned using On-linemeasurement data if available. Otherwise Matrix Factorization is a common method forsource modelling. xi(t) is modelled as sum of the products of the main characteristic’smatrix and activation matrix as formulated in Equation 2.5.xi(t) =∑ a(t) ∗ b(t) (2.5)In which xi(t) is the power consumption of the ith appliance at time instant t. a(t) containsmain features of appliance i and matrix b is the activation matrix. Figure 2.7 shows thethree way tensor. Note that, different matrices x1 to xk across the devices direction, aresource models for appliances 1 to k. This decomposition extracts information about thethree directions ( time of the day, days of the week, type of electrical appliance). RMSE(Root-Mean-Square-Error) is the evaluation metric in this approach which measures theerror between predicted and actual measured power consumption as formulated in Equa-tion 2.6.The restriction of this model is that each source is independently modelled, ignoring thedependency characteristics existing between different appliances. In other words, if twodifferent appliances are turned on simultaneously, the algorithm, is not able to distinguishbetween them.RMSE(X¯, ˆ¯X) =√√√√√ T∑t=1 m∑d=1(X¯− ˆ¯X)2T ∗m (2.6)In which X¯ is the actual aggregated signal and ˆ¯X is the predicted value. m is the num-ber of days for training. The other popular approach to load disaggregation problem isadopting probabilistic methods.Kim Etal was the first person to use probabilistic method for solving load disaggrega-22Figure 2.7: Load three dimensional characteristics[36].tion problem. He applied FHMM as the main method. He also, considered modellingthe coupling between dependant appliances. His method, incorporates a non-geometricparametric distribution of discrete-time ON durations. Even though his algorithm takesadvantage of the time-ON statistics information to address the overlapping problem, butit still suffers form the computational complexity which increases exponentially with thenumber of appliances [34]. In [28] Zeifman, applied probabilistic approaches to meet thesmart meter’s requirements as following:1. Load’s selected features have to be compatible with the smart meters.2. Algorithm is able to detect near real time appliances status.3. Disaggregation accuracy should not drastically drop when the number of appliancesis increasing.In this paper, Semi-Markov model is selected as the mathematical solution for load dis-aggregation problem. The candidate of appliances are chosen between the one which hasthe overlap in the power draw for the sake of decreasing the complexity of problem. Theidea is, to consider a group of two appliances i and i+1, with different time off durations.The probability of that an observed power change corresponds to the appliance i, dependson whether the appliance i+1 is turned on or off and for how long. Briefly explained, thealgorithm is able to determine the starting time and the ending time of an appliance, withhaving an aggregated metered data. For example, Figure 2.8 shows the measured powerconsumption versus the reconstructed power pattern of a refrigerator and a dish-washerusing Probabilistic method.23Figure 2.8: Metered vs reconstructed consumption of a dish washer and refrigeratorusing probabilistic method [28].In article [26], researches use “V-I” trajectory, (the mutual locus of instantaneous volt-age and current waveforms) towards separating an aggregated power data to its eigen-loads as shown in Figure 2.9.They also tested different numbers of learning algorithms such as feed forward ArtificialNeural Network (ANN), Hybrid learning algorithm and Support Vector Machine (SVM)to evaluate the precision and robustness of different classification algorithms. Energyconsumption data is needed for this method to extract instantaneous voltage and currentwaveforms information. This method is considered as the high frequency approach sinceit needs to sample voltage and current signals at a rate of 100 samples per second. Ex-tracted data is for the purpose of training which means, system learns different load’sbehaviours based on these sampling points.This method, suffers from its requirement for a huge learning data set. It also, needsclustering methods, to pre-process the dataset in order to apply V-I trajectory methodon it . Jin-Wen University researchers, proposed analysis of load’s transient operationalcharacteristics named as transient time and transient energy to detect power demand andload operation. Study applies DWT (Discrete Wavelet Transform) to the frequency-timedomain features vectors to analyse the behaviour of a specific type of loads. This researchasserts that using steady state features is less accurate in compare with employing tran-sient features of the loads. Because, in the event based methods, (steady state information24Figure 2.9: Graphical illustration of voltage-current trajectories from six different ap-pliances from REDD household 3 [26].are captured after an event like switching ON/OFF happens)whenever two different ap-pliances have the same switching time or when the power of a certain load equals the totalpower of the other load, it is not possible to distinguish between the two loads. FollowingEquation 2.7 shows how to calculate the transient energy function of an aggregated load.Vk = vk − vk−1 Ik = ik + ik−12 Ut =U1phtransient =K∑k=0Vk Ik (2.7)in which VkandVk−1 are derived from the transient voltage for Kth and (k+1)th samples.I(k) is the average of transient current received from smart meter measurement. Thisstudy, compares the results of the Short Time Fourier Transform (STFT) and DiscreteWavelet Transform (DWT) of the turn on transients of three different loads. Figure 2.10shows the result of applying short term Fourier transfer on transient current of three dif-ferent categories of loads. The transient response time of each load can be identified fromthe ending time and the starting time based on Figure 2.10.25Figure 2.10: The STFT results on turn on current transients. (a) 160 horse powerinduction motor, (b) 123 horse power induction motor driven by variable volt-age driver. (c) A bank of loads supplied by a six-pulse thyristor rectifier for ACpower [37].In order to compare STFT method with DWT method, Wavelet transform of same threewaveforms are captured and pictured as following: DWT is more robust since it shouldnot use fixed window frame same as STFT. Size of the window is scalable in DWT whichenables this method to capture all information from starting point to ending point (samestarting time and end time as load’s signal).In Figures 2.11, 2.12 and 2.13, first picture is instantaneous power waveforms of turn-on transients for a 160 hp induction motor, a 123 hp induction motor driven by variable-voltage drivers and a bank of loads supplied by a six-pulse thyristor respectively and2nd plot demonstrates discrete wavelet coefficients which are distinguishable for different26Figure 2.11: The DWT results of turn on instantaneous power transients for a 160 hpinduction motor [37].appliances.High frequencies harmonics envelopes are the other strong descriptor for detecting aspecific type of load in particular, power electronic types. Many loads, inherently drawnon-sinusoidal, distorted current due to their physical characteristics. Taking Fouriertransform of their current waveforms, aims to compute “spectral envelopes” which uniquelycharacterises the specific type of the load. In [38] a computer and an incandescent lightbulb got detected in the ∆P, ∆Q and ∆3rd harmonic, coordinate system as shown in Figure2.14.In [25] matrix pencil method is used to represent the current in time value in termsof conjugate poles and residues. 1st, 3rd and 5th current harmonics are studied for over 9different load classes including, incandescent Lamp, halogen Lamp, economy Lamp, wa-ter heater, electric convector, oven, hot Plate (one and two burners), television, Computer)and PC. There is a data set of 900 samples prepared to be passed in to a selected classifier.The drawn current by each appliance is the sum of M complex valued sinusoid signalsweighted by complex residues expressed in Equation 2.8.it =M∑m1rm.exp(am + j2pi. fm)t+ b(t) (2.8)27Figure 2.12: The DWT results of turn on instantaneous power transients for a 123 hpinduction motor driven by variable voltage driver [37].The discrete current signal is expressed in Equation 2.9:ik =M∑m1rm.zkm + b(k) k = 1,2, · · · ,N (2.9)Wherezm = exp(am + j2pi. fm)ts m = 1,2, · · · ,M (2.10)In Equation 4.1, am is the attenuation factor, fm is frequency, ts is the sampling timeand rm is the residue of mth component.The other innovative approach to load disaggregation arises from harmonic energy coef-ficients. This study, carried out different measurements to deduce that steady-state har-monics contents, are more repeatable than transients harmonics contents. Therefore, theyconcluded that steady state harmonics are more reflective of the loads behaviour. For thisreason, just before and after an event, steady state harmonics contents, should be exam-ined. This method measures currents harmonics content of the electrical loads, from athree-phase environment. A vector of length 27 corresponds to the first 8 harmonics plus28Figure 2.13: The DWT results of turn on instantaneous power transients for a bank ofloads supplied by a six-pulse thyristor rectifier for AC power [37].the fundamental frequency content, for all three phases is developed as the load’s fea-ture’s database. An observed new load is identified with comparing its harmonics energycontent to the reported library’s harmonics vectors and minimizing the match error [24].In [22] researchers assert that a load’s current waveform polluted with harmonics can berepresented as a normalized energy vector consisting of five elements. They took sam-ples of three loads ( personal computer, fluorescent lamp and dimmer for incandescent)current waveforms at the rate of 10 KHz with total period time of 51.2 ms. then, 5 leveldecomposition of wavelets has been applied for analysis. It is shown that for each type ofload the wavelet coefficients curves give distinctive signatures and they used these signa-tures to identify each specific type of load. The other innovative approach, considers thesteady-state starting noise of each appliance as a distinctive load feature. approach usessingle plug in sensor, to detect both the abrupt noise created with an abrupt switchingoperation, and also the noise created by the certain devices in operation. The single plug-in model connects to a PC and PC records the generated noise caused by an appliance’sswitching ON or OFF as shown in Figure 2.15. Learning classifier, learns the certain char-acteristics of the noise to detect what appliance is turned ON or OFF in future. Authorsin this paper, categorise the noise into two categories named as transient noise and steady29Figure 2.14: A computer and an incandescent light bulb in the ∆ P, ∆ Q and ∆ 3rdharmonic coordinate system [38].state noise. transient noise pulse is created while a large motor load such as a fan is turnedon while, transient noises last only for microseconds and contain rich spectrum of the fre-quency components which is ranged from 100 Hz-10 KHz.This study assumes noise signature of a particular device depends on both type of thedevice and the transmission line behaviour of the interconnecting power line. They cap-ture both contributors in the load identification process. Three main classes of electricalresidential loads, are considered in this research.1. Pure resistive loads such as light bulbs and stoves, do not create detectable amountof electrical noise.2. Motor loads like blender and fan which produce transient voltage noise synchronousto 60 Hz.3. Power electronic appliances, such as microwave, PC and chargers which emit the30Figure 2.15: Frequency spectrum of a particular light switch being toggled (on andoff events). Graphs indicate amplitude at each frequency level. Events at (b) werecaptured two days after (a). (c) events are captured one week after (a). On and offevents are different enough to be distinguished.steady-state noise synchronous to the internal oscillator.This study was able to classify the various electrical events using SVM learning algorithmwith accuracies ranging from 85-90% [23].2.4 Chapter SummaryThis chapter, summarized previous research and studies in the area of load disaggrega-tion. Each method has its own advantages and disadvantages as explained. The mutuallimitation in all above algorithms, is that in most of them, we need to have a large accu-rate load library, including different types of loads in order to train a classifier which willdetect a new type of load based on seen signatures in future. Also, none of the methods,consider disaggregating the total load from the distribution feeder level, which is moreefficient from an electrical grid point of view.31Chapter 3TheoryAs discussed in Chapter 2, researchers did not come up with an efficient algorithm to de-tect all power grid connected loads. Most of the studies focus on load identification basedon a supervised learning method, meaning that a learning system observes training dataset’s behaviour, thereafter it recognizes a new observed load type. This sort of approach isrestricted from the fact that if a new observed load is not operationally similar to the oneswhich system already is trained with, it will not detect the correct load’s type. The othershortcoming arises from the fact that non of recent researches look at load disaggregationproblem from distribution feeder point of view. In other words, the ability to detect totalload’s components connected to a specific feeder at each time instant. As yet, people whoare working in this domain, are looking at each customer’s entry in particular residentialloads.This thesis proves that looking at aggregated load from distribution feeder instead ofeach customer’s entry point will benefit both the power system engineers and also thecustomers. It helps Distribution and Transmission engineers to capture a larger picture ofconnected load’s structure. Because all consumer’s load information is gathered togetherand each feeder’s data is sent back to control centre for analysis.The thesis introduces an algorithm to obtain feeder level load structures and electrical pa-rameters using EMTP discrete time based solutions. This is the first time that there is aconnection between Load Disaggregation problem and EMTP analysis solution.We developed a method which enables engineers at control centre to know exactly howmany percent of the aggregated load is motor type, purely resistive or purely inductive.We are able to report the load type for all feeders. This will help in a sense that most32powerful tools in power systems need to know the exact type of load at a large picture.Estate estimation, power flow studies, stability related studies, etc. all need to know totalload’s major components.This thesis employs EMTP technique that helps predict the order of total load. Addi-tionally, the thesis adopts network synthesis methods for the sake of recognizing exactparameters of aggregated load constitutes.In Section 3.3 different electromagnetic transient program integration (discretization) rulesare described and compared. In Section 3.5 Foster and Cauer methods are introduced andconsidered to construct load’s electrical circuit based on Z (aggregated load impedance)and Laplace transform function formulated from EMTP coefficients.But before starting the theory, I give a brief overview of the two power analysis tools,EMTP and PSCAD in Sections 3.1 and 3.2 respectively, because these two power tools arethe foundation of my thesis.3.1 EMTP SoftwareThis section gives an overview of the EMTP solution with the purpose to provide the suf-ficient background to reach a full understanding of the proposed method.EMTP (Electro Magnetic Transient Program) is a general purpose computer program forsimulating transient events in power system. The program features wide variety of mod-elling capabilities. The beauty and uniqueness of EMTP solution are rooted in the natureof the discretization process. It focuses on discrete domain and offers solution methodsencompassing Trapezoidal, Backward Euler, Forward Euler, etc. Herman Dommel devel-oped EMTP in late 1960s at Bonneville Power Administration (BPA) which considered theprogram as a digital computer [39].EMTP solution, first, defines discretization of each elements of the network instead of dis-cretizing the circuit’s full sets of differential equations. In other word, in the EMTP everycomponent’s relationship between voltage and current is modelled with a discreet timeequivalents circuits consisting of resistance and source combinations. The values dependon the integration method and sampling ratio.This thesis assumes every single load a combination of passive elements R, L and C, whichare connected together to make a load’s circuit. Therefore, we behave a load like an elec-trical circuit consisting of R, L and C components which can be solved via EMTP discrete-time solutions. The initiative idea behind our solution for load disaggregation problem33is in form of EMTP-like solution techniques. EMTP can solve networks consisting of in-terconnections of resistances, inductances, capacitances, distributed- parameter lines andother certain elements. Another reason motivated employing EMTP as the main algo-rithm is small discretizing time step used in simulations. Small discretzation helps, tracethe instantaneous value of our signals in real time. Smart meters are based on the sametechnique, i.e., they send out real time voltage and current values to the control centre.Popularity of smart meters enthused this thesis toward employing EMTP technique.3.2 PSCADSystem’s behaviour, can be studied in time domain or frequency domain. EMTDC (elec-tromagnetic transients including DC) is the most powerful electro-magnetic transientssimulation engine for time domain simulations [40].This engine started becoming a power tool in 1975, at Manitoba Hydro by Dennis Wood-ford (Executive Director of the Centre 1986 - 2001), which at the time was sufficientlypowerful and flexible to study the Nelson River HVDC power system.PSCAD software is the advanced graphical user interface (GUI) of EMTDC which en-ables user to model the circuit, run it and analyse data extracted from the software inother types of environment such as Matlab. EMTDC serves time domain electromagnetictransient solution engines for the family of PSCAD software. PSCAD is widely used toanalyse, model and study power system including DC systems, AC ones and also powerelectronics.PSCAD is the main simulator in our research since it is accurate, user friendly, and com-patible with Matlab. More importantly it is written based on EMTP solution techniques.The main idea appears in 1969 by Dr. Hermann Dommel published at [39]. Electromag-netic transient solution methods considers a fixed or variable time step (in our case isfixed) to solve differential equations representing network behaviour in time domain. Thesolution is the step by step proceeding along the time axis. Each step depends on its previ-ous value in time, i.e., history parameters. In other words, state of the system is measuredduring t = 0, t = ∆t, t = 2∆t, · · · , t = n∆t. At each state we need to examine previous timesteps as well.343.3 Discretization MethodsThere are various methods applied to solve the differential equations using different inte-gration rules. The common benefit in all these methods is the capability of converting dif-ferential equation to a simple algebraic one including voltage, current and some historyvalues. EMTP software converts electrical circuit to the simple equivalent resistive andthen applies integration rules to calculate voltage and injected current. The main purposeis to model lumped passive components (resistance, inductance and capacitance) whichis easily solved in discrete time domain afterwards.Integration rules transform differential equations describing steady state and transientbehaviours of the system’s components to the discrete time models at finite time incre-ments. A complete network is then formulated to solve for the system conditions at eachtime step. Following Figure 3.1 shows the table of different integration methods.Figure 3.1: Different integration methods description.353.3.1 Forward EulerForward Euler is based on calculating the area under polygons between two endpoints ofthe integral. In other words, forward Euler approximated the integral value by productionof value of function at xn i.e., f (xn)times h = xn − xn−1. Figure 3.2 illustrates the methodit more clearly [41].In fact it is the expansion of Taylor series in two terms [41].xn+1=xn+ h/1! f (xn, tn) + h2/2! f ′(xn, tn) + h3/3! f (2)(xn, tn) + ...−→ xn+1=xn+ h ∗ f (xn, tn)Figure 3.2: Approximation of x(t) by Forward-Euler method [41].3.3.2 Backward EulerBackward Euler is so similar to Forward Euler in a sense that, both considers the areaunder polygons embedded by integral endpoints, but in Backward Euler the integralvalue equals the production of h = xn − xn+1 into f (xn+1) as you can see in Figure3.3.xn+1=xn + h ∗ f (xn+1, tn+1)36Figure 3.3: Approximation of x(t) by Backward Euler method [41].A quick review of these two methods shows that forward Euler gives larger value thanthe real one while backward Eulers output is less than real value of function.3.3.3 TrapezoidalHere is where Trapezoidal rule kicks in by taking average of two values f (xn+1) and f (xn).In other words the value of integral between xn and f (xn+1) is the area under trapezoidABCD shown in Figure 3.4.xn+1=xn +f (xn+1,tn+1)+ f (xn,tn)2 ∗ h37Figure 3.4: Approximation of x(t) by Trapezoidal method [41].In fact in each integration step the average value of the intervals beginning and end istaken into account. It is the most accurate 2nd order integration method.3.3.4 Simpson’sThe other integration rule is called Simpsons rule which is Newton-Cotes formula usingsecond order parabolas in estimating integral value. It uses quadratic polynomial insteadof piecewise linear segments as in trapezoidal method. As can be followed in Figure 3.5that the approximation for the integral function is the following:xn+1=xn−1 + h/3 ∗ ( f (xn−1, tn−1 + 4 f (xn, tn) + f (xn+1, tn+1)All these methods instead of giving continues solution of transients, aggregates snapshotsof all discrete intervals with width of ∆t as the solution. These integration methods de-scribed in Sections 3.3.1, 3.3.2 and 3.3.4 cause truncation error leading to numerical insta-bility. Trapezoidal integration rule as described in section 3.3.3, is selected as my adoptedmethod to distinguish load’s circuits behaviour because it is simple, accurate and numer-ically stable [39].38Figure 3.5: Approximation of x(t) by Simpson methodIn this thesis our main focus is predicting load types based on their unique equivalentcircuit. Motor types, lighting types and heating ones are our main three classes. All typesof residential, commercial and industrial loads are categorized under motor, heating orlighting categories. Motor class represents induction motors. Lighting is symbolic ofinductance types of loads. And lastly, heating is illustrative of resistive one. Our firstassumption is that lighting types are modelled by “L”, heating ones are represented by“R” and motor are modelled using equivalent model of an induction motor. Further in-vestigation shows that we can assume load circuits as different combinations of L, R andin lumped form. Then we apply Trapezoidal rule to solve the time domain equation.Trapezoidal Equivalent ModelsIn the following, I first, provide RLC trapezoidal equivalent models analysis to show howthree main components of our study (R, L and C) are modelled in EMTP using Trape-zoidal integration method. Then, we will calculate discrete time coefficients for complexcombinations of RLC elements in the form of an induction motor equivalent circuit andetc.39Inductance ModelFor the case of an inductor, it is modelled as demonstrated in Figure 3.6. The relationbetween voltage and current across is given by following calculations:Figure 3.6: Inductance equivalent EMTP model.V12(t) = L. di12dtV12(t).dt = Ldi12(t)∫ tt−∆tV12(t)dt = L.∫ tt−∆t di12(t)∫ tt−∆tV12(t)dt = L.[i12(t)− i12(t− ∆t)]Where ∆t is the descretization time step.TrapezoidalDiscretization−−−−−−−−−−−−−→ V12(t)+V12(t−∆)2 .∆t = L.[i12(t)− i12(t− ∆t)]40Then the voltage across the inductor is:−→ v12(t) = 2L∆t︸︷︷︸Equivalent Impedancei12(t)− 2L∆t i12(t− ∆t)− v12(t− ∆t)︸ ︷︷ ︸History Terms(eh(t))(3.1)Figure 3.7: Voltage across the inductor [42].Equation 3.1 describes the equivalent impedance value as a function of inductance (L)and time step ∆t. Also, voltage value across the inductance at time t depends on historyterm which includes current value at previous time step and voltage value at t=t-∆tCapacitance ModelFor the case of a capacitor it is modelled as demonstrated in Figure 3.8. Equation 3.2describes the equivalent impedance value as a function of capacitance (C) and time step∆t. Also, voltage value across the capacitance at time t depends on history term whichincludes current value at previous time step and voltage value at t=t-∆t.41Figure 3.8: Capacitance equivalent EMTP model.i12(t) = Cdv12(t)dt∫ tt−∆t i12(t)dt = C.∫ tt−∆t dv12(t)dt∫ tt−∆t i12(t)dt = C.[v12(t)− v12(t− ∆t)]TrapezoidalDiscretization−−−−−−−−−−−−−→ i12(t)+i12(t−∆t)2 .∆t = C.[v12(t)− v12(t− ∆t)]Then the voltage across the capacitance is:−→ v12(t) = ∆t2C︸︷︷︸Equivalent Impedancei12(t)+∆t2Ci12(t− ∆t) + v12(t− ∆t)C︸ ︷︷ ︸History Terms(eh(t))(3.2)42Figure 3.9: Voltage across the capacitance [42].Resistance ModelFor the case of a resistor it is modelled as demonstrated in Figure 3.10. The relation be-tween voltage and current across is described in Equation 3.3.Figure 3.10: Voltage across the resistance [42].For the case of a resistance the relation between voltage and current across it is givenby:V(t) = R.i(t) (3.3)433.4 The Thesis Method DiscussionIn this thesis hereafter, all the current coefficients are annotated by letter ‘a’ and all voltagecoefficients are annotated by letter ‘b’. Starting from i(t), the first coefficient is a1, thenfor i(t − ∆t) the corresponding coefficients is a2 and so on for i(t − n.∆t) which is an.Same for voltage coefficient’s notation, we have b1 for v(t − ∆t), b2 for v(t − 2∆t) andbn corresponding to v(t − n∆t). Substituting each passive element with its trapezoidalmodel leads to the following well known EMTP formulation which relates voltage andcurrent time series values:a0i(t) + a1i(t− ∆t) + a2i(t− 2∆t) + · · ·+ ani(t− n∆t) =b0v(t) + b1v(t− ∆t) + b2v(t− 2∆t) + · · ·+ bnv(t− n∆t)(3.4)For t = 0,∆t,2∆t · · · ,n∆tIn order to solve Equation (3.4), all terms are brought to one side of the equation exceptterm v(t). The new form of equation is as following:v(t) = a′0i(t) + a′1i(t− ∆t) + a′2i(t− 2∆t) + · · ·+ a′ni(t− n∆t)+b′0v(t− ∆t) + b′1v(t− 2∆t) + · · ·+ b′nv(t− (n+ 1)∆t)(3.5)For t = 0,∆t,2∆t, · · · ,n∆tIn order to find a′0, a′1, · · · , a′n,b′0,b′1, · · · ,b′n values, we need instantaneous voltage and cur-rent values at (t = 0,∆t, 2∆t,· · · ,n∆t) as the known parameters. Voltage and current co-efficients will be the unknown variables. Before solving Equation 3.6, we need to deter-mine the order of the load circuit. The load’s order is determined with solving z(t) =v(t)i(t) equation. Trapezoidal discretization rule will be applied on z(t) same as in Sec-tions 3.3.4 and 3.3.4. Integration of the impedance function using trapezoidal rule, wewill derive an equation in the same form as Equation (3.4). The largest ∆t index wherex= 1, · · · ,n, will be the load order. Equation (3.5) is in the form of AX= B or (AX− B= 0)where two main matrices A and B include instantaneous values of voltage and currentat time t, t − ∆t, t − 2∆t, · · · , t − n∆t. The rows include i(t − ∆t) · · · , i(t − n∆t)and alsov(t−∆t) · · · ,v(t−n∆t) terms for one specific sampling scenario while number of columnsis determined by the order of the load and interprets the total number of sampling points.Here is a symbolic representation of A and B matrices:44A=is1(t) is1(t− ∆t) . . . is1(t− n∆t) vs1(t− ∆t) vs1(t− 2∆t) . . . vs1(t− n∆t)is2(t) is2(t− ∆t) . . . is2(t− n∆t) vs1(t− ∆t) vs2(t− 2∆t) . . . vs2(t− n∆t)........................isn(t) isn(t− ∆t) . . . isn(t− n∆t) vs1(t− ∆t) vsn(t− 2∆t) . . . vsn(t− n∆t)Where ”n” is the order of our load. Depending on the order matrix size is defined.s1, s2, . . . , sn defines sample numbers, i.e., first sampling point, second sampling point andso on. Rows correspond to the samples and columns correspond to different time stepsvalues. B=vs1(t)vs2(t)...vsn(t)B matrix has n rows where n is again number of samples (or order of circuit) and it al-ways has one column including voltage values at time t. Inverse (A)*B gives the solutionmatrix for Equation (3.6) which describes the numerical relationship between voltage andcurrent in time.is1(t) is1(t− ∆t) · · · is1(t− n.∆t) vs1(t− ∆t) vs1(t− 2∆t) · · · vs1(t− n.∆t)is2(t) is2(t− ∆t) · · · is2(t− n.∆t) vs2(t− ∆t) vs2(t− 2∆t) · · · vs2(t− n.∆t)..................isn(t) isn(t− ∆t) · · · isn(t− n.∆t) vsn(t− ∆t) vsn(t− 2∆t) · · · vsn(t− n.∆t) ∗a0a1a2...anb0b1...bn=vs1(t)vs2(t)...vsn(t)(3.6)Basically, the result of Inv(A) ∗ B is a matrix with (2n+ 1) rows and one column wherethe first (n+ 1) coefficients correspond to current characteristics and the rest n coefficientsare numerical descriptor of voltage terms. Note that order of the circuit takes a decisionof number of sampling points required. For example if we have a circuit with order=n,numbers of sampling points required to calculate the matrices A and B is (2n+ 1).453.4.1 Load’s Electrical Circuit Order Determination Knowing the Structure ofthe LoadHere is a general sequence summarizing how the order of an electrical load is obtainedknowing the structure of the load’s circuit:1. Equivalent Z (impedance) of load’s circuit is calculated in the form of Laplace do-main (S domain). (Writing down the impedance of load in the form Z(s) = V(s)I(s) ,then expanding the formula, we get V(s) = Z(s).I(s).Resulted S domain expanded equation is converted to a discrete time domain ex-panded equation by substituting S parameter with d/dt, s2 with d/dt2 up to re-placement of sn with d/dtn.2. Integrating on both sides of the main equation and applying Trapezoidal discretiza-tion, results in a same format equation same as Equation (3.4).3. the maximum coefficient of ∆t appeared among history terms, determines the orderof the circuit.Up to this point, we discussed how to determine a load network’s order with having itsstructure. Now it is time to determine the order of load’s circuit with having no infor-mation about the load’s topology. Because the main aim of this research is to develop analgorithm detecting feeder connected load’s structures.3.4.2 Load’s Electrical Circuit Order Determination Unknowing the Structureof the LoadWe summarize the process of predicting an electrical circuit’s order having no knowledgeabout the topology of the circuit.There are two distinct features about resultant discrete time coefficient which makes theman excellent candidates to be chosen as for the load distinctive signatures.1. Since we assume an aggregated load as one single electrical circuit at each time step,for different harmonics scenarios these resultant matrix from Inv(A) ∗ B are equal.In other words, these time series coefficients are only function of circuit parameters(load passive elements) which are constant for one specific type of load, and also∆t which is the same for different harmonics as far as we use the same samplingfrequency being able to capture major signal information of that circuit.462. For different sampling scenarios (of course which samples to take is important)as discussed previous statement, corresponding coefficients are function of constantcircuit passive elements and ∆t.Again, choosing the same sampling frequency with selection of same initial samplepoints with different sampling windows, or different initial samples but with thesame window length we get the same values as for these coefficients.In this thesis, both approaches are evaluated on different combinations of loads in PSCADsimulation which lead to the same results. The first approach is more simple and practicalfor types of loads which produce harmonics such as power electronics types. Extractingat least two different frequency content both from voltage and current signals and calcu-lating discrete time coefficients for each set, will give the order of the connected load. Butlooking an ideal power system world, goal is to reduce the harmonics effects, nowadays.Accordingly, second approach is considered as the main validation method for the restof this study. In favour of predicting a connected load’s order based on our approach,an initial order is considered, starting from order=0 (it happens for resistive loads). Atleast two distinguishing groups of coefficients derived from different sampling scenariosas demonstrated in Figures 3.11 and 3.12 are taken into account. If the obtained discretetime coefficient values for the selected scenarios, were equal, it means that guessed orderwas correct. Otherwise, order is increased by one and the same process is followed forthe new order of the circuit. The process is repeated until for an order, the output ma-trix (Inv(A) ∗ B) shows unequal values for different sets of sampling points. We startedvalidating this idea based on simulations at PSCAD. In order to create different test casescenarios being able to compare different sets of coefficients, three sampling schemes areintroduced:1. Various width windows having same start point. (Sliding windows along timeaxis). We achieve this situation by changing the distance between sampling pointsas shown in Figure 3.11.47Figure 3.11: Different length windows having the same initial point as one of “volt-age/current sampling points” selection method.2. Windows with different starting point while having the same length as shown inFigure 3.12.Figure 3.12: Same length windows having different initial point as one of “volt-age/current sampling points” selection method.483. Analysing two 1st and 2nd voltage and current harmonics as an input, calculatingthe corresponding Matrix for both sets which is not employed in this thesis.Flowchart in Figure 3.13 exhibits this research’s algorithm. The first step is receiving timeseries voltage and current data. The next step is prediction of load’s order and finally,network synthesis methods help calculating load’s parameters. In the following section,START Voltage and Current signals are derived from smart meters.  Start from order=n for (n=0… Inf):     n=n+1  Start from order=n for (n=0… Inf):      n=n+1 Selecting at least two different sets of time sampling windows (number of samples=2(order of circuit) +1) Selecting at least two different sets of harmonics voltage and current time series signals (number of samples=2(order of circuit (n)) +1) Calculate impedance coefficients in the form of equation 1 Impedance coefficients are equal for two different sampling sets Form Voltage and Current as following format:    Convert time series equation to S domain Foster/Cauer network synthesis method are applied to calculate exact parameters (R,L and C) for the load circuit  No Yes Figure 3.13: Proposed method: Process from smart-meter data capturing to loadidentification.examples of parametric discrete time coefficients, describing the relation between voltageand current are calculated for few major load types.493.4.3 Examples of Parametric Equivalent RLC Circuits Trapezoidal ModelsCase One: Inductive Load-First Order ExampleThis example is already analysed in equation 3.1. Maximum ∆t appeared in the expandedequation form is one, which confirms that a purely inductive load is a first order systemwith following coefficients: a0 = 2L∆ta1 = − 2L∆tb0 = −1 This example calculates time series voltage and current coefficients for an ”In-ductor” . This category of load represents lighting eigen-loads.Case Two: RL Parallel Load-First Order ExampleFigure 3.14: Parallel RL case.V(s)I(s) =RL(s)R+LS −→RV(S) + LSV(s) = RLsI(s) convert−to−time−domain−−−−−−−−−−−−−→RV(t) + L(dV/dt) = RL.(di/dt) −→∫ tt−∆tRV(t) + L(dV/dt)dt =∫ tt−∆tRL.(di/dt)dtTrapezoidal−integration−−−−−−−−−−−−→ R[V(t)+V(t−∆t)2 .∆t] + L[V(t)−V(t− ∆t)] = RL[i(t)− i(t− ∆t)] −→V(t)[R.∆t2 + L] + [R.∆t2 − L].V(t− ∆t) = RLi(t)− RLi(t− ∆t)50−→ V(t) = [ 1R.∆t2 +L].[RLi(t)− RL i(t− ∆t) ]− [ R.∆t2 −L][ 1R.∆t2 +L]. V(t− ∆t) ] −→Voltage and current coefficeints are :a0 = RLR.∆t2 +La1 = RLR.∆t2 −Lb0 =[ R.∆t2 −L]R.∆t2 +LIn this case highlighted terms in Red are the largest ∆t history terms. Based on our algo-rithm, first order is determined as the correct order for his type of load. For a resistanceload, no ∆t term will appear since impedance is a constant value of R, no history termsexist in this situation. Therefore, for all resistive loads, order=0 is assumed. Adding aresistive load in parallel or series with other loads will not change the original order ofcircuit.Case Three: LC Parallel With A Resistance Load-Second Order ExampleFigure 3.15: A parallel LC in series with a resistance case.51V(s)I(s) =RLCS2+LS+RLCS2+1 −→LCs2.V(s) +V(s) = RLCs2.I(s) + LSI(s) + RI(s) convert−to−time−domain−−−−−−−−−−−−−→V(t) + LC.d2V/dt2 = RLC.d2i/dt2 + L.di/dt+ Ri(t) −→∫ tt−∆t[V(t) + LC.d2V/dt2]dt =∫ tt−∆t[RLC.d2i/dt2 + L.di/dt+ Ri(t)]dttrapezoidal−integration−−−−−−−−−−−−→LC[dV/dt− dV(t−∆t)/dt] + [V(t)+V(t−∆t)2 .∆t] = RLC[di/dt− di(t−∆t)/dt] + L[i(t)− i(t−∆t)] + R.[ i(t)+i(t−∆t)2 .∆t]trapezoidal−integration−−−−−−−−−−−−→ RLC.[i(t)− i(t− ∆t)]− RLC.[i(t− ∆t)− i(t− 2∆t)] + (L+ R.∆t2 ).[ i(t)+i(t−∆t)2 .∆t]+ (R.∆t2 − L).[ i(t−∆t)+i(t−2∆t)2 .∆t] = LC.[V(t)−V(t− ∆t)]− LC.[V(t− ∆t)−V(t− 2∆t)] + ∆t2 .[(V(t)+V(t−∆t)2 .∆t)+ (V(t−∆t)+V(t−2∆t)2 .∆t)] −→i(t)[RLC+ (L+ R.∆t2 ).∆t2 ] + i(t− ∆t).[−2RLC+ (L+ R.∆t2 ).∆t2 + (−L+ R.∆t2 ).∆t2 ]+ i(t− 2∆t).[RLC+ (−L+ R.∆t2 ).∆t2 ] = v(t).[LC+ ∆t24 ] + v(t− ∆t).[−2LC+ ∆t22 ] + v(t− 2∆t).[LC+ ∆t24 ] −→i(t).[R.∆t24 + L∆t2 + RLC] + i(t − ∆t).[R.∆t22 − 2RLC] + i(t− 2∆t) .[R.∆t24 − L∆t2 + RLC] = v(t).[LC + ∆t24 ] + v(t −∆t).[−2LC+ ∆t22 ] + v(t− 2∆t) .[LC+ ∆t24 ] −→ Voltage and current coefficeints are :a0 =[R. ∆t24 +L∆t2 +RLC][LC+ ∆t24 ]=R.[ ∆t24 +LC]+L∆t2[LC+ ∆t24 ]= R+ L∆t2[LC+ ∆t24 ]a1 =[R. ∆t22 −2RLC][LC+ ∆t24 ]a2 =[R. ∆t24 −L ∆t2 +RLC][LC+ ∆t24 ]=R.[ ∆t24 +LC]−L ∆t2[LC+ ∆t24 ]= R+ −L∆t2[LC+ ∆t24 ]b0 = − [−2LC+∆t22 ][LC+ ∆t24 ]b1 = − [LC+∆t24 ][LC+ ∆t24 ]= −152Highlighted terms in Red show 2∆t terms which determines this circuit as being the2nd order system.Case Four: Motor Load-Second Order ExampleThis example calculates time series Voltage and Current coefficients for an “InductionMotor” . Figure 3.16 represents all motor types load which is one of 3 major eigen-loadsin this study. Highlighted history terms have the largest index for ∆t which confirmsmotor loads as being 2nd order systems.Figure 3.16: Induction motor type.53V(s)I(s) = [(R2 + L2 ‖ LmS] + [R1 + L1S] = R1.R2+S.[R2.L1+R1.L2+R1.Lm+R2.Lm]+S2.[L1.L2+L1.Lm+L2.Lm]R2+S.[Lm+L2]−→ R2.V(s)+ [Lm+ L2].SV(S) = R1.R2.I(s)+ [R2.L1 +R1.L2 +R1.Lm+R2.Lm].SI(s)+ [L1.L2 + L1.Lm+ L2.Lm].S2 I(S)convert− f romS−domain−to−t−domain−−−−−−−−−−−−−−−−−−−−−→ R2.v(t) + [Lm + L2]. dv(t)dt = R1.R2.i(t) + [R2.L1 + R1.L2 + R1.Lm + R2.Lm]. di(t)dt+ [L1.L2 + L1.Lm + L2.Lm].d2i(t)d2t −→∫ t−∆tt [R2.v(t) + [Lm + L2].dv(t)dt ]dt=∫ t−∆tt [R1.R2.i(t) + [R2.L1 + R1.L2 + R1.Lm + R2.Lm].di(t)dt + [L1.L2 + L1.Lm + L2.Lm].d2i(t)d2t ]dtTrapezoidal−integration−−−−−−−−−−−−→ R2.( v(t)+v(t−∆t)2 .∆t) + [Lm + L2].[v(t)− v(t− ∆t)]= R1.R2.( i(t)+i(t−∆t)2 .∆t) + [R2L1+ R1.L2+ R1Lm+ R2Lm].[i(t)− i(t− ∆t)] + [L1.L2+ L1.Lm+ L2.Lm].[ di(t)d(t)− di(t−∆t)dt ] −→ [R2.∆t2 + [Lm+ L2]].V(t) + [R2.∆t2 − [Lm+ L2]]V(t− ∆t)= [R1.R2.∆t2 + R2.L1+ R1.L2+ R1.Lm+ R2.Lm].i(t) + [R1.R2.∆t2 − [R2.L1+ R1.L2+ R1.Lm+ R2.Lm]].i(t− ∆t) + [L1.L2 + L1.Lm + L2.Lm].[ di(t)dt − di(t−∆t)dt ]Trapezoidal−integration−−−−−−−−−−−−→ [R2.∆t2 + [Lm+ L2]].[V(t)+V(t−∆t)2 .∆t]+ [R2.∆t2 − [Lm+ L2]].[V(t−∆t)+V(t−2∆t)2 .∆t] = [R1.R2.∆t2 + R2.L1+ R1.L2+ R1.Lm+ R2.Lm].[ i(t)+i(t−∆t)2 .∆t]+ [R1.R2.∆t2 −]R2.L1+ R1.L2+ R1.Lm+ R2.Lm]].[ i(t−∆t)+i(t−2∆t)2 .∆t]+ [L1.L2 + L1.Lm+ L2.Lm].[i(t)− i(t−∆t)]− [L1.L2 + L1.Lm+ L2.Lm].[i(t−∆t)− i(t− 2∆t)]−→ [R2.∆t2 +[L2+Lm]].∆t2 .v(t)+ R2.∆t22 .v(t−∆t)+[R2. ∆t2 −[L2+Lm]].∆t2 . v(t− 2∆t) = [R1.R2. ∆t2 +[R2.L1+R1.L2+R1.Lm+R2.Lm].∆t2 +[L1.L2+ L1.Lm+ L2.Lm]].i(t)−2[L1.L2 + L1.Lm + L2.Lm].i(t− ∆t) + [R1.R2.∆t2 −[R2.L1+R1.L2+R1.Lm+R2.Lm].∆t2+ [L1.L2 + L1.Lm + L2.Lm]]. i(t− 2∆t) −→ a0 = [R1.R2. ∆t2 +[R2.L1+R1.L2+R1.Lm+R2.Lm].∆t2 +[L1.L2+L1.Lm+L2.Lm]][R2. ∆t2 +[L2+Lm]].∆t2a1 =−2[L1.L2+L1.Lm+L2.Lm][R2. ∆t2 +[L2+Lm]].∆t2a2 =[R1.R2. ∆t2 −[R2.L1+R1.L2+R1.Lm+R2.Lm].∆t2 +[L1.L2+L1.Lm+L2.Lm]][R2. ∆t2 +[L2+Lm]].∆t2b0 = −R2.∆t22[R2. ∆t2 +[L2+Lm]].∆t2b1 = −[R2. ∆t2 −[L2+Lm]].∆t2[R2. ∆t2 +[L2+Lm]].∆t254Discussion Based on Parametric CalculationsAs earlier discussed in this chapter, our new developed model recognizes three majortypes of loads as eigen-loads which includes induction motor, inductive and resistive.The reason of choosing these three categories of loads is that each of these types havetheir own voltage dependency characteristic. UBC Power Systems Lab invented a newLPF (Linear Power Flow) formulation instead of non-linear power flow solution. In theirmethod, loads are modelled considering their voltage dependence characteristics whichallows for a load representation with a constant-impedance and a constant-current syn-thesis (Z-I model) [11]. Above mentioned voltage dependence is determined with know-ing type of feeder connected loads. This fact triggered this thesis to investigate about loaddisaggregation more deeply and come up with an efficient pattern. As a result, the algo-rithm is able to detect motor, inductive and resistive contribution in the total load. Basedon few example calculations performed for a motor, an inductor and a resistor, we figuredout that “Resistive loads” do not change the circuit. “Inductive loads” are increasing thetotal order by one. Motor loads are increasing the total load’s order by two. The onlylimitation of this method is that, for the case of an inductor and a resistor, we cannot sayhow many inductors or resistances are in parallel to result in the achieved aggregatedload. Following Table 3.1 is a summary of different category of loads with correspondingorders:Table 3.1: Different loads order prediction.Type of Load Order of Load)Resistance 0Inductance 1Motor 2Resistance parallel with Inductance 1Resistance parallel with Motor 2Inductance parallel with Motor 3n parallel Resistance 0n parallel inductance 1n parallel Motor 2nSo far, we have obtained the order of load circuit. Now it is time to calculate the55electrical parameters of the load. For this reason, we select network synthesis to realizethe load network having its equivalent impedance transfer function.3.5 Network SynthesisIn this section network synthesis is explained. Two general methods namely, Foster andCauer are discussed for the synthesis of networks having two kinds of elements (LC, RC,or RL).Network synthesis is the opposite way of network analysis. In any network synthesisproblem engineers deal with three quantities namely, input, output and network.Note that, we use the calculated discrete time impedance coefficients of a load, in theform of the impedance “Z” domain transfer function. The next step is to covert the dis-crete time transfer function to the continuous time transfer function. Applying networksynthesis techniques to impedance Laplacian transfer function, outcomes electrical pa-rameters of the load.The question is, what is the network topology for a given impedance or admittance trans-fer function.3.5.1 Network Synthesis DescriptionNetwork synthesis is a block which is fed with a driving point such as impedance oradmittance as the input and the output is the equivalent network composed of passiveelements. In other words, it simplifies the interpretation of a closed loop driving pointimpedance or admittance.There are several techniques, which can be found in electrical engineering literature forthe synthesis of two terminal network from the driving point ( impedance/admittancefunction) [43] [44] [45] [46]. We employ the main two common network synthesis meth-ods, Foster forms and Cauer forms in this thesis [47] [48] [49].The oldest and most widely used synthesis forms, are implemented on LC, RC and RLnetworks. LC networks topologies can be generalized to the RL and RC forms as well.The first method introduced in Section 3.6 is Foster method which is based on partialfraction expansion of the impedance or admittance transfer function [50]. in the follow-ing section, the synthesis of LC one port networks (i.e., LC networks having a prescribeddriving-point impedance/admittance) will be described.563.6 Foster FormsFirst Foster FormFirst foster method realizes ZLC in the form of circuit in Figure 3.17.Realization of Zs in the form of Figure 3.17 was first introduced by Foster. As a conse-quence, this form of realization is named the First Foster form or F1 form.Figure 3.17: First Foster form-partial fraction expansion in the impedance form. Thecomponents are connected in series [51].The equivalent impedance which is seen from the input port of the circuit in Fig-ure 3.17 can be formulated as :Z(s) = a1s+a2s+j−1∑i=32.aiss2 + w2i+ zk(s) (3.7)Equation 3.7 can be reformulated in a general form as Equation 3.8.Z(s) =H.(s2 + w12)(s2 + w32) · · · (s2 + wm2)(s2 + w22)(s2 + w42) · · · (s2 + wr2) (3.8)Figure 3.17 is the result of applying 1st Foster technique on impedance Equation 3.8.Second Foster FormSecond foster method realizes YLC in the form of circuit in Figure 3.18.57Figure 3.18: Second Foster form-partial fraction expansion in the admittance form.The components are connected in parallel [51].The equivalent admittance which is seen from the input port of the circuit in Fig-ure 3.18 can be formulated as :Y(s) = a1s+a2s+j−1∑i=32.aiss2 + w2i+Yj(s) (3.9)3.7 Cauer FormsCauer methods apply Continued Fraction on impedance or admittance transfer functionsto realize the LC circuits. First Cauer and Second Cauer topologies are shown in Figures3.20 and 3.21 respectively.Cauer network configurations, which are in the form of the ladder circuits, were first in-troduced by Cauer and named as Cauer forms.The expansion of ZLC in the first Cauer form is accomplished through the alternate use ofF1 (First Foster) and F2 (second Foster) forms. Cauer configurations are obtained throughexpansion of ZLC in a continued fraction by the process of continued division.There are continued Fraction expansion about infinity and continued fraction expansionabout zero which form these two Cauer forms.Continued Fraction expansion about infinity process follows the following steps:1. Z(s) is written in a form that there is a pole at infinity. This pole is removed by longdivision leaving Z1(s).582. Z1(s) is inverted, yielding Y1(s) which has a pole at infinity.3. Each subsequent cycle consists of an inversion and the removal of a pole at s-infinity.4. The resultant is a circuit with capacities in series and inductors in parallel whichdescribes the 2nd Cauer form.Continued fraction expansion about zero is similar to the continued Fraction expansionabout infinity except that at each subsequent it consists of an inversion and the removal ofa pole at zero. The resulting network has the same structure but with capacitors in seriesand the inductors in parallel which is describes the 1st Cauer form.To get a better insight about continued fraction process toward Cauer forms synthesis,consider Figure 3.19 which has the series arms as the impedances and the parallel arms asthe admittances.Figure 3.19: Ladder network topology [52].To compute driving point of equivalent impedance for the circuit in Figure 3.19, wefollow the following steps:1. The impedance of the last term Z6 = 1Y62. The impedance of the last two arms Z5+ 1Y63. The impedance of the last three arms 1y4+ 1Z5+ 1y6Proceeding in this manner results: Z = z1 + 1y2+ 1z3+ 1y4+ 1z5+ 1y6Realization of LC driving point functions as in first and second Cauer forms are demon-strated in Figures 3.20 and 3.21.59Figure 3.20: First Cauer form continued fraction expansion about infinity. When ap-plied to the impedance of a positive LC network which results in a chain of seriesinductors and parallel capacitors [51].Z(s) = L1.s+1C2.s+ 1L3.s+ 1C4.s+···(3.10)Figure 3.20 shows the obtained circuit network resulted from 1st Cauer synthesis methodbeing applied on Equation 3.10.Figure 3.21 shows the obtained circuit network resulted from 2nd Cauer synthesis methodbeing applied on Equation 3.11.Figure 3.21: Second Cauer form continued fraction expansion about zero. When ap-plied to the impedance of a positive LC network which results in a chain of seriescapacitors and parallel inductors [51].60Z(s) =1C1.s+11L2.s+ 11C3.s+ 1L4.s+···(3.11)Figure 3.22 summarizes Foster and Cauer synthesis methods explained so far:Figure 3.22: First Foster, 2nd Foster, first Cauer and 2nd Cauer forms method sum-mary.Table 3.2 gives the general forms of the expansions for all 4 types of Foster and Cauermethods.Table 3.2: Foster and Cauer impedance/admittance expansionspartial fraction of impedance function Z(s) = K0 + K1s−s1 +K3s−s3 + · · ·+ Kns−snpartial fraction of admittance function Y(s) = k∞s+ k0 + K2s−s2 +K4s−s4 + · · ·+ Kms−smcontinued fraction expansion about a point of infinity Z = a1 +1b2s+ 1a3+ 1b4s+ 1a5+···continued fraction expansion about the origin Z =1a1s+1b2 + 1a3s+ 1b4+ 1a5s+···Up to here, we explained how to derive a LC circuit employing Foster and Cauer meth-ods. In Section 3.8 we will discuss about applying network synthesis methods toward RLand RC circuits derivation.613.8 RL-RC Network SynthesisThe synthesis of RC and RL networks is essentially identical. With the change of variables,RC and RL topologies will be converted to LC topology and vice versa. Either Foster orCauer techniques can be applied to RL and RC networks.3.8.1 Realization of ZRC in Foster 1 Form:Figure 3.23: Realization of ZRC in F1 form [53].Evidently, to obtain the F1 form of an RC network, it is necessary to expand ZRC in thepartial fraction form as:ZRC(s) = A0s +A2Ps+w22+ A4Ps+w24+ · · ·+ A2n−2Ps+w22n−2 + A2nAll residues A2k(k = 0,1, · · · ,n) and all poles [w22k](k = 1,2, · · · ,n − 1) are positive andall poles are simple. Hence, all residues of ZRC are positive and all its poles occur ats = −w22k(k = 0,1, · · · ,n) where s = −w22k is real and negative. Therefore all poles of ZRCare simple, lie on the negative real axis ( the -σaxis) and have positive residues.Realization of ZRC in F1 form is summarized as following:1. Expand ZRC into partial fractions.2. Synthesis ZRC as a series connection of impedance terms in the partial fraction form3.8.2 Realization of ZRC in Foster 2 Form:to obtain the F2 form of a RC network, it is necessary to expand YRC in the partial fractionform as:62YRC(s) = A1 + A3ss+w23+ A5ss+w25+ A7ss+w27+ · · ·+ A2n−1ss+w22n−1 + A2n+1sAll A’s and all w’s are positive so that all poles of YRC are simple and lie on -σaxis.Expansion of YRC leads to the circuit configuration in Figure 3.24.Realization of ZRC in F2 form is summarized as following:1. Expand 1s .YRC into partial fractions.2. Form YRC by multiplying each term of the partial fraction expansion by s.3. Synthesis ZRC as a parallel connection of admittance terms in the partial fractionformFigure 3.24: Realization of ZRC in F2 form [53].3.8.3 Realization of ZRC in Cauer 1 Form:Expansion of ZRC leads to the circuit configuration in Figure 3.25.Realization of ZRC in C1 form is summarized as following:1. Arrange the numerator and denominator polynomials in descending order.2. perform continued division3. Assign the quotient of the divisions as branch admittance in the C1 ladder networkFigure 3.25: Realization of ZRC in C1 form [53].633.8.4 Realization of ZRC in Cauer 2 Form:Expansion of ZRC leads to the circuit configuration in Figure 3.26.Realization of ZRC in C2 form is summarized as following:1. Arrange the numerator and denominator polynomials in ascending order.2. perform continued division3. Assign the quotient of the divisions as branch admittance in the C2 ladder networkFigure 3.26: Realization of ZRC in C2 form [53].3.8.5 Realization of ZRL in Foster 1 Form:Expansion of ZRl leads to the circuit configuration in Figure 3.27.Realization of ZRL in F1 form is summarized as following:1. Expand 1s .ZRL into partial fractions.2. Form ZRL by multiplying each term of the partial fraction expansion by s.3. Synthesis ZRL as a series connection of impedance terms in the partial fraction formFigure 3.27: Realization of ZRL in F1 form [53].643.8.6 Realization of ZRL in Foster 2 Form:Expansion of ZRL leads to the circuit configuration in Figure 3.28.Realization of ZRL in F2 form is summarized as following:1. Expand YRL into partial fractions.2. Synthesis ZRL as a parallel connection of admittance terms.Figure 3.28: Realization of ZRL in F2 form [53].3.8.7 Realization of ZRL in Cauer 1 Form:Expansion of ZRL leads to the circuit configuration in Figure 3.29.Realization of ZRL in C1 form is summarized as following:1. Expand YRL into partial fractions.2. Synthesis ZRL as a parallel connection of admittance terms.1. Arrange the numerator and denominator polynomials in descending order.2. perform continued division.3. Assign the quotient of the divisions as branch admittance in the C1 network.Figure 3.29: Realization of ZRL in C1 form [53].651. Arrange the numerator and denominator polynomials in ascending order.2. perform continued division.3. Assign the reciprocals of the quotients of the divisions as branch admittance in theC2 network.3.8.8 Realization of ZRL in Cauer 2 Form:Expansion of ZRL leads to the circuit configuration in Figure 3.30.Realization of ZRL in C2 form is summarized as following:Figure 3.30: Realization of ZRL in C2 form [53].As shall be discussed in Chapter 4, the Cauer method is adopted for the load circuit syn-thesis, which results in finding electrical parameters of the load.66Chapter 4SimulationIn this chapter, we describe simulations conducted on different combination of loads. Asexplained in Chapter 3, three major eigen-loads, heating (modelled as a resistive), light-ing (modelled as an inductive) and motor types (simulated as an equivalent model of aninduction motor) are considered in this work. We show that, with having time series infor-mation of the voltage and current of the total load, we are able to guess, the contributionsof the three eigen-loads, (1) heating load, (2) lighting load and (3) motor load in the totalload. As the result, we simulated all different topologies of the load circuit, consideringthe three mentioned eigen-loads. The significant strength of our method is that, it detectsbetween resistive, inductive and motor loads based on mostly steady Our algorithm de-veloped distinguishable steady state impedance coefficients for different load types. Forlower order loads, the steady state section of the voltage and current signals are sampled.While, for higher order loads, the transient section of the voltage and current signals aresample as well. The selection of sampling points, is important because, impedance coef-ficients are calculated for the sampled portion of voltage and current signals. Samplingrate should be sufficient to capture the essential characteristics of that circuit. Also, thesampling point distance, needs to be at least one cycle of 60 Hz period.PSCAD as an accurate reliable study tool is selected for this study for the following rea-sons:1. PSCAD is the most powerful power system transient simulation tool. For higherorder load circuits, we need accurate transient data for the voltage and the current.2. The main idea behind this research is taking advantage of the discrete time impedance67coefficients. These coefficients are derived from implementing the trapezoidal inte-gration rule on the voltage and current time series data. PSCAD is coded based ontrapezoidal method, i.e., the output data from PSCAD is a determined resource ofdata set, to validate our theory with.This work proves that by looking at the time series voltage and current signals, at realtime, we can say how many percent of the total load, is from motor, resistive and induc-tive type. Simulation based data are examined as described in Sections 4.2.1, 4.2.3, 4.2.4,4.2.6, 4.2.8, 4.2.10, 4.2.11, 4.2.12, and 4.2.14. Nine different load topologies, including dif-ferent combination of eigen-loads are simulated in PSCAD and the Output voltage andcurrent time series signals are imported into Matlab. The beauty of our algorithm is itscapability of deriving the electrical parameters of the load. this will exactly tell, how manypercent of the total impedance (or admittance) is from each of the eigen-loads. we employnetwork synthesis concept to achieve this goal. In order to calculate the exact values ofelectrical parameters for a load, we need to have whether the impedance or the admit-tance transfer function. Cauer (Foster) methods are stepped in to design a network whichrealizes the given Z(s) or Y(s).Note that we use admittance transfer function in developing our algorithm. We presenta method to construct a transfer function for the total impedance of the aggregated load.We have the discrete domain function of impedance by having the values of calculated co-efficients as will be discussed in Sections 4.2.1, 4.2.3, 4.2.4, 4.2.6, 4.2.8, 4.2.10, 4.2.11, 4.2.12,and 4.2.14. Obtained impedance coefficients, from the impedance transfer function asshown in Equation 4.1.Zz =an.Zn + an−1.Zn−1 + an−2.Zn−2 + · · ·+ a0bm.Zm + bm−1.Zm−1 + Bm−2.Zm−2 + · · ·+ b0 (4.1)The first step is to convert Equation 4.1 which is in Z domain to S domain as in Equa-tion 4.2.Zs =a′n.Sn + a′n−1.Sn−1 + a′n−2.Sn−2 + · · ·+ a′0b′m.Sm + b′m−1.Sm−1 + B′m−2.Sm−2 + · · ·+ b′0(4.2)Cauer method synthesizes admittance Laplacian function Ys = 1/Zs in the form of resis-tances and inductance (RL circuit)or RC circuit.684.1 Discrete-Continues TransformationIn this section, we briefly explain how to convert the derived impedance discrete do-main function to the continue domain function. There are five most common ”discrete-continuous” mapping techniques namely: backward difference method, forward differ-ence method, bilinear transform method, bilinear transform with pre-warping methodand matched-Z method. Since we applied trapezoidal to obtain our impedance discretecoefficients, we employ bilinear transform method.4.1.1 Bilinear Transform MethodThis method is also known as trapezoidal substitution method or Tustin method. Therelation between s and z is given by:s← 2T . z−1z+1 z = e(sT) = 1+sT/21−sT/2 where T is the time interval between samples of thediscrete-time system. The process consists of using a mapping function to replace ev-ery z in the F(z) with the function of s to obtain F(s). The highlighted portion of s plane ismapped inside the unit circle in z-plane as shown in Figure 4.1.Figure 4.1: Mapping of the s-plane to the z-plane with the bilinear transform method.The advantage of selected method (Tustin mapping)is that, it does not suffer stabilitylimitations associated with higher order and more complex numerical methods. There-fore, discrete-time derived coefficients(discrete-time digital filter) will be related to thecontinuous-time (analog filter)through the bilinear transform process. ”Network synthe-sis” block obtains analog coefficients. It applies whether 1st Cauer or 2nd Cauer methods,depending on the load topology determining the type of the loads and the exact values69of the load’s parameters. Sections 4.2.2, 4.2.2, 4.2.2, 4.2.2, 4.2.2, 4.2.2 will demonstrate theprocess of calculating each simulation load parameters.4.2 Simulation Results4.2.1 Simulation One: Parallel Resistive and Inductive LoadAs the first example, the simulation case is a parallel steady state RL circuit( representingResistive load in parallel with an Inductive load) with sampling rate of ∆t= 625µs (sam-pling frequency =1600 Hz) and maximum simulation time of tmax = 0.0625s. As stated,loads are classified into major three groups specified as resistor (R), lamps ( modelled asa single inductor), and motors ( model as an induction motor). This simulation case, de-scribes the behaviour of a resistive load paralleled with an inductive load. There are 100samples in total. Voltage magnitude is 230 KV to mimic the city electricity magnitude. Notransient sample is selected since the order of the circuit is one and the steady state por-tions of the voltage and current signals are sufficient to discriminate between inductanceand resistance values. Three samples are chosen from the voltage and the current signals.The reason is that this load case is a first order system. (Note that corresponding voltageand current meters are inserted right after the voltage source, because we need to getthe aggregated load’s information.)70Figure 4.2: Inductive load paralleled with a resistive load.As proved in Chapter 3, the magnitude of coefficients, a0, a1, · · · , an,b0,b1, · · · ,bn willbe equal regardless of the sampling points selection pattern. Therefore, for all nine simu-lation cases, we consider at least two different sampling scenario, (in some cases we illus-trated three different sampling scenarios). For each of the sampling schemes, a0, a1, · · · , anb0,b1, · · · ,bn coefficients are calculated. We prove that, as for the correct order of the loadcircuit, the magnitude of these coefficients will be equal. For the wrong choice of loadorder, the magnitude of these coefficients are diverse.Simulation One: First Sampling ScenarioIn this sampling scheme, the first sampling point is the 76th point of voltage and currentwaveforms. The distance between each samples is 10 samples. Matrix A describes thevalues for all current terms and all voltage history terms (refer to Equation 3.4 in Chap-ter 3). Matrix B contains voltage (Vt) values. First order impedance discrete time seriescoefficients are calculated and resulted as:[a0, a1,b0] = [12.3077− 12.3077− 0.2308]71Table 4.1: Values of current at time t and t − ∆t, voltage at time t − ∆t for the firstsampling scenario in the case of a parallel RL load.selected sample i(t) i(t− ∆t) v(t− ∆t)sample 76th -62.9551 -72.1146 -50.8833sample 86th 95.73191 89.20309 -191.188sample 96th -37.726 -19.3333 321.2645Table 4.2: Values of voltage at time t for the first sampling scenario in the case of aparallel RL load.selected sample v(t)sample 76th 124.475sample 86th 124.475sample 96th -300.509Simulation One: Second Sampling ScenarioIn this sampling scheme, the first sampling point is the 20th point of voltage and currentwaveforms. The distance between each samples is 20 samples. First order discrete timeTable 4.3: Values of current at time t and t− ∆t, voltage at time t− ∆t for the secondsampling scenario in the case of a parallel RL load.selected sample i(t) i(t− ∆t) v(t− ∆t)sample 20th -77.1399 -73.7391 147.6691sample 40th 14.39903 34.66287 289.8169sample 60th 97.46917 94.06845 -147.669impedance coefficients are calculated and resulted as:[a0, a1,b0] = [12.3077− 12.3077− 0.2308]For both sets of sampling schemes, derived coefficients {a0, a1,b0}, have the same val-ues. It confirms that this load case is a first order system. The same process is practicefor second order coefficients. Five samples are selected from the voltage and the currentsignals. [a0, a1, a2,b0,b1] coefficient values were different for two individual group of sam-pling points. It proves that, this load, does not full fill the second order requirements.72Table 4.4: Values of voltage at time t for the second sampling scenario in the case ofa parallel RL load.selected sample v(t)sample 20th -75.9326sample 40th -316.282sample 60th 75.93257Figure 4.3: First order impedance coefficients for the first and second scenarios of aparalleled RL load.4.2.2 Simulation One: Paralleled RL Load Parameters CalculationsFigure 4.4: Paralleled RL Load73Equation 4.3is resulted from the discrete time impedance coefficients that are derived fromPSCAD.V(z)I(z)= 12.3077I(z)− 12.3077I(z).z−1 − 0.2308V(z).z−1 TustinTrans f rom−−−−−−−−−→ Z(s) = s+ 200020s(4.3)Tustin method converts Equation 4.3 to Equation 4.4.Z(s) =s+ 200020s(4.4)Applying 2nd Cauer method (Continued fraction on Z(s)) determines load’s constituent’sparameters as following:Figure 4.5: Paralleled RL parameters calculation4.2.3 Simulation Two: Resistor in Series with a Parallel Inductance andCapacitanceFigure 4.6 is predicted to be a second order load based on calculations(refer to Chapter 3).We implement the same process as for Section 4.2.1 to calculate the discrete impedancecoefficients for this load as well. We selected 15 of the sample points from slow transientportion of the voltage and current signals. Note that slow transient is simulated withadding the ”ramp up time” value to the voltage source setting. The signal will reach to thesteady state condition after the ramp up time. Using a switch after the voltage source willalso produce the slow transient situation in PSCAD. In this simulation case, sampling rateis ∆t= 625µs(sampling frequency =1600 Hz), maximum simulation time is tmax = 0.0625s74and ramp up time=0.00625. Up to sample number 10, it is all slow transient portion ofthe signal. Steady state section starts from the 11th sample. One sample is selected fromthe slow transient samples. All the rest four samples are chosen from the steady statesamples. Note that, in Chapter 3, it is proved, knowing the circuit topology, we can guessthe order of the load circuit. Thereafter, network synthesis tool as explained in Chapter3,Section 3.5, will solve the load disaggregation problem. It will detect the load’s type anddetermine the loads’ parameters. But, In a real system, load’s circuit’s configuration isunknown, therefore the first and key step is to predict the order of the load’s circuit.Figure 4.6: Resistive load in series with a paralleled LC.Figure 4.6, represents a parallel inductance and capacitance which are placed in serieswith a resistance. Two different sampling scenarios are designated which both have thesame initial point, (both start from 8th point of the voltage/current signal )with the dif-ferent sample points distances.The process to determine the order of a load circuit is such that, (a0, a1, a2, · · · an,b0,b1,b2, · · · ,bn)coefficients are calculated for orders=(1,2,...,n) up to the correct order. for the correct pre-dicted order, resulting matrix of inv(A) ∗ B carries the same values for at least two differ-ent sampling patterns.Simulation Two: First Sampling ScenarioIn this sampling scheme, the first sampling point is the 8th point of voltage and currentwaveforms. The distance between each samples is 20 samples. First order discrete time75Table 4.5: Values of current at time t and t − ∆t, voltage at time t − ∆t for the firstsampling scenario in the case of a resistance in series with a paralleled LC load.selected sample i(t) i(t− ∆t) v(t− ∆t)8th sample 0.221231 1.549666 -30.5328th sample 15.56777 14.67393 -321.26548th sample -2.20467 -5.80907 50.8833Table 4.6: Values of voltage at time t for the first sampling scenario in the case of aresistance in series with a paralleled LC load.selected sample v(t)8th sample -17.864228th sample 324.266448th sample 25.52032impedance coefficients are calculated and resulted as:[a0, a1,b0] = [36.3046− 19.3411− 0.1335]Simulation Two: Second Sampling ScenarioIn this sampling scheme, the first sampling point is the 8th point of voltage and currentwaveforms. The distance between each samples is 35 samples.Table 4.7: Values of current at time t and t− ∆t, voltage at time t− ∆t for the secondsampling scenario in the case of a resistance in series with a paralleled LC load.selected sample i(t) i(t− ∆t) v(t− ∆t)8th sample 0.221231 1.549666 -30.5343th sample -15.3608 -15.8779 316.281978th sample 9.564251 6.340717 -191.18876Table 4.8: Values of voltage at time t for the first sampling scenario in the case of aresistance in series with a paralleled LC load.selected sample v(t)8th sample -17.86443th sample -289.81678th sample 247.336Discrete time impedance coefficients are calculated and resulted as:[a0, a1,b0] = [36.079 − 19.0608 − 0.1209] So far, we calculated the first order impedancecoefficients for two different sampling scenarios. The magnitude of the two impedancematrices are equal within less than 10% error. Since this load satisfies the 1st order criteria,we will continue to the second order calculation. As explained in Chapter 3, for the secondorder situation, we will engage five sampling points from the voltage and the currentsignals.Simulation Two: First Sampling Scenario For the Second Order StatusIn this sampling scheme, the first sampling point is the 8th point of voltage and currentwaveforms. The distance between each samples is 20 samples.Table 4.9: Values of current at time t, t−∆t and t− 2∆t, voltage at time t−∆t and t−2∆t for the first sampling scenario in the case of a resistance in series with a paralleledLC load.selected sample i(t) i(t− ∆t) i(t− 2∆t) v(t− ∆t) v(t− 2∆t)8th sample 0.221231 1.549666 2.336653 -30.53 -62.237628th sample 15.56777 14.67393 13.02486 -321.265 -300.50948th sample -2.20467 -5.80907 -9.0913 50.8833 124.475168th sample -15.7265 -14.7778 -13.0126 321.2645 300.509568th sample 2.200985 5.811509 9.100879 -50.8833 -124.47577Table 4.10: Values of voltage at time t for the first sampling scenario in the case of aparallel RL load.selected sample v(t)8th sample -17.864228th sample 324.266448th sample 25.5203268th sample -324.26688th sample -25.5203Discrete time impedance coefficients are: [a0, a1, a2,b0,b1] = [22.8470− 32.8826 17.1530 −1.6441 1.0000]Simulation Two: Second Sampling Scenario For the Second Order StatusIn this sampling scheme, the first sampling point is the 8th point of voltage and currentwaveforms. The distance between each samples is 10 samples.Table 4.11: Values of current at time t, t− ∆t and t− 2∆t, voltage at time t− ∆t andt − 2∆t for the second sampling scenario in the case of a resistance in series with aparalleled LC load.selected sample i(t) i(t− ∆t) i(t− 2∆t) v(t− ∆t) v(t− 2∆t)8th sample 0.221231 1.549666 2.336653 -30.53 -62.237618th sample -13.4113 -15.2336 -15.9516 263.1482 300.509528th sample 15.56777 14.67393 13.02486 -321.265 -300.50938th sample -9.59291 -6.34855 -2.74361 191.1884 124.475148th sample -2.20467 -5.80907 -9.0913 50.8833 124.475178Table 4.12: Values of voltage at time t for the second sampling scenario in the case ofa resistance in series with a paralleled LC load.selected sample v(t)8th sample -17.864218th sample -211.24528th sample 324.266438th sample -247.33748th sample 25.52032Discrete time impedance coefficients are calculated and resulted as:[a0, a1, a2,b0,b1] = [22.8470 − 32.8826 17.1530 − 1.6441 1.0000] Second order de-rived impedance coefficients are identical. It justifies that this circuit satisfies the secondorder characteristics. Up to this end, first order and 2nd order discrete time coefficientsvalues for different group of sampling points are equivalent. This forces us to continuethe same process for third order.Simulation Two: First Sampling Scenario For the Third Order StatusIn this sampling scheme, the first sampling point is the 5th point of voltage and currentwaveforms. The distance between each samples is 15 samples. For the third order situa-tion, we need at least 7 sampling points.Table 4.13: Values of current at time t, t − ∆t, t − 2∆t and t − 3∆t, voltage at timet− ∆t, t− 2∆t and t− 3∆t for the first sampling scenario in the case of a resistance inseries with a paralleled LC load.selected sample i(t) i(t− ∆t) i(t− 2∆t) i(t− 3∆t) v(t− ∆t) v(t− 2∆t) v(t− 3∆t)5th sample 2.688084 2.664771 2.253407 1.384349 -74.201 -57.9634 -31.628220th sample -7.19437 -10.6446 -13.4113 -15.2336 147.6691 211.2454 263.148235th sample 1.014964 4.708302 8.121243 11.05883 -25.5203 -100.514 -169.95350th sample 5.171216 1.524038 -2.20467 -5.80907 -100.514 -25.5203 50.883365th sample -10.5285 -7.46292 -3.98515 -0.28728 211.2454 147.6691 75.9325780th sample 14.27682 12.25925 9.564251 6.340717 -289.817 -247.337 -191.18895th sample -15.8509 -15.1888 -13.6873 -11.4295 324.2664 309.3493 277.337579Table 4.14: Values of voltage at time t for the first sampling scenario in the case of aresistance in series with a paralleled LC load.selected sample v(t)5th sample 76.4753620th sample -75.932635th sample -50.883350th sample 169.952665th sample -263.14880th sample 316.281995th sample -321.265Discrete time impedance coefficients are calculated and resulted as:[a0, a1, a2, a3,b0,b1,b2] = [537.1875 − 775.8906 407.2500 − 2.0454 37.210 − 61.4375 32.0879]Simulation Two: Second Sampling Scenario For the Third Order StatusIn this sampling scheme, the first sampling point is the 5th point of voltage and currentwaveforms. The distance between each samples is 10 samples.Table 4.15: Values of current at time t, t − ∆t, t − 2∆t and t − 3∆t, voltage at timet−∆t, t− 2∆t and t− 3∆t for the second sampling scenario in the case of a resistancein series with a paralleled LC load.selected sample i(t) i(t− ∆t) i(t− 2∆t) i(t− 3∆t) v(t− ∆t) v(t− 2∆t) v(t− 3∆t)5th sample 2.688084 2.664771 2.253407 1.384349 -74.201 -57.9634 -31.628215th sample -15.5555 -14.2037 -12.1923 -9.87563 324.2664 309.3493 277.337525th sample 10.6904 7.756909 4.333147 0.562823 -211.245 -147.669 -75.932635th sample 1.014964 4.708302 8.121243 11.05883 -25.5203 -100.514 -169.95345th sample -11.8722 -14.0006 -15.3608 -15.8779 247.3366 289.8169 316.281955th sample 15.85452 15.19178 13.68826 11.42731 -324.266 -309.349 -277.33865th sample -10.5285 -7.46292 -3.98515 -0.28728 211.2454 147.6691 75.9325780Table 4.16: Values of voltage at time t for the second sampling scenario in the case ofa resistance in series with a paralleled LC load.selected sample v(t)5th sample 76.4753615th sample -321.26525th sample 263.148235th sample -50.883345th sample -191.188455th sample 321.264565th sample -263.148Discrete time impedance coefficients are calculated and resulted as:[a0, a1, a2, a3,b0,b1,b2] = [6.9746 − 38.082 45.5977 − 21.0547 − 4.0745 4.9507 − 2.2212]Dif-ferent obtained impedance coefficient matrices confirm that this load is not a third orderload. The largest order, which resulted to the equal impedance values is two. Note that,Both inductance and capacitance, generate first order history terms and they are set inparallel. Second history terms will be generated when we implement trapezoidal inte-gration on Zz. We analyse this simulation load case with a different Deltat value as well.The new value of ∆t is 0.0022 seconds. tmax = 0.4s and ramp up time=0.044(slow transientends at sample 22). [a0, a1, a2,b0,b1] = [24.9774 3.8009 15.0226 0.1900 1.0000] For the sec-ond order conditions, we considered two different sampling scenarios. The two scenariosresulted in two identical impedance matrices for t = t, t − ∆t, t − 2∆t time stamps. thisproves our statement in Chapter 3 that says, different ∆t values will result in equal discreteimpedance coefficients. In other word, the sampling rate will not impact our algorithm.81Figure 4.7: First order impedance coefficients for the first and second scenarios of aresistance in series with a paralleled LC load.Figure 4.8: Second order impedance coefficients for the first and second scenarios ofa resistance in series with a paralleled LC load.82Figure 4.9: Third order impedance coefficients for the first and second scenarios of aresistance in series with a paralleled LC load.Figures 4.7 and 4.7 demonstrate first order and 2nd order impedance coefficient valuesfor both scenarios. This confirms that two graphs are exactly matched while in Figure 4.9,we see two different patterns of coefficients.4.2.4 Simulation Three: Motor LoadFigure 4.10: Motor load case.83This simulation case describes the behaviour of a motor type load. We assume motor loadsas the second order systems based on calculations in Chapter 3. Simulation data confirmsthis fact. In this simulation case, sampling rate is ∆t= 2200µs, maximum simulation timeis tmax = 0.4seconds and ramp up time is 0.066 seconds. Following analysis starts withorder=1. It calculates (a0, a1,b0) coefficients for at least two different sets of voltage andcurrent samples. the same process is continued for higher orders until we get the correctorder of the circuit. We selected one sampling point from transient portion of the voltageand the current signal.Simulation Three: First Sampling ScenarioIn this sampling scheme, the first sampling point is the 25th point of voltage and currentwaveforms. The distance between each samples is 20 samples.Table 4.17: Values of current at time t and t− ∆t, voltage at time t− ∆t for the firstsampling scenario for a single motor loadselected sample i(t) i(t− ∆t) v(t− ∆t)25th sample 13.1152 7.969328 -11.622745th sample -11.9759 -16.7281 6.97520665th sample -1.39197 11.37702 6.267753Table 4.18: Values of voltage at time t for the first sampling scenario for a single motorload.selected sample v(t)25th sample 6.13203145th sample 5.54446865th sample -14.7334Discrete time impedance coefficients are calculated and resulted as:[a0, a1,b0] = [2.3049− 1.56431.0007]84Simulation Three: Second Sampling ScenarioIn this sampling scheme, the first sampling point is the 25th point of voltage and currentwaveforms. The distance between each samples is 32 samples.Table 4.19: Values of current at time t and t−∆t, voltage at time t−∆t for the secondsampling scenario for a single motor loadselected sample i(t) i(t− ∆t) v(t− ∆t)25th sample 13.1152 7.969328 -11.622757th sample 4.449092 14.92126 0.97679189th sample -15.2177 -5.10202 15.47781Table 4.20: Values of voltage at time t for the second sampling scenario for a singlemotor load.selected sample v(t)25th sample 6.13203157th sample -12.1189th sample -11.6041Discrete coefficients are calculated and resulted as:[a0, a1,b0] = [2.1643− 1.51450.8789]Simulation Three: Third Sampling ScenarioIn this sampling scheme, the first sampling point is the 40th point of voltage and currentwaveforms. The distance between each samples is 32 samples.Table 4.21: Values of current at time t and t− ∆t, voltage at time t− ∆t for the thirdsampling scenario for a single motor loadselected sample i(t) i(t− ∆t) v(t− ∆t)40th sample 16.33137 8.214105 -15.477872th sample 6.438923 15.75991 -0.97679104th sample -14.2184 -3.06126 15.1600885Table 4.22: Values of voltage at time t for the third sampling scenario for a singlemotor load.selected sample v(t)40th sample 9.30121972th sample -10.7907104th sample -12.8112Discrete coefficients are calculated and resulted as:[a0, a1,b0] = [2.1643− 1.51450.8789] Looking at the three obtained matrices for [a0, a1,b0],all are equivalent. We will continue with second order analysis:Simulation Three: First Sampling Scenario For the Second Order StatusIn this sampling scheme, the first sampling point is the 25th point of voltage and currentwaveforms. The distance between each samples is 20 samples.Table 4.23: Values of current at time t, t− ∆t and t− 2∆t, voltage at time t− ∆t andt− 2∆t for the first sampling scenario for a single motor load.selected sample i(t) i(t− ∆t) i(t− 2∆t) v(t− ∆t) v(t− 2∆t)25th sample 13.1152 7.969328 -1.69501 -11.6227 -9.3948545th sample -11.9759 -16.7281 -10.6144 6.975206 14.9656465th sample -1.39197 11.37702 16.75935 6.267753 -6.2677585th sample 13.75697 2.231712 -10.7419 -14.9656 -6.97521105th sample -16.1424 -14.2184 -3.06126 12.81116 15.1600886Table 4.24: Values of voltage at time t for the first sampling scenario in the case of asingle motor loadselected sample v(t)25th sample 6.13203145th sample 5.54446865th sample -14.733485th sample 13.23837105th sample -2.14351Discrete time series coefficients are calculated and resulted as:[a0, a1, a2,b0,b1] = [2.2981− 3.8341 1.5405 0.0112− 0.9888]Simulation Three: Second Sampling Scenario For the Second Order StatusIn this sampling scheme, the first sampling point is the 25th point of voltage and currentwaveforms. The distance between each samples is 32 samples.Table 4.25: Values of current at time t, t− ∆t and t− 2∆t, voltage at time t− ∆t andt− 2∆t for the second sampling scenario for a single motor load.selected sample i(t) i(t− ∆t) i(t− 2∆t) v(t− ∆t) v(t− 2∆t)25th sample 13.1152 7.969328 -1.69501 -11.6227 -9.3948557th sample 4.449092 14.92126 15.70543 0.976791 -10.790789th sample -15.2177 -5.10202 8.327259 15.47781 9.301219121th sample -9.39709 -16.5789 -12.9948 4.057742 13.81616153th sample 12.16258 -0.2891 -12.5526 -14.1579 -4.8071887Table 4.26: Values of voltage at time t for the second sampling scenario in the case ofa single motor loadselected sample v(t)25th sample 6.13203157th sample -12.1189th sample -11.6041121th sample 8.335509153th sample 14.31545Discrete time series coefficients are calculated and resulted as:[a0, a1, a2,b0,b1] = [2.2980− 3.83391.54040.0112− 0.9887]Simulation Three: Third Sampling Scenario For the Second Order StatusIn this sampling scheme, the first sampling point is the 15th point of voltage and currentwaveforms. The distance between each samples is 40 samples.Table 4.27: Values of current at time t, t− ∆t and t− 2∆t, voltage at time t− ∆t andt−−2∆t for the third sampling scenario for a single motor load.selected sample i(t) i(t− ∆t) i(t− 2∆t) v(t− ∆t) v(t− 2∆t)15th sample -3.3824 -6.74192 -5.57201 1.429158 5.37579855th sample 15.70543 6.29245 -7.2055 -15.5514 -10.214195th sample 2.804772 14.07973 16.21294 2.530041 -9.61164135th sample -16.7539 -11.5663 1.132217 14.60327 13.81616175th sample 3.476008 -9.74306 -16.6352 -8.0028 4.43385988Table 4.28: Values of voltage at time t for the third sampling scenario in the case of asingle motor loadselected sample v(t)15th sample 4.19296455th sample 10.7907195th sample -13.0289135th sample -5.90798175th sample 15.24297Discrete time series coefficients are calculated and resulted as:[a0, a1, a2,b0,b1] = [2.2981− 3.83411.54050.0112− 0.9888]Motor circuit full fills 2nd order principals because all three matrices[A = 2.2981 −3.8341 1.5405 0.0112 −0.9888][B = 2.2980 −3.8339 1.5404 0.0112 −0.9887][C = 2.2980 −3.8339 1.5404 0.0112 −0.9887]are equivalent as shown in figure 4.12.Additionally, since it is determined that motor circuit is satisfying the second order charac-teristics, we analyse the case with the third order criteria. It resulted in different impedancecoefficients. This interprets that third order and higher orders cannot be the correct orderfor the single motor load.89Figure 4.11: First order impedance coefficients for three different scenarios of a singlemotor load.Figure 4.12: Second order impedance coefficients for three different scenarios of asingle motor load.904.2.5 Simulation Three: Single Motor Load Parameters CalculationsFigure 4.13: Single motor load.Y(s)is the results from converting the discrete time coefficients to the continuous domain.Figure 4.14: Motor load parameter calculation.As shown in Figure 4.14, in every step of the division, the quotients are R and L valuesof the load circuit. The first quotient is Zs which is a series arm and is equivalent to R1.Second quotient is Ys which indicates shunt arm and is equivalent to Lm. Third quotientshows R2 and the last quotient is the value of L2. L1 value is not synthesized.91We captured four parameters of a single motor load. Note that we have to continue thedivision as yet we get zero in the remainder.4.2.6 Simulation Four: Paralleled Resistive and a MotorThis simulation, evaluates the behaviour of an aggregated load including a single motorload in parallel with a resistive load. We will consider all combinations of heating (resis-tive), lighting (inductive) and motor loads , because these are the three major eigneloadswe considered in this work. We will prove that we can distinguish between these threetypes In this simulation case, sampling rate is ∆t = 2200µs(sampling frequency =1600Figure 4.15: Paralleled resistive and a motor.Hz), maximum simulation time is tmax = 0.8s(totalsamples : 367) and ramp up time=0.08.Theoretically, this circuit is also supposed to be 2nd order. Because, resistance does notgenerate any history terms. Motor load is also confirmed to be a second order system inSection ??. It may seem that this algorithm is not able to discern between a case of thesingle motor load and a case of the paralleled motor and resistance since both systems aresecond order. But in fact, due to different obtained discrete coefficients, network synthesiswill result in two different topologies which determines whether there is a single motoror a paralleled motor and resistance.92Simulation Four: First Sampling ScenarioIn this sampling scheme, the first sampling point is the 25th point of voltage and currentwaveforms. The distance between each samples is 20 samples.Table 4.29: Values of current at time t and t− ∆t, voltage at time t− ∆t for the firstsampling scenario in the case of a paralleled resistance and a motor.selected sample i(t) i(t− ∆t) v(t− ∆t)25th sample 11.07299 7.054133 -9.5887545th sample -11.6826 -17.0501 6.97520665th sample -2.13397 11.05827 6.267753Table 4.30: Values of voltage at time t for the first sampling scenario in the case of aparalleled resistance and a motor.selected sample v(t)25th sample 5.05892545th sample 5.54446865th sample -14.7334Discrete time series coefficients are calculated and resulted as:[a0, a1,b0] = [2.0939− 1.4113 0.8522]Simulation Four: Second Sampling ScenarioIn this sampling scheme, the first sampling point is the 25th point of voltage and currentwaveforms. The distance between each samples is 32 samples.Table 4.31: Values of current at time t and t−∆t, voltage at time t−∆t for the secondsampling scenario in the case of a paralleled resistance and a motor.selected sample i(t) i(t− ∆t) v(t− ∆t)25th sample 11.07299 7.054133 -9.5887557th sample 3.838176 14.86706 0.97679189th sample -15.8025 -5.88056 15.4778193Table 4.32: Values of voltage at time t for the second sampling scenario in the case ofa paralleled resistance and a motor.selected sample v(t)25th sample 5.05892557th sample -12.1189th sample -11.6041Discrete time series coefficients are calculated and resulted as:[a0, a1,b0] = [2.0772− 1.4058 0.8369]Simulation Four: Third Sampling ScenarioIn this sampling scheme, the first sampling point is the 30th point of voltage and currentwaveforms. The distance between each samples is 40 samples.Table 4.33: Values of current at time t and t− ∆t, voltage at time t− ∆t for the thirdsampling scenario in the case of a paralleled resistance and a motor.selected sample i(t) i(t− ∆t) v(t− ∆t)30th sample -7.47218 -12.8154 3.98664270th sample 15.44856 5.059682 -15.3798110th sample 4.256216 15.0745 0.586322Table 4.34: Values of voltage at time t for the third sampling scenario in the case of aparalleled resistance and a motor.selected sample v(t)30th sample 5.83964270th sample 12.11003110th sample -11.8608Discrete time series coefficients are calculated and resulted as:[a0, a1,b0] = [2.1281− 1.4220 0.8824]94Considering three calculated coefficient matrices:[A = 2.0667 −1.4027 0.8271][B = 2.0717 −1.4043 0.8317][C = 2.1110 −1.4168 0.8672]They are all equivalent. We will continue with order=2.Simulation Four: First Sampling Scenario For The Second Order ConditionIn this sampling scheme, the first sampling point is the 25th point of voltage and currentwaveforms. The distance between each samples is 40 samples.Table 4.35: Values of current at time t, t− ∆t and t− 2∆t, voltage at time t− ∆t andt− 2∆tfor the first sampling scenario in the case of a parallel resistance and a motorresistanceselected sample i(t) i(t− ∆t) i(t− 2∆t) v(t− ∆t) v(t− 2∆t)25th sample 11.07299 7.054133 -1.01084 -9.58875 -7.7507565th sample -2.13397 11.05827 17.06735 6.267753 -6.26775105th sample -16.2537 -14.8632 -3.82351 12.81116 15.16008145th sample 8.217491 -5.49593 -15.6423 -11.0689 0.586322185th sample 13.16801 16.91667 9.679471 -8.66295 -15.3798Table 4.36: Values of voltage at time t for the first sampling scenario in the case of aparalleled resistance and a motorselected sample v(t)25th sample 5.05892565th sample -14.7334105th sample -2.14351145th sample 15.5367185th sample -3.67906Discrete time series coefficients are calculated and resulted as:[a0, a1, a2,b0,b1] = [2.0612− 3.4389 1.3817− 0.1619− 0.8178]95Simulation Four: Second Sampling Scenario For The Second Order ConditionIn this sampling scheme, the first sampling point is the 15th point of voltage and currentwaveforms. The distance between each samples is 50 samples.Table 4.37: Values of current at time t, t− ∆t and t− 2∆t, voltage at time t− ∆t andt − 2∆tfor the second sampling scenario in the case of a parallel resistance and amotor resistance.selected sample i(t) i(t− ∆t) i(t− 2∆t) v(t− ∆t) v(t− 2∆t)15th sample -2.61752 -5.62104 -4.81866 1.179055 4.43503365th sample -2.13397 11.05827 17.06735 6.267753 -6.26775115th sample 11.68301 -1.30108 -13.4423 -13.4395 -3.29806165th sample -16.7807 -8.9643 4.671463 15.47781 11.60412215th sample 15.46045 15.79731 5.875347 -11.6041 -15.4778Table 4.38: Values of voltage at time t for the second sampling scenario in the case ofa paralleled resistance and a motor.selected sample v(t)15th sample 3.45919565th sample -14.7334115th sample 14.85421165th sample -9.30122215th sample 0.195482Impedance coefficients are calculated and resulted as:[a0, a1, a2,b0,b1] = [2.0612− 3.4389 1.3817− 0.1618− 0.8178]Simulation Four: Third Sampling Scenario For The Second Order ConditionIn this sampling scheme, the first sampling point is the 20th point of voltage and currentwaveforms. The distance between each samples is 50 samples.96Table 4.39: Values of current at time t, t− ∆t and t− 2∆t, voltage at time t− ∆t andt− 2∆tfor the third sampling scenario in the case of a parallel resistance and a motorresistance.selected sample i(t) i(t− ∆t) i(t− 2∆t) v(t− ∆t) v(t− 2∆t)20th sample -4.67116 2.140791 6.900771 5.479104 -0.2741170th sample 15.44856 5.059682 -8.6172 -15.3798 -8.66295120th sample -16.7856 -13.6895 -1.70625 13.81616 14.60327170th sample 11.70028 17.07949 11.367 -6.97521 -14.9656220th sample -2.15386 -13.9538 -16.6941 -2.53004 9.61164Table 4.40: Values of voltage at time t for the third sampling scenario in the case of aparalleled resistance and a motor.selected sample v(t)20th sample -8.1179370th sample 12.11003120th sample -4.05774170th sample -5.54447220th sample 13.02888Discrete time series coefficients are calculated and resulted as:[a0, a1, a2,b0,b1] = [2.0604− 3.4383 1.3819− 0.1629− 0.8171] Motor in parallel with Resis-tance also full fills second order principals, because all three obtained matrices:[A = 2.0615 −3.4394 1.3820 −0.1618 −0.8180][B = 2.0613 −3.4396 1.3823 −0.1622 −0.8179][C = 2.0652 −3.4422 1.3811 −0.1566 −0.8214]are equivalent. This is shown in Figures 4.16 and 4.17.97Figure 4.16: First order impedance coefficients for three scenarios of a paralleled mo-tor and resistanceFigure 4.17: Second order impedance coefficients for three scenarios of a paralleledmotor and resistance984.2.7 Simulation One: Paralleled Motor and Resistance Load ParametersCalculationsFigure 4.18: Paralleled motor and resistance.This case determines a motor type in parallel with resistance. As described in Section 4.2.6,this load is a second order load.Figure 4.19: Motor load in parallel with a resistive load parameter calculation.Note that, PSCAD discrete time coefficients are different for a Motor load and a Mo-tor+Resistance load cases, although, the circuit order and the number of coefficients arethe same.99Following case products five steps of division prior to getting a remainder of zero whereasin Section 4.2.5, we get zero in remainder after continuing the long division for four con-secutive stages. The first four quotients, determine motor’s parameters, R1,Lm,R2&L2and the fifth quotient is the value of parallel resistance.4.2.8 Simulation Five: Single Motor in Parallel with One Resistance and OneInductanceFigure 4.20: Single motor in parallel with one resistance and one inductance load.In this simulation case, sampling rate is ∆t = 2200µs(sampling frequency =1600 Hz),maximum simulation time is tmax = 0.8s and ramp up time=0.132s. This case describes anaggregated load combining of motor, inductance and resistance. As discussed in Chap-ter 3, induction motor load is a 2nd order system. Inductance is a first order system.Therefore, a parallel circuit including both will have the characteristics of a third ordersystem. We calculated impedance discrete time coefficients for first and second order sit-uation. Resulted discrete time coefficients were within 10% equivalence. Thus, we needto move on to the third order analysis. In Following you will see the results earned fromorder=3 analysis.100Simulation Five: First Sampling Scenario For The Third Order StatusIn this sampling scheme, the first sampling point is the 35th point of voltage and currentwaveforms. The distance between each samples is 10 samples.Table 4.41: Values of current at time t, t − ∆t, t − 2∆t and t − 3∆t, voltage at timet− ∆t, t− 2∆t and t− 3∆t for the first sampling scenario in the case of a motor loadin parallel with an inductance and a resistance.selected sample i(t) i(t− ∆t) i(t− 2∆t) i(t− 3∆t) v(t− ∆t) v(t− 2∆t) v(t− 3∆t)35th sample -3.6274 4.025515 8.683997 7.590059 5.286402 -1.34935 -6.7315945th sample -8.59834 -12.2633 -7.98461 1.151465 4.998898 10.47595 9.04621855th sample 14.64574 6.203539 -5.87806 -13.7587 -13.7371 -8.8522 1.49228865th sample -1.65254 11.63228 17.44466 12.04937 6.267753 -6.26775 -14.733475th sample -14.8389 -16.5727 -7.54313 6.387411 10.21408 15.55144 10.7907185th sample 14.35952 2.584053 -10.8687 -17.2634 -14.9656 -6.97521 5.54446870th sample 2.614842 14.37684 16.80403 8.320351 2.530041 -9.61164 -15.5122Table 4.42: Values of voltage at time t for the first sampling scenario in the case of amotor load in parallel with an inductance and a resistance.selected sample v(t)35th sample -8.7902245th sample 4.06594355th sample 9.71164265th sample -14.733475th sample 1.75563385th sample 13.2383770th sample -13.0289Discrete time series coefficients are calculated and resulted as:A=[2.2080 −4.3249 2.5497 −0.4297 −0.2943 −0.9923 0.3021]101Simulation Five: Second Sampling Scenario For The Third Order StatusIn this sampling scheme, the first sampling point is the 20th point of voltage and currentwaveforms. The distance between each samples is 15 samples.Table 4.43: Values of current at time t, t − ∆t, t − 2∆t and t − 3∆t, voltage at timet − ∆t, t − 2∆t and t − 3∆t for the second sampling scenario in the case of a motorload in parallel with an inductance and a resistance.selected sample i(t) i(t− ∆t) i(t− 2∆t) i(t− 3∆t) v(t− ∆t) v(t− 2∆t) v(t− 3∆t)20th sample -2.73869 1.437729 4.278048 4.175772 3.320669 -0.16612 -3.1628935th sample -3.6274 4.025515 8.683997 7.590059 5.286402 -1.34935 -6.7315950th sample -3.44124 7.511785 13.2344 10.30246 6.402241 -3.47319 -10.726865th sample -1.65254 11.63228 17.44466 12.04937 6.267753 -6.26775 -14.733480th sample 0.451432 13.05818 17.18663 10.15616 4.433859 -8.0028 -15.24395th sample 2.614842 14.37684 16.80403 8.320351 2.530041 -9.61164 -15.5122110th sample 4.737216 15.46906 16.15683 6.353933 0.586322 -11.0689 -15.5367Table 4.44: Values of voltage at time t for the first sampling scenario in the case of amotor load in parallel with an inductance and a resistance.selected sample v(t)20th sample -4.9199635th sample -8.7902250th sample -12.448465th sample -14.733480th sample -13.991595th sample -13.0289110th sample -11.8608Discrete time series coefficients are calculated and resulted as:B=[2.2080 −4.3250 2.5498 −0.4298 −0.2944 −0.9923 0.3021]Matrices A and B equal. Thus, third order properties match with this load’s properties.Worth mentioning, for higher order systems, to be able to detect the correct load’s com-ponents, more samples from the slow transient section of the voltage and the current102signals need to be selected.Order=4 is also verified not to be the correct order for this load circuit, since it results innon-equal impedance coefficients, for the two different sampling scenarios as following:Simulation Five: First Sampling Scenario For The Fourth Order StatusIn this sampling scheme, the first sampling point is the 30th point of voltage and currentwaveforms. The distance between each samples is 10 samples. Impedance coefficientsare:A=[1.4781 −1.8820 −0.2716 0.8721 −0.1931 −0.6313 0.4679 −1.1464 0.4666]Simulation Five: Second Sampling Scenario For The Fourth Order StatusIn this sampling scheme, the first sampling point is the 20th point of voltage and currentwaveforms. The distance between each samples is 15 samples.B=[0.8023 −3.2489 4.2122 −2.0933 0.3266 −2.5863 2.2977 −0.8562 0.4071]Figure 4.21: Third order impedance coefficients for three different scenarios of a mo-tor load in parallel with an inductance and a resistance.103Figure 4.22: Fourth order impedance coefficients for three different scenarios of a mo-tor load in parallel with an inductance and a resistance.4.2.9 Simulation Five: Single Motor in Parallel with One Resistance and OneInductance Parameters CalculationsFigure 4.23: Single motor in parallel with one resistance and one inductanceThis case determines an induction motor type in parallel with an inductance. FollowingYs is derived from the third order discrete time coefficients which have been analysedin Section 4.2.8. Figure 4.24 indicates six stages of continued division, prior to getting aremainder of zero.104Figure 4.24: Motor in parallel with an inductive and a resistive load parameter calcu-lation.The first four quotients correspond to the motor’s parameters, R1,Lm,R2&L2 and 5thand 6th quotients, are the values of individual resistance and inductance.1054.2.10 Simulation Six: Two Series Motor in Parallel with a Resistive and anInductive LoadFigure 4.25: Two series motor in parallel with a resistive and an inductive load.In this simulation case, sampling rate is ∆t = 625µs(sampling frequency =1600 Hz), max-imum simulation time is tmax = 0.1875s, ramp up time=0.0375s (slow transient ends atsample number 60). This simulation describes two series motor loads in parallel with aresistance and an inductance. Although, in real power system, loads are connected in par-allel from a distribution feeder. But we analysed a higher order case in this section. Some-times, some power electronic loads have more complex circuit structures which leads tothe higher orders. We proved that our algorithm is capable of detecting the higher orderloads as well.Simulation Six: First Sampling Scenario For The Third Order StatusIn this sampling scheme, the first sampling point is the 20th point of voltage and currentwaveforms. The distance between each samples is 10 samples.106Table 4.45: Values of current at time t, t − ∆t, t − 2∆t and t − 3∆t, voltage at timet − ∆t, t − 2∆t and t − 3∆t for the first sampling scenario in the case of two seriesmotor in parallel with a resistive and an inductive Load.selected sample i(t) i(t− ∆t) i(t− 2∆t) i(t− 3∆t) v(t− ∆t) v(t− 2∆t) v(t− 3∆t)20th sample -1.08422 -1.18047 -1.20521 -1.1643 2.11873 2.862528 3.35609430th sample 1.847889 1.670556 1.409224 1.084638 -6.90432 -6.97878 -6.6580940th sample -1.63121 -1.07295 -0.48892 0.086761 8.778512 7.294637 5.48627650th sample -0.58222 -1.30023 -1.92157 -2.41571 -3.84574 -0.95609 1.86572160th sample 3.152638 3.541393 3.719616 3.684601 -6.82702 -9.59789 -11.746370th sample -4.04172 -3.75978 -3.27442 -2.6127 14.79497 15.50839 15.3648280th sample 2.216767 1.375249 0.45409 -0.49584 -13.8608 -11.8291 -9.14379Table 4.46: Values of voltage at time t for the first sampling scenario in the case oftwo series motor in parallel with a resistive and an inductive Load.selected sample v(t)20th sample -1.1499930th sample 6.41091840th sample -9.8322450th sample 6.63800660th sample 3.57103270th sample -13.26480th sample 15.12653Discrete time series coefficients are calculated and resulted as:A=[10.7906 −28.2620 24.2639 −6.7923 −1.7882 0.6953 0.0954]Simulation Six: Second Sampling Scenario For The Third Order StatusIn this sampling scheme, the first sampling point is the 10th point of voltage and currentwaveforms. The distance between each samples is 15 samples.107Table 4.47: Values of current at time t, t − ∆t, t − 2∆t and t − 3∆t, voltage at timet− ∆t, t− 2∆t and t− 3∆t for the second sampling scenario in the case of two seriesmotor in parallel with a resistive and an inductive Load.selected sample i(t) i(t− ∆t) i(t− 2∆t) i(t− 3∆t) v(t− ∆t) v(t− 2∆t) v(t− 3∆t)10th sample -0.07954 0.03305 0.109863 0.15053 0.640957 0.142396 -0.2433525th sample 0.339286 -0.03408 -0.37887 -0.67713 -3.87283 -2.58956 -1.2710540th sample -1.63121 -1.07295 -0.48892 0.086761 8.778512 7.294637 5.48627655th sample 3.023137 2.447016 1.754769 0.989035 -13.6991 -12.8223 -11.274470th sample -4.04172 -3.75978 -3.27442 -2.6127 14.79497 15.50839 15.3648285th sample 3.889271 3.97408 3.835769 3.481958 -11.8291 -13.8608 -15.1265100th sample -3.39011 -3.83106 -4.0635 -4.07461 7.062435 10.10304 12.58535Table 4.48: Values of voltage at time t for the second sampling scenario in the case oftwo series motor in parallel with a resistive and an inductive Load.selected sample v(t)10th sample -1.2192325th sample 5.0341440th sample -9.8322455th sample 13.8283470th sample -13.26485th sample 9.143793100th sample -3.63156Discrete time series coefficients are calculated and resulted as:B=[10.7411 −27.0579 22.4343 −6.1140 −1.7042 0.6530 0.0644]It is important to note that, since this circuit is a higher order circuit, more transient sam-ples are selected in compare with the lower order loads. In this specific example, foursampling points are selected from point 1 to 60, which are considered as the transientdata. We will analyse 4th order coefficients since third order coefficients satisfy our as-sumptions.108Simulation Six: First Sampling Scenario For The Fourth Order Statusn this sampling scheme, the first sampling point is the 20th point of voltage and currentwaveforms. The distance between each samples is 10 samples.Table 4.49: Values of current at time t, t− ∆t, t− 2∆t, t− 3∆t and t− 4∆t, voltage attime t− ∆t, t− 2∆t, t− 3∆t and t− 4∆t for the first sampling scenario in the case oftwo series motor in parallel with a resistive and an inductive Load.selected sample i(t) i(t− ∆t) i(t− 2∆t) i(t− 3∆t) i(t− 4∆t) v(t− ∆t) v(t− 2∆t) v(t− 3∆t) v(t− 4∆t)20th sample -1.08422 -1.18047 -1.20521 -1.1643 -1.06742 2.11873 2.86E+00 3.356094 3.59304830th sample 1.847889 1.670556 1.409224 1.084638 0.719978 -6.90432 -6.97878 -6.65809 -5.9884140th sample -1.63121 -1.07295 -0.48892 0.086761 0.622167 8.778512 7.294637 5.486276 3.47267550th sample -0.58222 -1.30023 -1.92157 -2.41571 -2.76055 -3.84574 -0.95609 1.865721 4.46486860th sample 3.152638 3.541393 3.719616 3.684601 3.445646 -6.82702 -9.59789 -11.7463 -13.174570th sample -4.04172 -3.75978 -3.27442 -2.6127 -1.81158 14.79497 15.50839 15.36482 14.3721980th sample 2.216767 1.375249 0.45409 -0.49584 -1.42207 -13.8608 -11.8291 -9.14379 -5.9531690th sample 0.68201 1.587901 2.402648 3.081205 3.586052 4.807176 1.220537 -2.43355 -5.95316100th sample -3.39011 -3.83106 -4.0635 -4.07461 -3.8638 7.062435 10.10304 12.58535 14.37219Table 4.50: Values of voltage at time t for the first sampling scenario in the case oftwo series motor in parallel with a resistive and an inductive Load.selected sample v(t)20th sample -1.1499930th sample 6.41091840th sample -9.8322450th sample 6.63800660th sample 3.57103270th sample -13.26480th sample 15.1265390th sample -8.12817100th sample -3.63156Discrete time series coefficients are calculated and resulted as:A=[5.4039 −17.4479 20.7260 −10.7033 2.0216 −3.1295 3.7334 −2.0526 0.4507]109Simulation Six: Second Sampling Scenario For The Fourth Order StatusIn this sampling scheme, the first sampling point is the 30th point of voltage and currentwaveforms. The distance between each samples is 5 samples.Table 4.51: Values of current at time t, t− ∆t, t− 2∆t, t− 3∆t and t− 4∆t, voltage attime t− ∆t, t− 2∆t, t− 3∆t and t− 4∆t for the second sampling scenario in the caseof two series motor in parallel with a resistive and an inductive Load.selected sample i(t) i(t− ∆t) i(t− 2∆t) i(t− 3∆t) i(t− 4∆t) v(t− ∆t) v(t− 2∆t) v(t− 3∆t) v(t− 4∆t)30th sample 1.847889 1.670556 1.409224 1.084638 0.719978 -6.90432 -6.97878 -6.65809 -5.9884135th sample 1.08932 1.46565 1.73506 1.888555 1.92441 -0.6713 -2.56383 -4.19955 -5.540th sample -1.63121 -1.07295 -0.48892 0.086761 0.622167 8.778512 7.294637 5.486276 3.47267545th sample -2.94331 -2.96092 -2.8198 -2.53509 -2.1294 8.477551 9.702566 10.33646 10.370950th sample -0.58222 -1.30023 -1.92157 -2.41571 -2.76055 -3.84574 -0.95609 1.865721 4.46486855th sample 3.023137 2.447016 1.754769 0.989035 0.195007 -13.6991 -12.8223 -11.2744 -9.1666760th sample 3.152638 3.541393 3.719616 3.684601 3.445646 -6.82702 -9.59789 -11.7463 -13.174565th sample -0.91593 0.023743 0.953767 1.819716 2.568313 10.10304 7.062435 3.631558 -4.10E-1370th sample -4.04172 -3.75978 -3.27442 -2.6127 -1.81158 14.79497 15.50839 15.36482 14.37219Table 4.52: Values of voltage at time t for the second sampling scenario in the case oftwo series motor in parallel with a resistive and an inductive Load.selected sample v(t)30th sample 6.41091835th sample -1.3790140th sample -9.8322445th sample -6.7054550th sample 6.63800655th sample 13.8283460th sample 3.57103265th sample -12.585470th sample -13.264Discrete time series coefficients are calculated and resulted as:B=[5.2638 −16.9958 20.1887 −10.4259 1.9692 −3.1492 3.7865 −2.1002 0.4650]Resulted coefficients matrices as for order =4 are equivalent which confirms that this load110is a 4th order load.Simulation Six: First Sampling Scenario For The Fifth Order StatusIn this sampling scheme, the first sampling point is the 20th point of voltage and currentwaveforms. The distance between each samples is 5 samples. Discrete time series coeffi-cients are calculated and resulted as:[a0, a1, a2, a3, a4, a5,b0,b1,b2,b3,b4] =[14.74 −31.096 3.26 34.092 −27.1717 6.173 −0.5759 −2.310 2.0379 0.2863 −0.4369]Simulation Six: Second Sampling Scenario For the Fifth Order StatusIn this sampling scheme, the first sampling point is the 10th point of voltage and currentwaveforms. The distance between each samples is 10 samples. Discrete time series coeffi-cients are calculated and resulted as:[a0, a1, a2, a3, a4, a5,b0,b1,b2,b3,b4] =[8.89 −26.292 26.25 −8.26 −1.502 0.912 −2.582 2.537 −1.470 0.710 −0.192]Coefficients values for 5th order mismatch. Therefore, the largest order which results inequal coefficients is four.1114.2.11 Simulation Seven: Three Series Motor in Parallel with a Resistive andan Inductive LoadFigure 4.26: Three series motor in parallel with resistive and inductive load.In this simulation case, sampling rate is∆t= 2200µs. {tmax= 0.7seconds and Rampuptime=0.05 (slow transient ends at sample number 22). This simulation describes two series mo-tor loads in parallel with a resistance and an inductance.Simulation Seven: First Sampling Scenario For the Fifth Order StatusIn this sampling scheme, the first sampling point is the 10th point of voltage and currentwaveforms. The distance between each samples is 5 samples.112Table 4.53: Values of current at time t, t − ∆t, t − 2∆t, t − 3∆t, t − 4∆t and t − 5∆t,voltage at time t − ∆t , t − 2∆t, t − 3∆t, t − 4∆t and t − 5∆t for the first samplingscenario in the case of three series motors in parallel with a resistive and an inductiveloadselected sample i(t) i(t− ∆t) i(t− 2∆t) i(t− 3∆t) i(t− 4∆t) i(t− 5∆t) v(t− ∆t) v(t− 2∆t) v(t− 3∆t) v(t− 4∆t) v(t− 5∆t)10th sample 1.193384 1.112585 0.358973 -0.45581 -0.82461 -0.65301 -5.14035 -4.25538 -1.07124 1.833812 2.69565615th sample -0.31091 -1.69494 -1.88942 -0.94974 0.368754 1.193384 1.886488 7.096053 7.189436 2.914371 -2.3395620th sample -2.31235 -0.18528 1.803115 2.431806 1.463945 -0.31091 8.766567 -0.43857 -8.35002 -10.1862 -5.5347125th sample 3.403337 3.010995 0.694947 -1.90806 -3.11962 -2.31235 -15.1601 -12.4012 -1.98061 8.726071 12.9886830th sample -0.95577 -3.11054 -3.2285 -1.22808 1.58982 3.403337 5.177457 14.31545 14.15793 4.807176 -7.6650435th sample -2.29308 0.39524 2.838765 3.451327 1.836957 -0.95577 9.61164 -2.53004 -13.0289 -15.0676 -7.3224340th sample 3.46756 2.74645 0.250337 -2.3994 -3.48173 -2.29308 -15.4778 -11.6041 -0.19548 11.34009 15.5121545th sample -1.38373 -3.29701 -3.0628 -0.83291 1.945039 3.46756 6.975206 14.96564 13.23837 2.914969 -9.3012250th sample -1.95381 0.819001 3.065566 3.327291 1.434418 -1.38373 8.002801 -4.43386 -13.9915 -14.4639 -5.5444755th sample 3.503459 2.446064 -0.19479 -2.70418 -3.45255 -1.95381 -15.5514 -10.2141 1.755633 12.58535 15.2429760th sample -1.77804 -3.41712 -2.83301 -0.40488 2.290679 3.503459 8.662952 15.37981 12.11003 0.976791 -10.7907Table 4.54: Values of voltage at time t for the first sampling scenario in the case of athree series motors in parallel with a resistive and an inductive loadselected sample v(t)10th sample 2.33955915th sample 5.53471220th sample -12.988725th sample 7.66503930th sample 7.32243535th sample -15.512240th sample 9.30121945th sample 5.54446850th sample -15.24355th sample 10.7907160th sample 3.679061Discrete time series coefficients are calculated and resulted as:[a0, a1, a2, a3, a4, a5,b0,b1,b2,b3,b4] =A=[2.261 −4.796 1.66 2.1376 −1.330 0.067 −1.918 1.553 −0.484 0.050 0.1637]113Simulation Seven: Second Sampling Scenario For The Fifth Order StatusIn this sampling scheme, the first sampling point is the 15th point of voltage and currentwaveforms. The distance between each samples is 5 samples.Table 4.55: Values of current at time t, t − ∆t, t − 2∆t, t − 3∆t, t − 4∆t and t − 5∆t,voltage at time t− ∆t , t− 2∆t, t− 3∆t, t− 4∆t and t− 5∆t for the second samplingscenario in the case of three series motors in parallel with a resistive and an inductiveloadselected sample i(t) i(t− ∆t) i(t− 2∆t) i(t− 3∆t) i(t− 4∆t) i(t− 5∆t) v(t− ∆t) v(t− 2∆t) v(t− 3∆t) v(t− 4∆t) v(t− 5∆t)15th sample -0.31091 -1.69494 -1.88942 -0.94974 0.368754 1.193384 1.886488 7.096053 7.189436 2.914371 -2.3395620th sample -2.31235 -0.18528 1.803115 2.431806 1.463945 -0.31091 8.766567 -0.43857 -8.35002 -10.1862 -5.5347125th sample 3.403337 3.010995 0.694947 -1.90806 -3.11962 -2.31235 -15.1601 -12.4012 -1.98061 8.726071 12.9886830th sample -0.95577 -3.11054 -3.2285 -1.22808 1.58982 3.403337 5.177457 14.31545 14.15793 4.807176 -7.6650435th sample -2.29308 0.39524 2.838765 3.451327 1.836957 -0.95577 9.61164 -2.53004 -13.0289 -15.0676 -7.3224340th sample 3.46756 2.74645 0.250337 -2.3994 -3.48173 -2.29308 -15.4778 -11.6041 -0.19548 11.34009 15.5121545th sample -1.38373 -3.29701 -3.0628 -0.83291 1.945039 3.46756 6.975206 14.96564 13.23837 2.914969 -9.3012250th sample -1.95381 0.819001 3.065566 3.327291 1.434418 -1.38373 8.002801 -4.43386 -13.9915 -14.4639 -5.5444755th sample 3.503459 2.446064 -0.19479 -2.70418 -3.45255 -1.95381 -15.5514 -10.2141 1.755633 12.58535 15.2429760th sample -1.77804 -3.41712 -2.83301 -0.40488 2.290679 3.503459 8.662952 15.37981 12.11003 0.976791 -10.790765th sample -1.57773 1.236617 3.251941 3.159695 1.019851 -1.77804 6.267753 -6.26775 -14.7334 -13.6321 -3.67906Table 4.56: Values of voltage at time t for the second sampling scenario in the case ofa three series motors in parallel with a resistive and an inductive load.selected sample v(t)15th sample 5.53471220th sample -12.988725th sample 7.66503930th sample 7.32243535th sample -15.512240th sample 9.30121945th sample 5.54446850th sample -15.24355th sample 10.7907160th sample 3.67906165th sample -14.7334114Discrete time series coefficients are calculated and resulted as:[a0, a1, a2, a3, a4, a5,b0,b1,b2,b3,b4] =B=[2.2236 −4.714 1.632 2.1014 −1.3074 0.066 −1.9253 1.576 −0.4994 0.050 0.1667]Two derived impedance coefficient matrices are equivalent for the 5th order analysis.Also, 6th order calculations yield to two different impedance matrices. Thus, this loadis a 5th order system.4.2.12 Simulation Eight: Two Parallel Motors LoadsFigure 4.27: Two parallel motors loads.This example illustrates an aggregated load constituting of two parallel motor loads.In this simulation case, sampling rate is ∆t = 625µs(sampling frequency =1600 Hz), max-imum simulation time is tmax = 0.0625s and ramp up time=0.03(slow transient ends atsample number 48).Simulation Eight: First Sampling Scenario For The Third Order StatusIn this sampling scheme, the first sampling point is the 30th point of voltage and currentwaveforms. The distance between each samples is 4 samples.115Table 4.57: Values of current at time t, t − ∆t, t − 2∆t and t − 3∆t, voltage at timet− ∆t, t− 2∆t and t− 3∆t for the first sampling scenario in the case of two parallelmotors.selected sample i(t) i(t− ∆t) i(t− 2∆t) i(t− 3∆t) v(t− ∆t) v(t− 2∆t) v(t− 3∆t)30th sample 10.61917 9.017144 7.040091 4.826744 -8.6304 -8.72347 -8.3226134th sample 11.15864 12.0443 12.22007 11.72101 -3.20478 -5.24944 -6.87538th sample 1.289478 4.522545 7.326548 9.56945 6.857844 4.340844 1.72376442th sample -12.313 -9.21498 -5.77706 -2.20721 12.96362 12.2903 10.9731446th sample -17.7378 -17.7059 -16.7179 -14.8728 8.38181 10.59694 12.1282150th sample -8.25894 -11.9136 -14.8025 -16.7668 -4.80718 -1.19511 2.33215154th sample 9.265115 4.953669 0.424436 -4.06529 -14.795 -13.264 -11Table 4.58: Values of voltage at time t for the first sampling scenario in the case oftwo parallel motors.selected sample v(t)30th sample 8.01364734th sample 0.83911938th sample -9.118342th sample -12.920646th sample -5.5810850th sample 8.1281754th sample 15.50839Discrete time series coefficients are calculated and resulted as:A=[5.0677 −11.4741 8.1121 −1.7024 −0.5052 −0.9970 0.51539]Simulation Eight: Second Sampling Scenario For The Third Order StatusIn this sampling scheme, the first sampling point is the 40th point of voltage and currentwaveforms. The distance between each samples is 4 samples.116Table 4.59: Values of current at time t, t − ∆t, t − 2∆t and t − 3∆t, voltage at timet−∆t, t− 2∆t and t− 3∆t for the second sampling scenario in the case of two parallelmotors.selected sample i(t) i(t− ∆t) i(t− 2∆t) i(t− 3∆t) v(t− ∆t) v(t− 2∆t) v(t− 3∆t)40th sample -5.77706 -2.20721 1.289478 4.522545 10.97314 9.118296 6.85784444th sample -16.7179 -14.8728 -12.313 -9.21498 12.12821 12.92057 12.9636248th sample -14.8025 -16.7668 -17.7378 -17.7059 2.332151 5.581085 8.3818152th sample 0.424436 -4.06529 -8.25894 -11.9136 -11 -8.12817 -4.8071856th sample 16.28371 13.11448 9.265115 4.953669 -15.3648 -15.5084 -14.79560th sample 19.32273 20.15932 19.91034 18.59294 -7.06243 -10.103 -12.585464th sample 6.816886 11.01555 14.62358 17.44378 7.062435 3.631558 -4.10E-13Table 4.60: Values of voltage at time t for the second sampling scenario in the case oftwo parallel motors.selected sample v(t)40th sample -12.290344th sample -10.596948th sample 1.19510952th sample 13.2639756th sample 14.3721960th sample 3.63155864th sample -10.103Discrete time series coefficients are calculated and resulted as:B=[5.0644 −11.4930 8.1528 −1.7211 −0.5136 −0.9904 0.5159]Impedance coefficient values for t = t, t − ∆t, t − 2∆t, t − 3∆t are equal for the two sam-pling scenarios. Therefore, we will proceed with order=4 coefficient calculation. Threedifferent sampling sets are selected.Simulation Eight: First Sampling Scenario For The Fourth Order StatusIn this sampling scheme, the first sampling point is the 20th point of voltage and currentwaveforms. The distance between each samples is 5 samples.117Table 4.61: Values of current at time t, t−∆t, t− 2∆t and t− 3∆t and t− 4∆t, voltageat time t− ∆t, t− 2∆t, t− 3∆t and t− 4∆t for the first sampling scenario in the caseof two parallel motors.selected sample i(t) i(t− ∆t) i(t− 2∆t) i(t− 3∆t) i(t− 4∆t) v(t− ∆t) v(t− 2∆t) v(t− 3∆t) v(t− 4∆t)20th sample -6.79403 -7.00161 -6.82318 -6.31469 -5.54857 2.648413 3.57816 4.195117 4.4913125th sample 0.258391 -1.83327 -3.64976 -5.11123 -6.16645 -4.84104 -3.23695 -1.58881 1.74E-1330th sample 10.61917 9.017144 7.040091 4.826744 2.52011 -8.6304 -8.72347 -8.32261 -7.4855235th sample 9.56945 11.15864 12.0443 12.22007 11.72101 -0.83912 -3.20478 -5.24944 -6.87540th sample -5.77706 -2.20721 1.289478 4.522545 7.326548 10.97314 9.118296 6.857844 4.34084445th sample -17.7059 -16.7179 -14.8728 -12.313 -9.21498 10.59694 12.12821 12.92057 12.9636250th sample -8.25894 -11.9136 -14.8025 -16.7668 -17.7378 -4.80718 -1.19511 2.332151 5.58108555th sample 13.11448 9.265115 4.953669 0.424436 -4.06529 -15.5084 -14.795 -13.264 -1160th sample 19.32273 20.15932 19.91034 18.59294 16.28371 -7.06243 -10.103 -12.5854 -14.3722Table 4.62: Values of voltage at time t for the first sampling scenario in the case oftwo parallel motors.selected sample v(t)20th sample -1.4374925th sample 6.29267530th sample 8.01364735th sample -1.7237640th sample -12.290345th sample -8.3818150th sample 8.1281755th sample 15.3648260th sample 3.631558Discrete time series coefficients are calculated and resulted as:[a0, a1, a2, a3, a4,b0,b1,b2,b3] =[5.0678 −16.4878 19.4712 −9.7424 1.6912 −1.4958 −0.4912 1.4957 −0.5087]Simulation Eight: Second Sampling Scenario For The Fourth Order StatusIn this sampling scheme, the first sampling point is the 30th point of voltage and currentwaveforms. The distance between each samples is 5 samples.118Table 4.63: Values of current at time t, t−∆t, t− 2∆t and t− 3∆t and t− 4∆t, voltageat time t− ∆t, t− 2∆t, t− 3∆t and t− 4∆t for the second sampling scenario in thecase of two parallel motors.selected sample i(t) i(t− ∆t) i(t− 2∆t) i(t− 3∆t) i(t− 4∆t) v(t− ∆t) v(t− 2∆t) v(t− 3∆t) v(t− 4∆t)30th sample 10.61917 9.017144 7.040091 4.826744 2.52011 -8.6304 -8.72347 -8.32261 -7.4855235th sample 9.56945 11.15864 12.0443 12.22007 11.72101 -0.83912 -3.20478 -5.24944 -6.87540th sample -5.77706 -2.20721 1.289478 4.522545 7.326548 10.97314 9.118296 6.857844 4.34084445th sample -17.7059 -16.7179 -14.8728 -12.313 -9.21498 10.59694 12.12821 12.92057 12.9636250th sample -8.25894 -11.9136 -14.8025 -16.7668 -17.7378 -4.80718 -1.19511 2.332151 5.58108555th sample 13.11448 9.265115 4.953669 0.424436 -4.06529 -15.5084 -14.795 -13.264 -1160th sample 19.32273 20.15932 19.91034 18.59294 16.28371 -7.06243 -10.103 -12.5854 -14.372265th sample 2.257668 6.816886 11.01555 14.62358 17.44378 10.10304 7.062435 3.631558 -4.10E-1370th sample -17.2642 -14.5713 -11.0641 -6.93537 -2.41191 14.79497 15.50839 15.36482 14.37219Table 4.64: Values of voltage at time t for the second sampling scenario in the case oftwo parallel motors.selected sample v(t)30th sample 8.01364735th sample -1.7237640th sample -12.290345th sample -8.3818150th sample 8.1281755th sample 15.3648260th sample 3.63155865th sample -12.585470th sample -13.264Discrete time series coefficients are calculated and resulted as:[a0, a1, a2, a3, a4,b0,b1,b2,b3] =[5.0678 −16.4825 19.4589 −9.7335 1.6892 −1.4947 −0.4918 1.4947 −0.5081]Simulation Eight: Third Sampling Scenario For The Fourth Order StatusIn this sampling scheme, the first sampling point is the 40th point of voltage and currentwaveforms. The distance between each samples is 5 samples.119Table 4.65: Values of current at time t, t−∆t, t− 2∆t and t− 3∆t and t− 4∆t, voltageat time t− ∆t, t− 2∆t, t− 3∆t and t− 4∆t for the third sampling scenario in the caseof two parallel motors.selected sample i(t) i(t− ∆t) i(t− 2∆t) i(t− 3∆t) i(t− 4∆t) v(t− ∆t) v(t− 2∆t) v(t− 3∆t) v(t− 4∆t)40th sample -5.77706 -2.20721 1.289478 4.522545 7.326548 10.97314 9.118296 6.857844 4.34084445th sample -17.7059 -16.7179 -14.8728 -12.313 -9.21498 10.59694 12.12821 12.92057 12.9636250th sample -8.25894 -11.9136 -14.8025 -16.7668 -17.7378 -4.80718 -1.19511 2.332151 5.58108555th sample 13.11448 9.265115 4.953669 0.424436 -4.06529 -15.5084 -14.795 -13.264 -1160th sample 19.32273 20.15932 19.91034 18.59294 16.28371 -7.06243 -10.103 -12.5854 -14.372265th sample 2.257668 6.816886 11.01555 14.62358 17.44378 10.10304 7.062435 3.631558 -4.10E-1370th sample -17.2642 -14.5713 -11.0641 -6.93537 -2.41191 14.79497 15.50839 15.36482 14.3721975th sample -15.2861 -17.7612 -19.2498 -19.6691 -18.9951 1.220537 4.807176 8.12817 1180th sample 5.666001 1.092029 -3.53966 -7.97276 -11.9619 -13.8608 -11.8291 -9.14379 -5.95316Table 4.66: Values of voltage at time t for the third sampling scenario in the case oftwo parallel motors.selected sample v(t)40th sample -12.290345th sample -8.3818150th sample 8.1281755th sample 15.3648260th sample 3.63155865th sample -12.585470th sample -13.26475th sample 2.43354980th sample 15.12653Discrete time series coefficients are calculated and resulted as:[a0, a1, a2, a3, a4,b0,b1,b2,b3] =[5.0678 −16.4670 19.4235 −9.7082 1.6839 −1.4916 −0.4937 1.4919 −0.5066]Fifth order impedance coefficients are not equivalent for two selected sampling patterns.The largest order which results in the equal impedance coefficients is 4. Usually, in thecase of parallel motors, if they are not identical, each single motor will increase the totalorder of the circuit, by two. In other words, if n motor loads are connected to a feeder, the120aggregated load will be a 2n order system.Figure 4.28: Third order impedance coefficients for two different sampling scenariosof two parallel motor loads.121Figure 4.29: Fourth order impedance coefficients for three different sampling scenar-ios of two parallel motor loads.1224.2.13 Simulation Eight: Two Parallel Motors Load Parameters CalculationsFigure 4.30: Two parallel motors.Ys in Figure 4.31 is derived from the 4th order calculated discrete impedance coefficients.Following division products 8 steps of division prior to getting a remainder of zeroThe first eight quotients are indicating two motor’s parameters respectively, R1,Lm,R2&L2& R3,L2m,R4&L4.123Figure 4.31: Two parallel motors load parameter calculation.1244.2.14 Simulation Nine: Two Parallel Motors in Parallel with One ResistiveLoadFigure 4.32: Two parallel motor in parallel with one resistive load.In this simulation case, sampling rate is ∆t = 625µs(sampling frequency =1600 Hz), max-imum simulation time is tmax = 0.0625s and ramp up time=0.03 (slow transient ends atsample number 48). For this case, since resistor does not change the total order of theaggregated circuit, this load is a 4th order load same as Section 4.2.12.Simulation Nine: First Sampling Scenario For The Fourth Order StatusIn this sampling scheme, the first sampling point is the 30th point of voltage and currentwaveforms. The distance between each samples is 5 samples.125Table 4.67: Values of current at time t, t− ∆t, t− 2∆t, t− 3∆t and t− 4∆t, voltage attime t− ∆t, t− 2∆t, t− 3∆t and t− 4∆t for the first sampling scenario in the case oftwo parallel motors in parallel with a resistive load.selected sample i(t) i(t− ∆t) i(t− 2∆t) i(t− 3∆t) i(t− 4∆t) v(t− ∆t) v(t− 2∆t) v(t− 3∆t) v(t− 4∆t)30th sample 10.77944 9.189752 7.21456 4.993196 2.67E+00 -8.6304 -8.72347 -8.32261 -7.4855235th sample 9.534974 11.17542 12.1084 12.32506 11.85851 -0.83912 -3.20478 -5.24944 -6.87540th sample -6.02286 -2.42667 1.107112 4.385388 7.239731 10.97314 9.118296 6.857844 4.34084445th sample -17.8735 -16.9299 -15.1153 -12.5714 -9.47426 10.59694 12.12821 12.92057 12.9636250th sample -8.09637 -11.8174 -14.7786 -16.8134 -17.8494 -4.80718 -1.19511 2.332151 5.58108555th sample 13.42177 9.575283 5.249569 0.689716 -3.84529 -15.5084 -14.795 -13.264 -1160th sample 19.39536 20.30056 20.1124 18.84465 16.57115 -7.06243 -10.103 -12.5854 -14.372265th sample 2.005961 6.614825 10.8743 14.55095 17.44378 10.10304 7.062435 3.631558 -4.10E-1370th sample -17.5295 -14.8672 -11.3743 -7.24266 -2.69935 14.79497 15.50839 15.36482 14.37219Table 4.68: Values of voltage at time t for the second sampling scenario in the case oftwo parallel motors in parallel with a resistive load.selected sample v(t)30th sample 8.01364735th sample -1.7237640th sample -12.290345th sample -8.3818150th sample 8.1281755th sample 15.3648260th sample 3.63155865th sample -12.585470th sample -13.264Discrete time series coefficients are calculated and resulted as:B=[4.6014 −14.9632 17.6626 −8.8337 1.5329 −1.6559 −0.0936 1.1800 −0.4305]Simulation Nine: Second Sampling Scenario For The Fourth Order StatusIn this sampling scheme, the first sampling point is the 40th point of voltage and currentwaveforms. The distance between each samples is 5 samples.126Table 4.69: Values of current at time t, t− ∆t, t− 2∆t, t− 3∆t and t− 4∆t, voltage attime t− ∆t, t− 2∆t, t− 3∆t and t− 4∆t for the second sampling scenario in the caseof two parallel motors in parallel with a resistive load.selected sample i(t) i(t− ∆t) i(t− 2∆t) i(t− 3∆t) i(t− 4∆t) v(t− ∆t) v(t− 2∆t) v(t− 3∆t) v(t− 4∆t)40th sample -6.02286 -2.42667 1.107112 4.385388 7.239731 10.97314 9.12E+00 6.857844 4.34084445th sample -17.8735 -16.9299 -15.1153 -12.5714 -9.47426 10.59694 12.12821 12.92057 12.9636250th sample -8.09637 -11.8174 -14.7786 -16.8134 -17.8494 -4.80718 -1.19511 2.332151 5.58108555th sample 13.42177 9.575283 5.249569 0.689716 -3.84529 -15.5084 -14.795 -13.264 -1160th sample 19.39536 20.30056 20.1124 18.84465 16.57115 -7.06243 -10.103 -12.5854 -14.372265th sample 2.005961 6.614825 10.8743 14.55095 17.44378 10.10304 7.062435 3.631558 -4.10E-1370th sample -17.5295 -14.8672 -11.3743 -7.24266 -2.69935 14.79497 15.50839 15.36482 14.3721975th sample -15.2374 -17.7856 -19.3459 -19.8316 -19.2151 1.220537 4.807176 8.12817 1180th sample 5.968531 1.369245 -3.30308 -7.78988 -11.8428 -13.8608 -11.8291 -9.14379 -5.95316Table 4.70: Values of voltage at time t for the second sampling scenario in the case oftwo parallel motors in parallel with a resistive load.selected sample v(t)40th sample -12.290345th sample -8.3818150th sample 8.1281755th sample 15.3648260th sample 3.63155865th sample -12.585470th sample -13.26475th sample 2.43354980th sample 15.12653Discrete time series coefficients are calculated and resulted as:B=[4.6014 −14.9667 17.6706 −8.8395 1.5341 −1.6567 −0.0930 1.1806 −0.4308]127Figure 4.33: Third order impedance coefficients for two different sampling scenariosof two parallel motors in parallel with one resistance load.Figure 4.34: Fourth order impedance coefficients for two different sampling scenariosof two parallel motors in parallel with one resistance load.128Figure 4.35: Fifth order impedance coefficients for two different sampling scenariosof two parallel motors in parallel with one resistance load.4.2.15 Simulation Nine: Two Parallel Motors in Parallel with A ResistiveLoad Parameters CalculationsFigure 4.36: Two parallel motors in parallel with a resistive load.129In Figure 4.37, Ys is derived from the 4th order impedance coefficients. In this case, we see9 stages of division before we reach zero in remainder. The first eight quotients are thevalues of the 2 motors parameters respectively, R1,Lm,R2&L2 & R3,L2m,R4&L4 and thelast quotient is the value of paralleled resistance.130Figure 4.37: Two parallel motors load parameter calculation.1314.2.16 Simulation Ten: Two Parallel Motors in Parallel with A Resistive andAn Inductive Load Parameters CalculationsIn this section, we describe two parallel motors in parallel with a resistive and inductiveload. The load is a 5th order load. Since inductance generates 1st order history terms, itincrease the total order of the aggregated load circuit by one. According to our algorithm,the order of an aggregated load combining of motor loads and inductive loads in parallelis formulated as:(number of inductive loads)*1+(number of motor loads *2)Figure 4.38: Two Parallel motors in parallel with a resistive and an inductive load.This case determines an inductor type in parallel with two induction motors. In Fig-ure 4.39, Ys is derived from the 5th order discrete time impedance coefficients. This loadscenario products ten steps of long division before reaching to zero in remainder.132Figure 4.39: Two parallel motors in parallel with a series RL load parameter calcula-tion. 133The first eight quotients are 2 motors parameters respectively namely, R1,Lm,R2&L2& R3,L2m,R4&L4. The 9th and the 10th quotients are the values of paralleled Resistorand inductor within some errors. In conclusion, having the electrical parameters of anaggregated load, and the types of its constituent eigenloads, will determine, what type ofloads are connected to the feeder at each time stamp. Figure 4.40 demonstrates the generaloverview of our load disaggregation algorithm.Figure 4.40: General overview of or work.In which u1,u2 and u3 are the contributions of each eigen-loads in the total load.134Chapter 5Experimental ResultsThe proposed “Load decomposition” algorithm is validated with real experimental dataaccomplished on a number of different appliance types at “UBC Power Lab”. As dis-cussed in Chapter 3, voltage and current signals in near-real time domain, are the inputsof the proposed algorithm. The main goal of this chapter is to verify the proposed methodwith real measurements. A number of electric appliances were analysed in the lab andcurrent and voltage waveforms were recorded for a period of few seconds. Samplingfrequency is 10KHz. A current and voltage probe (different probe) were attached to anoscilloscope for measuring the voltage and current waveforms. Table 5.1 summarizes thelist of appliances which have been tested.5.1 Appliances Electrical Load Types and ClassificationIn this section, we briefly explain about four major category of appliances according totheir characteristics.5.1.1 Resistive LoadsLoads consisting of any type of heating elements such as conventional electric heating ap-pliances (without power electronic interfaces), iron, heater, toaster, oven, kettles, electricrice cookers, electric hair curlers and even incandescent lights are categorized as resistiveloads. These type of loads have negligible reactive power.135Table 5.1: Experimental appliances list.Appliance type order of the load/type of the loadCompact Florescent Lamp with magnetic ballast 1st order-inductiveFan 2nd order-Motor loadRefrigerator 2nd order-Motor LoadHair Dryer 2nd order-Motor LoadLCD Samsung power electronicsiphone charger power electronicsHTC charger power electronicsAsus charger power electronicsRazor 2nd order and MotorVacuum Cleaner 2nd order and Motor typeIron resistive loadincandescent lamp purely Resistive5.1.2 Inductive LoadsAn example of this group are florescent lamps which convert the electromagnetic radia-tion to the light.5.1.3 Motor LoadsInduction motors have the most appearance among the household appliances. Fan, vac-uum cleaner, dishwasher and refrigerator are examples of motor loads.5.1.4 Non-Linear LoadsThese loads do not draw the sinusoidal current pattern. The most predominant ones arepower electronic loads such as charger, TV and desktop computer. Figure 5.1 shows themeasured voltage at the outlet. To obtain the aggregated voltage value, we measuredthe voltage between the outlet and appliances. Note that, in order to mimic the similarsituation as in real power gird, all loads are directly connected to the 120 volt outlet. Weuse the same voltage signal for all the appliances.136Figure 5.1: Voltage measured at main outlet.5.2 Measuring General ObservationsSome loads such as the vacuum cleaners exhibit a fairly constant frequency over the time.While the magnitude of current has a decaying component (inrush current). Some loadshave multiple stages from the starting moment to the steady state. The current waveformof a HTC laptop charger as shown in Figure 5.5, has this characteristic. Current has thehigher frequency components at the middle stage and then starts damping after a fewcycles.5.2.1 FanFigure 5.2 indicates the current waveform for a table fan captured on-line while it wasconnected to the city electricity.137Figure 5.2: Current waveform in time for a fan.5.2.2 Fan Parameters CalculationAs the first step, we need to choose at least two different sampling scenarios or two differ-ent harmonic data sets. We extracted the 60 Hz instantaneous voltage and current signals.Additionally, two sampling patterns were selected due to calculating the impedance co-efficients. Both vt and it are sampled with sampling rate of ∆t = 1.6e-06 seconds. Asdiscussed in Chapter 3, from Equation 3.4, we derive the second order impedance coeffi-cients for two different sampling scenarios as following.First Sampling Scenario For the Second Order StatusA = [a0, a1, a2,b0,b1] =[257.1647 −499.2356 242.3976 −1.9700 0.9726]Second Sampling Scenario For the Second Order StatusB = [a0, a1, a2,b0,b1] =[259.2187 −499.2816 240.2983 −1.9631 0.9653]Matrices A and B are equivalent within 10% error. We also, analysed the third ordercoefficients. The derived values were not equivalent for two different sampling schemes.Implementing second Cauer synthesis method on the continuous admittance function ofthe fan leads to:R1 = 107ohmLm = 11.02mH138R2 = 419.43ohmL2 = 11.02mH5.2.3 Florescent LampFigure 5.3 demonstrates the current waveform of a magnetic ballast florescent lamp. Thereare some noise and also harmonics embedded in the current waveform. We extracted60Hz component and implemented our method on the instantaneous voltage and currentdata. Florescent lamp was confirmed to be a first order load based on real voltage andcurrent data.Figure 5.3: Current waveform in time for a florescent lamp5.2.4 Florescent Parameters CalculationAs the first step, we need to choose at least two different sampling scenarios or two dif-ferent harmonic data sets. We extracted the 60 Hz instantaneous voltage and current sig-nals. Additionally, two sampling patterns were selected due to calculating the impedancecoefficients. Both vt and it are sampled with sampling rate of ∆t = 8e-07 seconds. As dis-cussed by Equation 3.4, we derive the first order impedance coefficients for two differentsampling scenarios as following.First Sampling Scenario For the First Order StatusA = [a0, a1,b0] =[1567.7 −1567.8 −1]139Second Sampling Scenario For the First Order StatusB = [a0, a1, a2,b0,b1] =[1567.5 −1567.5 −1]Matrices A and B are equivalent within 10% error. We also, analysed the second ordercoefficients. The derived values were not equivalent for two different sampling schemes.Implementing second Cauer synthesis method on the continuous admittance function ofthe fan leads to:L = 0.637mH5.2.5 Fan in Parallel With Florescent Lamp Parameters CalculationWe put the fan and the florescent lamp in parallel for this case and aggregated voltageand current are measured.5.2.6 Parallel Florescent and Fan Parameters CalculationAs the first step, we need to choose at least two different sampling scenarios or two differ-ent harmonic data sets. We extracted the 60 Hz instantaneous voltage and current signals.Additionally, two sampling patterns were selected due to calculating the impedance co-efficients. Both vt and it are sampled with sampling rate of ∆t = 1.6e-06 seconds. Asdiscussed in Chapter 3, from Equation 3.4, we derive the second order impedance coeffi-cients for two different sampling scenarios as following.First Sampling Scenario For the Second Order StatusA = [a0, a1, a2,b0,b1] =[351.25 −703.5 −352 −2 1]Second Sampling Scenario For the Second Order StatusB = [a0, a1, a2,b0,b1] =[347.2160 −694.2740 347.0588 −1.9993 .9993]Matrices A and B are equivalent within 10% error. This confirms that this parallel loadtopology, at least satisfies the second order circuit requirements. Therefore the next pro-posed order is order number three. As shown by Equation 3.4, we derive the third orderimpedance coefficients for two different sampling scenarios as following.First Sampling Scenario For the Third Order Status Since the selected order is three, weneed at least seven coefficients. A = [a0, a1, a2, a3,b0,b1,b2] =[365 −1094.3 1093.6 −364.3 3 −3 1]140Second Sampling Scenario For the Third Order Status B = [a0, a1, a2, a3,b0,b1,b2] =[368.5 −1104.5 1103.6 −376.6 3 −3 1]Matrices A and B are equivalent within 10% error. This confirms that this parallel loadtopology, at least satisfies the third order circuit requirements. Therefore the next pro-posed order is order number four. Fourth order nine impedance coefficients including[a0, a1, a2, a3, a4,b0,b1,b2,b3] received different values for the two different group of 60 Hzvoltage and current samples. This proves that a parallel fan and florescent lamp load is athird order system. Note, as described in Chapter 3, motors are second order loads andinductive loads are first order loads. Whenever, we have an odd number for the aggre-gated load, there has to be some motor and some inductance. But, in the case of the evenorders, total load can either be a combination of motors, or parallel motor and resistances.We concluded that a parallel fan and florescent lamp is a third order load. This impliesthat there is a second order motor load and a first order inductive load in parallel. Sec-tions 5.2.1 and 5.2.3, already proved that a fan is a motor load and a florescent lamp is aninductive load type. Therefore, the parallel combination will be a third order load.5.2.7 Hair DryerFigure 5.4 shows the current waveform of a current waveform. The current signal is inphase with voltage signal which confirms that, hair dryer is a resistive load. We alsoverified this with our algorithm.141Figure 5.4: Current waveform in time for a hair-dryer.This thesis for the hair dryer, did not cover finding circuit parameters. However, basedour theory and simulation results in Section 3.4, we believe, hair dryer is a first order load.5.2.8 HTC ChargerHTC charger as discussed before is a power electronic load. For the power electronicloads, we extract at least two different frequency components from the voltage and thecurrent signals. Impedance coefficients are calculated for each set of two harmonics.142Figure 5.5: Current waveform in time for an HTC laptop charger.This thesis does not cover power electronic load analysis.5.2.9 RefrigeratorRefrigerator is an induction motor type of load. It is widely used for heavy duty applica-tions requiring high starting torques.Figure 5.6: Current waveform in time for a refrigerator.For the refrigerator, we did not cover finding circuit parameters. However, based ourtheory and simulation results in Section 3.4, we believe, refrigerator is a second order143load.5.3 Chapter SummaryIn this chapter, we validated our method with some real measurements, which have beenconducted on a few home electrical appliances. This proved that the method is imple-mentable to BC-Hydro distribution feeders.144Chapter 6Conclusions and Future Work6.1 SummaryIn this thesis, we devised an innovative solution for the load disaggregation challenges.So far, most researches in this domain, looked at load’s active and reactive power wave-forms to extract the distinctive patterns for individual load types. But this thesis exploitssmart meter voltage and current load information to discern between different type ofloads. Our algorithm is able to distinguish between motor, heating and lighting loads. Weproposed an EMTP based solution to solve the load disaggregation problem. The mostimportant breakthrough is that we are able to recognize the type and the value of the loadsat any level in distribution network. The result of this research is an input feeding the lin-ear power flow technique provided in paper [11]. LPF, assumes a linear voltage-powerdependency for all the load categories. Since, the dependency factor is different for eachtype of load, we need to know the load type before solving a LPF problem.First chapter, introduced the concept of load disaggregation and its applications in smartdistribution system. It also covered the definition of DMS, smart meter, and VVO. Chap-ter 2 summarized literature review. Most researchers employed machine learning tech-niques to classify different loads. Third chapter is the main body of this research. Itexplained, how to derive the discrete time values of the impedance matrix for each typeof load. Network synthesis is integrated into our algorithm to calculate load’s electricalcircuit parameters. In Chapter 4, we analysed some of PSCAD simulation load cases.Three major type of loads, (1) motor, (2) inductive, and (3) resistive are accounted in ourtechnique. Total loads are synthesized from different combinations of these eigen-loads.145For this reason, simulation examples covered different circuit topologies including R andL elements in the form of a motor, pure inductance and pure resistance. In Chapter 5, wevalidate our method with real measurement data.6.2 Future WorkThe following are suggested areas for future work and improvement to this thesis:1. This method is implementable on data received from BC-Hydro smart meters.2. This thesis can be coupled with Shifted Frequency Analysis (SFA) models and em-ployed by PSCAD. The shifted system is then numerically integrated to obtain dy-namic phasor solutions, which are more easily understood by power system oper-ators and planners than instantaneous time domain results. SFA allows the exactsimulation of frequencies within a band centred around a fundamental frequencyusing large time steps in a discrete-time EMTP type of solution environment [54].This powerful tool, gives us the freedom to choose the sampling rate from the loadvoltage and current signals.3. The other future’s improvement is to expand the selection of eigen-loads. One canconsider different category of induction motors, as an example. For instance, adopt-ing another circuit synthesis method, may result in more exact determination of loadcircuit parameters.4. The other extension can be the process of converting a discrete time impedancefunction to a continuous time impedance function. There is always a compromisebetween accuracy and stability in all the discrete-continuous mapping methods.5. One last limitation of this thesis is in its incapability of detecting individual induc-tance values and resistance components in total load. This work only reported theaggregated value for R and L. Working toward this challenge is another future work.146Bibliography[1] S.S. Venkata and H. Rudnick. Distribution systems [guest editorial]. Power andEnergy Magazine, IEEE, 5(4):16–22, July 2007. → pages 1[2] EDISON utlity. 2014 pv distribution system modeling workshop, May 2010. →pages xiv, 3[3] BC-Hydro. 35 kv and below interconnection requirements for power generators,May 2010. → pages 4[4] World energy needs and nuclear power, May 2015. → pages 5[5] key world energy statistics 2014, May 2014. → pages xiv, 5[6] H. Farhangi. The path of the smart grid. Power and Energy Magazine, IEEE,8(1):18–28, January 2010. → pages 6[7] Jixuan Zheng, D.W. Gao, and Li Lin. Smart meters in smart grid: An overview. InGreen Technologies Conference, 2013 IEEE, pages 57–64, April 2013. → pages 6[8] BC-Hydro. Epri technical executive discusses the importance of distributionmanagement, April 2012. → pages 7[9] Robert Uluski. Grid modernization, 2015. → pages xiv, 7[10] US department of energy report. Application of automated controls for voltage andreactive power management, December 2012. → pages 8[11] J.R. Marti, H. Ahmadi, and L. Bashualdo. Linear power-flow formulation based ona voltage-dependent load model. Power Delivery, IEEE Transactions on,28(3):1682–1690, July 2013. → pages 9, 10, 11, 55, 145147[12] A. Ellis, D. Kosterev, and A. Meklin. Dynamic load models: Where are we? InTransmission and Distribution Conference and Exhibition, 2005/2006 IEEE PES, pages1320–1324, May 2006. → pages 9[13] A. Dwyer, R.E. Nielsen, J. Stangl, and N.S. Markushevich. Load to voltagedependency tests at b.c. hydro. Power Systems, IEEE Transactions on, 10(2):709–715,May 1995. → pages 9[14] UBC. Local voltage stability assessment for variable load characteristics, May 2009.→ pages 11[15] Effects of reduced voltage on the operation and efficiency of electric loads), Sep1981. → pages 11[16] Non-intrusive load monitoring approaches for disaggregated energy sensing: Asurvey, 2012. → pages xiv, 12, 19, 20[17] Stephen Makonin. Nonintrusive load monitoring, December 2014. → pages 13[18] G.W. Hart. Nonintrusive appliance load monitoring. Proceedings of the IEEE,80(12):1870–1891, Dec 1992. → pages xiv, 13, 15[19] V. Amenta, S. Gagliano, G.M. Tina, G. Di Modica, and O. Tomarchio. Webinteractive non intrusive load disaggregation system for energy consumptionawareness. In AEIT Annual Conference, 2013, pages 1–6, Oct 2013. → pages xiv, 14[20] M. Zeifman and K. Roth. Nonintrusive appliance load monitoring: Review andoutlook. Consumer Electronics, IEEE Transactions on, 57(1):76–84, February 2011. →pages 15[21] A. Marchiori, D. Hakkarinen, Qi Han, and L. Earle. Circuit-level load monitoringfor household energy management. Pervasive Computing, IEEE, 10(1):40–48, Jan2011. → pages 15, 17[22] W.L. Chan, A.T.P. So, and L.L. Lai. Wavelet feature vectors for neural network basedharmonics load recognition. In Advances in Power System Control, Operation andManagement, 2000. APSCOM-00. 2000 International Conference on, volume 2, pages511–516 vol.2, Oct 2000. → pages 16, 18, 29148[23] At the flick of a switch: Detecting and classifying unique electrical events on theresidential power line (nominated for the best paper award), May 2010. → pages 16,17, 31[24] A. Cole and A. Albicki. Nonintrusive identification of electrical loads in athree-phase environment based on harmonic content. In Instrumentation andMeasurement Technology Conference, 2000. IMTC 2000. Proceedings of the 17th IEEE,volume 1, pages 24–29 vol.1, 2000. → pages 16, 29[25] something that i do not know. → pages 16, 17, 27[26] T. Hassan, F. Javed, and N. Arshad. An empirical investigation of v-i trajectorybased load signatures for non-intrusive load monitoring. Smart Grid, IEEETransactions on, 5(2):870–878, March 2014. → pages xiv, 16, 17, 18, 24, 25[27] A-Zoha. Non-intrusive load monitoring approaches for disaggregated energysensing: A survey, May 2010. → pages xiv, 16[28] M. Zeifman. Disaggregation of home energy display data using probabilisticapproach. Consumer Electronics, IEEE Transactions on, 58(1):23–31, February 2012. →pages xiv, 17, 23, 24[29] Ching-Lung Lin Hsueh-Hsien Chang, Hong-Tzer Yang. Load identification inneural networks for a non-intrusive monitoring of industrial electrical loads,December 2014. → pages 18[30] Yizheng Xu. Artificial intelligence based methodology for load disaggregation atbulk supply point, December 2014. → pages 18[31] Daniel Kroening. Computer aided verification, May 2014. → pages 18[32] D. Egarter, V.P. Bhuvana, and W. Elmenreich. Paldi: Online load disaggregation viaparticle filtering. Instrumentation and Measurement, IEEE Transactions on,64(2):467–477, Feb 2015. → pages 19, 20[33] M. Figueiredo, B. Ribeiro, and A. de Almeida. Electrical signal source separation vianonnegative tensor factorization using on site measurements in a smart home.Instrumentation and Measurement, IEEE Transactions on, 63(2):364–373, Feb 2014. →pages 19, 20149[34] Albeik. Unsupervised disaggregation of low frequency power measurements, May2010. → pages 19, 23[35] J.Zico Kolter and Matthew J. Johnson. Download page for GAISLER IP Cores, May2011. → pages xiv, 20, 21[36] M. Figueiredo, B. Ribeiro, and A. de Almeida. Electrical signal source separation vianonnegative tensor factorization using on site measurements in a smart home.Instrumentation and Measurement, IEEE Transactions on, 63(2):364–373, Feb 2014. →pages xiv, 23[37] Hsueh-Hsien Chang. Non-intrusive demand monitoring and load identification forenergy management systems based on transient feature analyses, May 2010. →pages xiv, xv, 26, 27, 28, 29[38] Christopher Laughman. Ieee power and energy magazine feature analyses, May2010. → pages xv, 27, 30[39] H.W. Dommel. Digital computer solution of electromagnetic transients insingle-and multiphase networks. Power Apparatus and Systems, IEEE Transactions on,PAS-88(4):388–399, April 1969. → pages 33, 34, 38[40] A comprehensive resource for emtdc transient analysis for pscad power systemsimulation. → pages 34[41] Charles A. Thompson. A study of numerical integration techniques for use in thecompanion circult method of transient circuit analysis. → pages xv, 36, 37, 38[42] Jorge Hollman. Step by step eigenvalue analysis with emtp discrete time solutions.→ pages xv, 41, 43[43] G. Groenewold. Optimal ladder filters. Circuits and Systems II: Express Briefs, IEEETransactions on, 56(2):147–151, Feb 2009. → pages 56[44] K.L. Su. Analog Filters. K.L. Su, K.L. Su, 2002. → pages 56[45] Henry Ruston. Analog Filters using MATLAB. ML. Wanhammar, L. Wanhammarl,2009. → pages 56[46] Active and passive filter synthesis using matlab. → pages 56150[47] K. Murthy and R.E. Bedford. Transformation between foster and cauer equivalentnetworks. Circuits and Systems, IEEE Transactions on, 25(4):238–239, Apr 1978. →pages 56[48] J.N. Davidson, D.A. Stone, and M.P. Foster. Required cauer network order formodelling of thermal transfer impedance. Electronics Letters, 50(4):260–262,February 2014. → pages 56[49] L. Codecasa. Canonical forms of one-port passive distributed thermal networks.Components and Packaging Technologies, IEEE Transactions on, 28(1):5–13, March 2005.→ pages 56[50] modern synthesis. → pages 56[51] The electrical engineering handbook. → pages 57, 58, 60[52] Synthesis of analogue circuits. → pages xv, 59[53] Henry Ruston. Electric Networks. McGraw-Hill, McGraw-Hill, 1966. → pages xv,xvi, 62, 63, 64, 65, 66[54] Peng Zhang. Shifted frequency analysis for emtp simulation of power systemdynamics, December 2009. → pages 146151

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0166755/manifest

Comment

Related Items