"Applied Science, Faculty of"@en .
"Electrical and Computer Engineering, Department of"@en .
"DSpace"@en .
"UBCV"@en .
"Chen, Bo"@en .
"2019-08-15T20:57:35Z"@en .
"2019"@en .
"Master of Applied Science - MASc"@en .
"University of British Columbia"@en .
"The widespread integration of distributed energy resources (DERs) in the future power system is likely to cause large and frequent excursions from its steady state. Such circumstances call for systematic methods to design DER controllers toward system-level goals and to estimate unexpected changes. We address these challenges via the development of a reduced-order power system dynamical model in this thesis.\r\n \r\nFirst, we develop a reduced second-order power-system dynamical model that accounts for locational effects of load disturbances on system frequency. The locational aspects are retained in the proposed model by incorporating linearized power-flow balance into differential equations that describe synchronous-generator dynamics. Individual synchronous-generator frequency dynamics are then combined into a single aggregate frequency state via weighting factors that can be tuned to maximize the accuracy of the reduced-order model. The proposed reduced-order model is general in the sense that its parameters are related to those of the original full-order model in analytical closed form, so that it can be constructed easily for different systems. Time-domain simulations demonstrate the accuracy of the reduced-order model with various choices of weighting factors and highlight the effect of load disturbance location on aggregate-frequency dynamics.\r\n \r\nThe reduced-order model offers the basic foundation on which DER controllers can be systematically designed to meet grid-level aims. Via a transfer-function representation of the reduced-order model, we relate, in analytical closed form, DER controller parameters to system steady-state and dynamic frequency response. Furthermore, time-domain simulations demonstrate that we can design DER controller parameter in order to meet system frequency-response performance requirements. \r\n \r\nBy leveraging the reduced-order power system dynamical model developed earlier, we propose an optimization-based method to detect the occurrence, estimate the magnitude, and identify the location of load changes. The proposed method relies on measurements of only frequency at the output of synchronous generators along with the reduced-order model that captures locational effects of load disturbances on generator frequency dynamics. The sparsity structure of load-change disturbances is leveraged so that fewer measurements are needed to estimate load changes. Furthermore, a convex relaxation of the problem ensures that it can be solved online in a computationally efficient manner."@en .
"https://circle.library.ubc.ca/rest/handle/2429/71307?expand=metadata"@en .
"A Reduced-order Power-system Dynamical Model andApplications in Design and Monitoring of Power SystemsbyBo ChenB. Eng., Sichuan University, 2015a 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)August 2019c\u00C2\u00A9 Bo Chen, 2019The following individuals certify that they have read, and recommend to the Faculty of Grad-uate and Postdoctoral Studies for acceptance, the thesis entitled:A Reduced-order Power-system Dynamical Model and Applications in De-sign and Monitoring of Power Systemssubmitted by Bo Chen in partial fulfillment of the requirements for the degree of MASTEROF APPLIED SCIENCE in Electrical and Computer Engineering.Examining Committee:Dr. Yu Christine Chen, Electrical and Computer EngineeringSupervisorDr. Juri Jatskevich, Electrical and Computer EngineeringSupervisory Committee MemberDr. Jose R. Marti, Electrical and Computer EngineeringSupervisory Committee MemberiiAbstractThe widespread integration of distributed energy resources (DERs) in the future power systemis likely to cause large and frequent excursions from its steady state. Such circumstances call forsystematic methods to design DER controllers toward system-level performance goals and toestimate unexpected changes in the system. We address these challenges via the developmentof a reduced-order power system dynamical model in this thesis.First, we develop a reduced second-order power-system dynamical model that accounts forlocational effects of load disturbances on system frequency. The locational aspects are retainedin the proposed model by incorporating linearized power-flow balance into differential equationsthat describe synchronous-generator dynamics. Individual synchronous-generator frequencydynamics are then combined into a single aggregate frequency state via weighting factors thatcan be tuned to maximize the accuracy of the reduced-order model. The proposed reduced-order model is general in the sense that its parameters are related to those of the originalfull-order model in analytical closed form, so that it can be constructed easily for differentsystems. Time-domain simulations demonstrate the accuracy of the reduced-order model withvarious choices of weighting factors and highlight the effect of load disturbance location onaggregate-frequency dynamics.The reduced-order model offers the basic foundation on which DER controllers can be sys-tematically designed to meet grid-level aims. Via a transfer-function representation of thereduced-order model, we relate, in analytical closed form, DER controller parameters to systemsteady-state and dynamic frequency response. Furthermore, time-domain simulations demon-strate that we can design DER controller parameter in order to meet system frequency-responseperformance requirements.iiiBy leveraging the reduced-order power system dynamical model developed earlier, we pro-pose an optimization-based method to detect the occurrence, estimate the magnitude, andidentify the location of load changes. The proposed method relies on measurements of onlyfrequency at the output of synchronous generators along with the reduced-order model thatcaptures locational effects of load disturbances on generator frequency dynamics. The sparsitystructure of load-change disturbances is leveraged so that fewer measurements are needed toestimate load changes. Furthermore, a convex relaxation of the problem ensures that it can besolved online in a computationally efficient manner.ivLay SummaryIn recent years, distributed energy resources (DERs) are assuming a more important role in themodern power system. These DERs pose significant challenges to operations and monitoringtasks that ensure overall system reliability. In this thesis, we develop an analytically tractablereduced-order power system dynamic model that retains attributes pertinent to the systemfrequency response. Using this model, we develop a method to design DER-controller param-eters in order to meet the system steady-state and dynamic-frequency response requirements.Furthermore, leveraging the reduced-order power-system model, we propose an optimizationmethod to detect the occurrence, estimate the magnitude, and identify the location of loadchanges with measurements of generator frequencies.vPrefaceThe majority of the research was conducted by the author of this thesis. My research advisor,Dr. Christine Chen, provided overall comments and editing during the process of conducting theresearch and writing this manuscript. Abdullah Al-Digs provided help in generating simulationresults in Chapters 2 and 4.A part of Chapter 2 has been published. B. Chen, A. Al-Digs, and Y. C. Chen, \u00E2\u0080\u009CA network-cognizant aggregate-frequency reduced-order power system dynamical model,\u00E2\u0080\u009D in Proceeding ofNorth American Power Symposium, Fargo, ND, September 2018.viTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiGlossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Background and Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2.1 Power-system Reduced-order Dynamical Models . . . . . . . . . . . . . . 31.2.2 Designing DER Controllers . . . . . . . . . . . . . . . . . . . . . . . . . . 31.2.3 Load-change Detection and Identification . . . . . . . . . . . . . . . . . . 41.3 Contributions of Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51.4 Thesis Organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 Reduced-order Power System Dynamical Model . . . . . . . . . . . . . . . . 7vii2.1 System Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.1.1 Power-flow Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.1.2 Synchronous-generator Dynamical Model . . . . . . . . . . . . . . . . . . 92.1.3 Power-system Dynamical Model . . . . . . . . . . . . . . . . . . . . . . . 102.2 Model-order Reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.2.1 Aggregate-frequency Model . . . . . . . . . . . . . . . . . . . . . . . . . . 142.2.2 Common-frequency Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.2.3 Aggregate-governor Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.3 Case Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.3.1 WECC 9-Bus Test System . . . . . . . . . . . . . . . . . . . . . . . . . . . 182.3.2 New England 39-bus Test System . . . . . . . . . . . . . . . . . . . . . . . 202.4 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233 DER Controller Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 243.1 Modified Reduced-order Model to Incorporate DERs . . . . . . . . . . . . . . . . 243.1.1 Frequency-responsive DER Model . . . . . . . . . . . . . . . . . . . . . . 253.1.2 Full State-space Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253.1.3 Model-order Reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 273.2 Transfer Function from Load Changes to Aggregate Frequency . . . . . . . . . . 293.3 Frequency Regulation Using DERs . . . . . . . . . . . . . . . . . . . . . . . . . . 293.3.1 Steady-state Response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313.3.2 Dynamic Response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313.3.3 Solution Strategy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 333.4 Case Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 343.4.1 Verifying Transfer Function . . . . . . . . . . . . . . . . . . . . . . . . . . 343.4.2 Designing Frequency Response . . . . . . . . . . . . . . . . . . . . . . . . 353.5 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 384 Load-change Detection Leveraging Frequency Measurements . . . . . . . . 394.1 Model Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 394.2 Load-change Disturbance Estimation and Identification . . . . . . . . . . . . . . 40viii4.2.1 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 404.2.2 Solution Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 424.3 Case Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 434.3.1 Simulation Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 444.3.2 Load-change Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 444.4 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 455 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 465.1 Contributions of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 465.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48A Model Details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54A.1 Detailing BGG , BGL, BLG , BLL Matrices . . . . . . . . . . . . . . . . . . . . . . . 54A.2 Showing H1G = 0G . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55B Parameters for Numerical Case Studies . . . . . . . . . . . . . . . . . . . . . 57ixList of TablesTable 2.1 Values of weighting factor c used in case studies involving the New England39-bus test system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21Table 3.1 Comparison of \u00E2\u0088\u0086\u00CF\u0089nadir [rad/s] and \u00E2\u0088\u0086\u00CF\u0089ss [rad/s]. . . . . . . . . . . . . . . . . . 37Table 3.2 Deff and Reff values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37Table 3.3 DD and MD values for designing \u00E2\u0088\u0086\u00CF\u0089ss and \u00E2\u0088\u0086\u00CF\u0089nadir. . . . . . . . . . . . . . . . 37xList of FiguresFigure 2.1 Network topology for two-bus system. . . . . . . . . . . . . . . . . . . . . . . 13Figure 2.2 Network topology for WECC 9-bus test system. . . . . . . . . . . . . . . . . 18Figure 2.3 WECC 9-bus test system: comparison of trajectories resulting from PSATsimulations and reduced-order models with various values of weighting factorsc and load disturbance at bus 5, \u00E2\u0088\u0086PL,5 = \u00E2\u0088\u00920.2 p.u. (a) Aggregate-frequencydeviations. (b) Aggregate-turbine-mechanical-power deviations. . . . . . . . . 19Figure 2.4 WECC 9-bus test system: comparison of trajectories resulting from PSATsimulations and reduced-order models with load disturbances at bus 5 or 8,which are distinguishable with c = c?, but not with c = MG . (a) Aggregate-frequency deviations. (b) Aggregate-turbine-mechanical-power deviations. . . 20Figure 2.5 Network topology for New Ebgland 39-bus test system. . . . . . . . . . . . . 21Figure 2.6 New England 39-bus test system: comparison of trajectories resulting fromPSAT simulations and reduced-order models with various values of weightingfactors c and load disturbance at bus 20, \u00E2\u0088\u0086PL,20 = \u00E2\u0088\u00921 p.u. (a) Aggregate-frequency deviations. (b) Aggregate-turbine-mechanical-power deviations. . . 22Figure 2.7 New England 39-bus test system: comparison of trajectories resulting fromPSAT simulations and reduced-order models with load disturbances at bus15 or 20, which are distinguishable with c = c?, but not with c = MG .(a) Aggregate-frequency deviations. (b) Aggregate-turbine-mechanical-powerdeviations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22Figure 3.1 System frequency response [1] . . . . . . . . . . . . . . . . . . . . . . . . . . . 30xiFigure 3.2 Modified WECC 9-bus test system with DER (marked as square). . . . . . . 34Figure 3.3 Comparison of trajectories resulting from transfer function and PSAT. . . . . 35Figure 3.4 Modified New England 39-bus test system with DERs (marked as squares). . 36Figure 3.5 \u00E2\u0088\u0086\u00CF\u0089nadir and \u00E2\u0088\u0086\u00CF\u0089ss design result. . . . . . . . . . . . . . . . . . . . . . . . . . . 36Figure 4.1 WECC test system: load change detection, estimation, and identificationalgorithm performance evaluation. (a)(b) Demonstrating that dynamic statetrajectories reconstructed from \u00E2\u0088\u0086\u00CF\u0089G [k] measurements (dashed traces) areaccurate as compared to actual trajectories (solid traces). . . . . . . . . . . . 44Figure 4.2 Estimated load changes due to a ramp-followed-by-step disturbance in loadat bus 6 using measurements \u00E2\u0088\u0086\u00CF\u0089G [k]. . . . . . . . . . . . . . . . . . . . . . . 45xiiGlossaryCOI Center-of-inertiaDAE Differential Algebraic EquationDER Distributed Energy ResourcePMU Phasor Measurement UnitVSG Virtual Synchronous GeneratorxiiiAcknowledgmentsI would like to express my deepest appreciation to my supervisor Professor Yu Christine Chen,for her guidance, support, and patience. She teaches me not only the knowledge but also theprofessional way to undertake research. It is my great honor to have such a knowledgeable andenthusiastic mentor.I also owe thanks to Abdullah Al-Digs, for his invaluable help in research. His attitudeand ability always inspire me when we worked together. Many thanks to Shuan Dong, for hisinsight and encouragement, and also for relaxing conversions to share daily life.Lastly, I owe special gratitude to my parents for their understanding, support, and love.xivChapter 1IntroductionThis chapter presents background, motivation, and related literature for developing reduced-order dynamical models of the power system and applications thereof. Also included in thischapter are the thesis contributions and structure.1.1 Background and MotivationPresently and further in the near future, distributed energy resources (DERs), such as rooftopsolar, small wind, and energy storage, are expected to gradually but widely displace conventionalcentralized fossil fuel-based generation [2]. While this transformation would help to alleviateenvironmental concerns associated with fossil fuels, it poses numerous challenges in the reliableand efficient operation of the integrated power grid [3]. For example, since renewable generationdepends on local weather conditions, such as cloud cover and wind speed, the power generatedis inherently intermittent, variable, and uncertain. These characteristics are at odds with theutility\u00E2\u0080\u0099s efforts to dispatch generation to exactly meet the customer electricity demand at alltimes. Moreover, intermittency and variability patently affect local voltage levels and maycause them to deviate from acceptable operating ranges [4]. On the other hand, furnished withappropriate strategies, power-electronic-based DER controllers offer more flexible capabilities tosupport and optimize the integrated power grid. To achieve this vision, the systematic design ofDER controllers is essential to promote static and dynamic power availability and quality [5].Moreover, the widespread integration of renewable generation will likely result in large andfrequent excursions away from the steady-state operating point. Particularly important under1these circumstances is the ability to estimate nodal-injection (or load) changes in a timelymanner. Indeed, the focus of this thesis is to develop a power-system model that enables oneto systematically design DER controllers and estimate load changes.The power system is a complex critical infrastructure that comprises dynamics arising fromnumerous device types (including synchronous generators, DERs, loads, etc.) coupled acrossa geographically expansive network. Thus, reduced-order models that capture only the mostrelevant attributes for a particular problem setting are prudent and often necessary due to lim-itations in analytical tool sets and computational resources. Model-order reduction has beenstudied extensively for a wide range of applications in the context of power system static anddynamic simulation, analysis, and control (see, e.g., [6] and references therein). For exam-ple, reduced-order models have been developed to study the power-flow problem [7], transientstability [8], and interarea oscillations [9], to name a few.In this thesis, we develop a reduced-order power system dynamical model that focuses ontime scales corresponding to system inertial and primary-frequency response. Such a modelis potentially useful to tackle a multitude of contemporary power-system problems, such asconducting dynamic contingency analysis and assessing frequency deviations due to renew-able generation variations. Particularly, we focus our attention to leverage the reduced-ordermodel in designing DER controllers to provide frequency regulation support and estimatingload changes with frequency measurements.1.2 Related WorkIn this thesis, we develop a reduced-order model that closely mimics the frequency dynamics ofthe actual system. The model helps us in designing DER controller parameters to meet systemdynamical and steady-state frequency requirements. We also propose a measurement-basedmethod to estimate and identify load changes by using reduced-order model. In this section, wereview previous work related to reduced-order power-system models and applications exploredin this thesis.21.2.1 Power-system Reduced-order Dynamical ModelsMost related prior work in power system dynamic model-order reduction can be broadly cat-egorized into algorithmic or numerical approaches. Algorithmic methods, such as repeatedtime-domain simulations [8], singular perturbation techniques [9], and weak-coupling identifi-cation [10], have been used to identify coherent groups of generators within a large-scale powersystem. These and application of numerical techniques, such as modal truncation [11], selectivemodal analysis [12], and Krylov-subspace methods [13] are generally deployed on a detailedhigher-order dynamical model. Typically, due to the algorithmic and numerical nature of thesemethods, they are unable to relate the parameters of the original full-order system to thoseof the reduced-order one once it is constructed. Recent work has addressed this issue by ana-lytically reducing the full-order model while retaining original-system parameters [14\u00E2\u0080\u009316]. Ourapproach begins with a similar set up as in [15] and [16], and the resulting reduced-order modelalso retains parameters of the original model in closed form. However, by embedding the net-work constraints within the reduced-order model, it can distinguish amongst dynamics arisingfrom load disturbances that occur at different locations in the network. Moreover, we introducea parameter that enables the flexibility of tuning the reduced-order model so that its dynamicsmore closely match those of the original full-order system.1.2.2 Designing DER ControllersVarious approaches have been proposed to address DER controller design problem to providesystem frequency support, which includes inertial- and primary-frequency response. Manypapers consider DER controller model according to the device type specifically. For exam-ple, [17, 18] propose methods to design inertial response for photovoltaic inverters. The effortsin [19\u00E2\u0080\u009321] involve inertial and primary-frequency design for double-fed induction generator inwind turbines. On the other hand, a more general controller structure is the so-called virtualsynchronous generator (VSG), which embeds the mathematical model of a synchronous gener-ator into the DER power-electronic converter [22\u00E2\u0080\u009324]. Particularly, the VSG is endowed withinertia and droop constants that are typically associated with synchronous generators. Partic-ularly, the DER synthetic-inertia constant is not constrained by the physical characteristics ofa rotating machine. Thus, the synthetic-inertia and damping constants can be designed so that3DERs help to achieve desired system-level operational requirements, such as steady-state anddynamic frequency excursions.Limited attention has been given to methods that rely on unveiling the analytical closed-form relationship that map DER inertia or damping constants to the steady-state and dynamicsystem frequency response. Constant-ramp-rate model for governors was applied in design-ing damping constants [25\u00E2\u0080\u009327]. The model can simplify the optimization problem, which thegenerator mechanical power changes linearly with frequency deviation. However, this modeldoes not accurately capture dynamic frequency response. Other efforts have been proposedto address the model accuracy problem in time-domain analysis [28]. Moreover, [29] proposean optimization program analytically derive the sensitivities of transient frequency overshootand damping ratio to inertia and damping. This method considers system eigenvalues andtheir sensitivities, which incurs large computation burden. A second-order aggregate-frequencymodel was leveraged to design inertia and damping constants in [1, 30]. This approach cap-tures the dynamics of synchronous generators, as well as frequency-responsive DERs endowedwith inertial and damping control. The method presented in this thesis follows a similar pathas in [1, 30], but by incorporating the network power-flow equations, we tailor the design ofDER-controller parameters for load disturbances at different locations in the system.1.2.3 Load-change Detection and IdentificationMost prior art has addressed the problem of detecting and identifying transmission-line out-ages [31\u00E2\u0080\u009334], power-quality disturbances [35\u00E2\u0080\u009337], and cascading events leading to blackouts [38\u00E2\u0080\u009341]. With respect to load estimation, nonintrusive load monitoring has been widely studied toanalyze load-state changes at the individual household appliance level [42\u00E2\u0080\u009344]. Unlike these,our work focuses on estimating load changes in the bulk transmission system by relying ononline measurements of synchronous-generator frequencies, which are commonly collected tocompute the system-wide frequency [45].Many countries are deploying phasor measurement units (PMUs) at the transmission levelto provide wide-area measurements for real-time monitoring [46]. The PMU measures voltageand current phasors at a very high frequency, and phasors measured at different locationsby different devices are time-synchronized [31]. A family of measurement-based methods is4based on Thevenin\u00E2\u0080\u0099s Theorem: local measurements at the monitored load bus or area areused to estimate a Thevenin equation approximating the rest of the system, i.e., a voltagesource connected through a Thevenin impedance; the power transfer to the monitored loadreaches its maximum when that Thevenin impedance has the same magnitude as the loadimpedance [47, 48].1.3 Contributions of ThesisSince we aim to capture inertial and primary-frequency response in a reduced-order power-system model, our approach begins by adopting a third-order model for each synchronousgenerator, capturing its rotor-angle position, electrical frequency, and governor dynamics. Wefurther formulate linearized power-flow equations and incorporate them into the generator dy-namical equations, thereby reducing the standard power-system differential-algebraic-equationmodel into one that contains only differential equations. In this way, the network topology isembedded within the original full-order system dynamical model, and the locational effects ofdifferent load disturbances on frequency response are captured. Next, by applying tuneableweighting factors to each synchronous-generator speed state, we obtain a single aggregate fre-quency that represents a weighted average of all generator speeds. Here, the weighting factorsmay be designed to, e.g., minimize the error between the reduced- and the full-order models.Next, considering DERs connected to the power system, we model the frequency-responsiveDER dynamics by a variant of the swing equation, which is similar to a simplified synchronousgenerator model. We further derive a reduced-order model, and generate the transfer functionthat reveals the system frequency response after different load disturbances. The transfer func-tion relates frequency dynamic overshoot and steady-state frequency offset with DER inertiaand damping constants. Then, leveraging the transfer function, we solve a system of nonlinearequations to compute DERs inertia and damping constants simultaneously. Previous workssolve the DER inertia and damping constants separately through trial and error process, whichincurs greater computational burden and a lot of time [49]. By tuning these DERs parame-ters, we can control the system steady-state and dynamic frequency response to meet a specificrequirement. The main advantage of this systematic solution method is that it results in accu-rate DER controller parameters quickly, so that DER parameters can be changed on the fly to5regulate system frequency in near-real time.Another application of the reduced-order dynamical model is to use it along with frequencymeasurements to detect and identify load changes. The proposed approach begins with a clas-sical perturbative power-system dynamical model around an initial operating point [50]. Usingdc power-flow assumptions, we embed the impact of load changes on synchronous-generatordynamics by algebraically manipulating and incorporating pertinent admittance-like matricesin the dynamical model [51]. Then, leveraging this model along with generator frequency mea-surements collected at a particular sampling instant, we solve a convex optimization problemto simultaneously estimate the magnitude and location of load changes in a computationallyefficient manner. The main advantage of the proposed method is that it relies on online mea-surements of synchronous-generator frequencies, which can be obtained from sensors, such asPMUs, equipped at buses connected to generators. Furthermore, unlike most previous workin power-system change detection and identification, we explicitly incorporate the underlyingsystem dynamics that give rise to measured quantities into the modelling framework, therebyaccomplishing the intended estimation task using measurements of dynamic states prior toreaching the post-disturbance steady state. We demonstrate the effectiveness of the proposedmethod using synthetic measurements periodically sampled from a time-domain simulation ofa detailed, lossy, and nonlinear power-system DAE model.1.4 Thesis OrganizationIn Chapter 2, we develop the reduced-order power system dynamic model which differentiatessystem frequency dynamics arising from load changes at different locations in the system.Model accuracy and performance are verified using the WECC 9-bus and New England 39-bus test systems. In Chapter 3, we incorporate DERs modelled as VSGs into the previousreduced-order model, which can be used to compute DER controller parameters and achievedesired steady-state and dynamic frequency response. Furthermore, using the reduced-orderdynamical model, we propose a method to estimate load changes with generator frequencymeasurements in Chapter 4. Finally, Chapter 5 offers concluding remarks by summarizing thetheoretical model development and application-driven results presented in this thesis.6Chapter 2Reduced-order Power SystemDynamical ModelIn this chapter, we propose a reduced second-order power system dynamical model that ac-counts for locational effects of load-change disturbances on system frequency dynamics overtime scales corresponding to inertial and primary-frequency response. The reduced-order modelretains attributes that are pertinent to system frequency response. As a direct consequence,the proposed model facilitates analytical development and reduces computational burden, es-pecially for large-scale power systems. A salient feature of the model is that it distinguishesfrequency dynamics arising from load-change disturbances at different locations in the system.2.1 System DescriptionConsider a power transmission network with N buses collected in the set N = {1, . . . , N} andtransmission lines in the set E \u00E2\u008A\u0082 N \u00C3\u0097N . Let Vk(t) = |Vk(t)|\u00E2\u0088\u00A0\u00CE\u00B8k(t) represent the voltage phasorat bus k and time t. Denote, by G = {1, . . . , G} \u00E2\u008A\u0082 N , the set of G buses that are connected toconventional turbine-based generators. Further denote, by L = N \G = {G+1, . . . , N}, the setof L load buses, which are connected to only non-frequency-sensitive loads, e.g., constant-powerloads. A transmission line with current flowing from bus k to ` is denoted by (k, `) \u00E2\u0088\u0088 E .72.1.1 Power-flow ModelTransmission line (k, `) is modelled using the lumped-element \u00CE\u00A0-model with series admittanceyk` = y`k = gk` + jbk` \u00E2\u0088\u0088 C \ {0} and shunt admittance yshk` = gshk` + jbshk` \u00E2\u0088\u0088 C \ {0} on both endsof the line. The active-power balance at bus k \u00E2\u0088\u0088 N is given by0 = Pk(t)\u00E2\u0088\u0092\u00E2\u0088\u0091`\u00E2\u0088\u0088NkPk`(t), (2.1)where Pk(t) is the net active-power injection at bus k, Nk is the set of buses electrically con-nected to bus k, and Pk`(t) represents the active-power flow in line (k, `).AC Power-flow FormulationIf bus k is a load bus, i.e., k \u00E2\u0088\u0088 L, thenPk(t) = PL,k(t), (2.2)where PL,k(t) is the non-frequency-sensitive active-power injection at bus k, which is a negativequantity if it corresponds to a constant-power load. On the other hand, if bus k is connectedto a synchronous generator, i.e., k \u00E2\u0088\u0088 G, thenPk(t) = PL,k(t) + Pek (t), (2.3)where, similar to (2.2), PL,k(t) is the non-frequency-sensitive active-power injection at bus k.As for P ek (t) in (2.3), using the classical synchronous-machine model for generator k \u00E2\u0088\u0088 G, whichcomprises a constant voltage Ek behind the transient reactance jX\u00E2\u0080\u00B2d,k, Pek (t) can be expressedas [50]P ek (t) =Ek|Vk(t)|X \u00E2\u0080\u00B2d,ksin(\u00CE\u00B4k(t)\u00E2\u0088\u0092 \u00CE\u00B8k(t)), (2.4)8where \u00CE\u00B4k(t) denotes the rotor electrical angular position of generator k. Finally, followingstandard power-flow computations, Pk`(t) in (2.1) is given byPk`(t) =|Vk(t)|2(gshk` + gk`)\u00E2\u0088\u0092 |Vk(t)||V`(t)|gk` cos(\u00CE\u00B8k(t)\u00E2\u0088\u0092 \u00CE\u00B8`(t))\u00E2\u0088\u0092 |Vk(t)||V`(t)|bk` sin(\u00CE\u00B8k(t)\u00E2\u0088\u0092 \u00CE\u00B8`(t)).(2.5)DC Power-flow ApproximationsIn transmission networks, the line resistance is much smaller than its reactance, i.e., for line(k, `), gk` << bk` and gshk` << bshk`, so we can approximate yk` = y`k \u00E2\u0089\u0088 jbk` and yshk` \u00E2\u0089\u0088 jbshk`. Also,under typical operating conditions, the voltage phase-angle differences between two buses aresmall, i.e., \u00CE\u00B8k` << 1, for all k, ` \u00E2\u0088\u0088 N , and \u00CE\u00B4g \u00E2\u0088\u0092 \u00CE\u00B8g << 1 for all g \u00E2\u0088\u0088 G. Under these assumptions,it follows that sin(\u00CE\u00B8k \u00E2\u0088\u0092 \u00CE\u00B8`) \u00E2\u0089\u0088 \u00CE\u00B8k \u00E2\u0088\u0092 \u00CE\u00B8` and sin(\u00CE\u00B4k \u00E2\u0088\u0092 \u00CE\u00B8k) \u00E2\u0089\u0088 \u00CE\u00B4k \u00E2\u0088\u0092 \u00CE\u00B8k. Finally, in the per-unit system,the numerical values of voltage magnitudes are typically near 1 p.u., i.e., |Vk| \u00E2\u0089\u0088 1, for all k \u00E2\u0088\u0088 N ,and Eg \u00E2\u0089\u0088 1, for all g \u00E2\u0088\u0088 G. Taken together, these so-called DC assumptions help to simplify (2.5)asPk`(t) = \u00E2\u0088\u0092bk`(\u00CE\u00B8k(t)\u00E2\u0088\u0092 \u00CE\u00B8`(t)). (2.6)Similarly, (2.4) simplifies asP ek (t) =1X \u00E2\u0080\u00B2d,k(\u00CE\u00B4k(t)\u00E2\u0088\u0092 \u00CE\u00B8k(t)). (2.7)2.1.2 Synchronous-generator Dynamical ModelDynamics of synchronous generator g \u00E2\u0088\u0088 G can be modelled as follows. For each generatorg \u00E2\u0088\u0088 G, let \u00CE\u00B4g(t) and \u00CF\u0089g(t) denote its rotor electrical angular position and frequency, respectively.Further let Pmg (t) denote its turbine mechanical power. Assume each generator initially operatesat the steady-state equilibrium point with \u00CF\u0089g(0) = \u00CF\u0089s = 2pi60 rad/s, the synchronous frequency.Then defining \u00E2\u0088\u0086\u00CF\u0089g := \u00CF\u0089g \u00E2\u0088\u0092 \u00CF\u0089s, dynamics of generator g \u00E2\u0088\u0088 G can be described by the following9third-order system:1\u00CE\u00B4\u00CB\u0099g(t) = \u00E2\u0088\u0086\u00CF\u0089g(t),Mg\u00E2\u0088\u0086\u00CF\u0089\u00CB\u0099g(t) = Pmg (t)\u00E2\u0088\u0092Dg\u00E2\u0088\u0086\u00CF\u0089g(t)\u00E2\u0088\u0092\u00E2\u0088\u0091`\u00E2\u0088\u0088NgPg`(t) + PL,g(t),\u00CF\u0084gP\u00CB\u0099mg (t) = Prg \u00E2\u0088\u0092 Pmg (t)\u00E2\u0088\u0092Rg\u00E2\u0088\u0086\u00CF\u0089g(t), (2.8)where Mg and Dg denote, respectively, the inertia and damping constants, and \u00CF\u0084g, Prg , andRg denote the governor time constant, reference power input, and inverse2 droop constant,respectively. The system in (2.8) consists of the classical synchronous-machine model combinedwith a speed-governor model [50].In (2.8), the generator dynamics are coupled to the network through the\u00E2\u0088\u0091`\u00E2\u0088\u0088Ng Pg`(t)term, thereby accounting for network effects on generator dynamics. Motivated by this, weaim to develop a corresponding reduced-order aggregate-frequency model for the entire system,which incorporates network effects. Before delving into these aspects below, we first combinesynchronous-generator dynamics in (2.8) with the approximate power-flow equations in (2.6)\u00E2\u0080\u0093(2.7).2.1.3 Power-system Dynamical ModelWith the individual-bus power balance and single-generator dynamic equations described thusfar, we next combine them into a full-order system dynamical model. We begin by defin-ing relevant vector variables as follows. For the set of synchronous generators G, collectrotor angular positions and variations in their speeds into vectors \u00CE\u00B4G = [\u00CE\u00B41, . . . , \u00CE\u00B4G]T and\u00E2\u0088\u0086\u00CF\u0089G = [\u00E2\u0088\u0086\u00CF\u00891, . . . ,\u00E2\u0088\u0086\u00CF\u0089G]T, respectively. Similarly, collect the voltage angles of generator busesand non-frequency-responsive injections at these buses into \u00CE\u00B8G = [\u00CE\u00B81, . . . , \u00CE\u00B8G]T and PG =[PL,1, . . . , PL,G]T, respectively. Analogously, collect these quantities for the set of load busesinto \u00CE\u00B8L = [\u00CE\u00B8G+1, . . . , \u00CE\u00B8N ]T and PL = [PL,G+1, . . . , PL,N ]T.1While more detailed models exist [50], we choose this simplified one because (i) it captures the dynamicbehaviour that we are interested in, i.e., inertial- and primary-frequency response, and (ii) it allows us to retainoriginal-system parameters in closed form in the model-order reduction process.2In most reference textbooks, Rg refers to the droop constant; in this paper, we deviate from the standard tocontain notational burden later.10Power-flow Algebraic EquationsUsing the notation established above, and by substituting (2.3) and (2.6) into (2.1), we canexpress the active-power balance at generator buses compactly in matrix form as30G = PG(t) + diag(1GX \u00E2\u0080\u00B2d)\u00CE\u00B4G(t)\u00E2\u0088\u0092BGG\u00CE\u00B8G(t)\u00E2\u0088\u0092BGL\u00CE\u00B8L(t), (2.9)where X \u00E2\u0080\u00B2d = [X\u00E2\u0080\u00B2d,1, . . . , X\u00E2\u0080\u00B2d,G]T, and matrices BGG and BGL are constructed appropriately us-ing (2.6) in conjunction with the network topology. Similarly, by substituting (2.2) and (2.6)into (2.1), the active-power balance at load buses can be compactly expressed as0L = PL(t)\u00E2\u0088\u0092BLG\u00CE\u00B8G(t)\u00E2\u0088\u0092BLL\u00CE\u00B8L(t), (2.10)where, again, entries of BLG and BLL are evaluated using (2.6) based on the network topology.To further reduce notational burden in (2.9) and (2.10), define \u00CE\u00B8 := [\u00CE\u00B8TG , \u00CE\u00B8TL ]T and P :=[PTG , PTL ]T. Then, we combine (2.9) and (2.10) and rearrange the resultant to get\u00CE\u00B8(t) = B\u00E2\u0088\u00921 (D\u00CE\u00B4G(t) + P (t)) , (2.11)whereB =\u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0BGG BGLBLG BLL\u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB , D =\u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0diag(1GX\u00E2\u0080\u00B2d)0L\u00C3\u0097G\u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB . (2.12)Entries in and structures of matrices BGG , BGL, BLG , and BLL are detailed in Appendix A.1.Note that by making use of partitioned matrix inversion, and assuming relevant sub-matricesare invertible, (2.11) can be simplified as [53]\u00CE\u00B8(t) = A\u00CE\u00B4G(t) +B\u00E2\u0088\u00921P (t), (2.13)3The M -dimensional vectors with all 0s and 1s are denoted by 0M and 1M , respectively. The diagonal matrixdiag(x) is formed with entries of the vector x stacked on the main diagonal; and diag(x/y) forms a diagonalmatrix with the ith diagonal entry given by xi/yi, where xi and yi are the ith entries of vectors x and y,respectively. The M by N matrix of all 0s is denoted by 0M\u00C3\u0097N .11whereA =\u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0 diag(1G)\u00E2\u0088\u0092B\u00E2\u0088\u00921LLBLG\u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB (BGG \u00E2\u0088\u0092BGLB\u00E2\u0088\u00921LLBLG)\u00E2\u0088\u00921diag(1GX \u00E2\u0080\u00B2d). (2.14)Synchronous-generator Dynamic EquationsFor the set of synchronous generators G, collect turbine mechanical power states and referencepower inputs into PmG = [Pm1 , . . . , PmG ]T and P rG = [Pr1 , . . . , PrG]T, respectively. Then, we canexpress generator dynamics in (2.8) compactly in matrix form as\u00CE\u00B4\u00CB\u0099G(t) = \u00E2\u0088\u0086\u00CF\u0089G(t), (2.15)diag(MG)\u00E2\u0088\u0086\u00CF\u0089\u00CB\u0099G(t) = PmG (t)\u00E2\u0088\u0092 diag(DG)\u00E2\u0088\u0086\u00CF\u0089G(t) +K\u00CE\u00B8(t) + PG(t), (2.16)diag(\u00CF\u0084G)P\u00CB\u0099mG (t) = PrG \u00E2\u0088\u0092 PmG (t)\u00E2\u0088\u0092 diag(RG)\u00E2\u0088\u0086\u00CF\u0089G(t), (2.17)whereMG = [M1, . . . ,MG]T, DG = [D1, . . . , DG]T, \u00CF\u0084G = [\u00CF\u00841, . . . , \u00CF\u0084G]T, andRG = [R1, . . . , RG]T.Also, in (2.16), the entries of K \u00E2\u0088\u0088 RG\u00C3\u0097L are formed by appropriately evaluating (2.6) based onthe network topology asK = \u00E2\u0088\u0092[diag(1GX\u00E2\u0080\u00B2d)+BGG BGL]. (2.18)Full-order Power-system ModelFinally, by substituting (2.11) into (2.16), we obtain the power-system dynamical model asfollows:\u00CE\u00B4\u00CB\u0099G(t) = \u00E2\u0088\u0086\u00CF\u0089G(t), (2.19)\u00E2\u0088\u0086\u00CF\u0089\u00CB\u0099G(t) = diag(MG)\u00E2\u0088\u00921(PmG (t)\u00E2\u0088\u0092 diag(DG)\u00E2\u0088\u0086\u00CF\u0089G(t) +H\u00CE\u00B4G(t) +WP (t)), (2.20)P\u00CB\u0099mG (t) = diag(\u00CF\u0084G)\u00E2\u0088\u00921 (P rG \u00E2\u0088\u0092 PmG (t)\u00E2\u0088\u0092 diag(RG)\u00E2\u0088\u0086\u00CF\u0089G(t)) , (2.21)12Figure 2.1: Network topology for two-bus system.where matrices H \u00E2\u0088\u0088 RG\u00C3\u0097G and W \u00E2\u0088\u0088 RG\u00C3\u0097N capture network effects and are expressed as,respectively,H = KB\u00E2\u0088\u00921D = KA, (2.22)W = KB\u00E2\u0088\u00921 +[diag(1G) 0G\u00C3\u0097L]. (2.23)Furthermore, with a view that entries of P (t) represent inputs to the power system, entries of Windicate the effects of load variations at different buses on the speed dynamics of each machine.Here, we should say explicitly that, even though there exist more detailed DAE models of thesystem, we will refer to the model in (2.19)\u00E2\u0080\u0093(2.21) as the full-order model4 in the rest of thethesis.Next, we illustrate the ideas presented above on forming the full-order power-system modelvia an example.Example 1 (Two-bus System) In this example, we consider the two-bus test system shownin Fig. 2.1, where each bus is connected to a synchronous generator as well as non-frequency-dependent injections. Using the notation established in (2.12), for this system, we haveB =\u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0 1X\u00E2\u0080\u00B2d,1 + 1Xl \u00E2\u0088\u0092 1Xl\u00E2\u0088\u0092 1Xl 1X\u00E2\u0080\u00B2d,1 +1Xl\u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB , D =\u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0 1X\u00E2\u0080\u00B2d,1 00 1X\u00E2\u0080\u00B2d,2\u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB .Furthermore, making use of (2.18), we get thatK =\u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0\u00E2\u0088\u0092 1Xl 1Xl1Xl\u00E2\u0088\u0092 1Xl\u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB .4We refer to the model in (2.19)\u00E2\u0080\u0093(2.21) as the \u00E2\u0080\u009Cfull-order\u00E2\u0080\u009D model as it will serve as the basis for subsequentmodel-order reduction. However, in the simulation case studies, more detailed machine models are utilized.13Then, using the expressions for B, D, and K above, we can evaluate (2.22) to obtainH =1X \u00E2\u0080\u00B2d,1 +Xl +X\u00E2\u0080\u00B2d,2\u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0\u00E2\u0088\u00921 11 \u00E2\u0088\u00921\u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB ,as expected via direct inspection of Fig. 2.1. Finally, using (2.23) for the two-bus system, weget thatW =1X \u00E2\u0080\u00B2d,1 +Xl +X\u00E2\u0080\u00B2d,2\u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0X \u00E2\u0080\u00B2d,2 +Xl X \u00E2\u0080\u00B2d,2X \u00E2\u0080\u00B2d,1 X\u00E2\u0080\u00B2d,1 +Xl\u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB .As shown above, under DC assumptions, the input matrix W is expressed as a function of onlythe network parameters. The entries in W highlight the effect of load variations at differentbuses on the speed dynamics of each synchronous machine. \u00042.2 Model-order ReductionIn this section, we sequentially present several reductions of the full-order system in (2.19)\u00E2\u0080\u0093(2.21). Our goal is to obtain an aggregate-frequency model that accounts for network effects.2.2.1 Aggregate-frequency ModelDefine system aggregate frequency \u00E2\u0088\u0086\u00CF\u0089\u00CB\u009C := cT\u00E2\u0088\u0086\u00CF\u0089G , where the gth entry in c \u00E2\u0088\u0088 RG representsthe contribution of generator g rotor angular speed to the aggregate frequency. From (2.20),we get that\u00E2\u0088\u0086 \u00CB\u0099\u00CB\u009C\u00CF\u0089(t) = cT\u00E2\u0088\u0086\u00CF\u0089\u00CB\u0099G(t) = cTdiag(MG)\u00E2\u0088\u00921(PmG (t)\u00E2\u0088\u0092 diag(DG)\u00E2\u0088\u0086\u00CF\u0089G(t) +H\u00CE\u00B4G(t) +WP (t)). (2.24)Also define scalar variable \u00CE\u00B4\u00CB\u009C := cTdiag(MG)\u00E2\u0088\u00921H\u00CE\u00B4G , and from (2.19), we have that\u00CB\u0099\u00CB\u009C\u00CE\u00B4(t) = cTdiag(MG)\u00E2\u0088\u00921H\u00E2\u0088\u0086\u00CF\u0089G(t). (2.25)The governor dynamics for all generators are retained as in (2.21). Later, in Section 2.2.3, wewill develop a reduced single-governor state for the entire system.Example 2 (Special Case\u00E2\u0080\u0093Canonical Two-area Reduced-order Model) In order to study14interarea oscillations in the canonical two-area power system, it is common practice to establisha two-machine interarea equivalent circuit [6]. This system\u00E2\u0080\u0099s structure is similar to the oneshown in Fig. 2.1, except that the loads PL,1 = PL,2 = 0. The swing equations for the two-machine equivalent can then be combined to form a reduced second-order system that is similarin form to the single-machine swing equation [6]. In this example, we show that, with appropri-ate choice of weighting factor c, we recover this reduced-order system, which is well establishedin the literature. To do so, choose c = [1,\u00E2\u0088\u00921]T, so that \u00E2\u0088\u0086\u00CF\u0089\u00CB\u009C = \u00E2\u0088\u0086\u00CF\u00891 \u00E2\u0088\u0092\u00E2\u0088\u0086\u00CF\u00892. For simplicity, inaddition to setting PL,1 = PL,2 = 0, also let damping constants D1 = D2 = 0. Then, (2.24) canbe expressed as\u00E2\u0088\u0086 \u00CB\u0099\u00CB\u009C\u00CF\u0089(t) = M\u00E2\u0088\u009211 Pm1 (t)\u00E2\u0088\u0092M\u00E2\u0088\u009212 Pm2 (t) + \u00CE\u00B4\u00CB\u009C(t), (2.26)and, based on (2.25), as well as H evaluated from Example 1,\u00CB\u0099\u00CB\u009C\u00CE\u00B4(t) = \u00E2\u0088\u0092 M\u00E2\u0088\u009211 +M\u00E2\u0088\u009212X \u00E2\u0080\u00B2d,1 +Xl +X\u00E2\u0080\u00B2d,2\u00E2\u0088\u0086\u00CF\u0089\u00CB\u009C(t). (2.27)Taken together, (2.27) and (2.26) form the proposed aggregate-frequency model with c = [1,\u00E2\u0088\u00921]T.Next, to show that this model is equivalent to the reduced second-order system model in [6],we define auxiliary angle variable\u00CE\u00B4 := \u00E2\u0088\u0092(M\u00E2\u0088\u009211 +M\u00E2\u0088\u009212X \u00E2\u0080\u00B2d,1 +Xl +X\u00E2\u0080\u00B2d,2)\u00E2\u0088\u00921\u00CE\u00B4\u00CB\u009C, (2.28)so that \u00CE\u00B4\u00CB\u0099 = \u00E2\u0088\u0086\u00CF\u0089\u00CB\u009C, (2.26) can then be expressed as\u00E2\u0088\u0086 \u00CB\u0099\u00CB\u009C\u00CF\u0089(t) = M\u00E2\u0088\u009211 Pm1 (t)\u00E2\u0088\u0092M\u00E2\u0088\u009212 Pm2 (t)\u00E2\u0088\u0092M\u00E2\u0088\u009211 +M\u00E2\u0088\u009212X \u00E2\u0080\u00B2d,1 +Xl +X\u00E2\u0080\u00B2d,2\u00CE\u00B4(t). (2.29)Furthermore, by defining aggregate quantitiesM :=M1M2M1 +M2, Pm:=M2Pm1 \u00E2\u0088\u0092M1Pm2M1 +M2, (2.30)we simplify (2.29) asM\u00E2\u0088\u0086 \u00CB\u0099\u00CB\u009C\u00CF\u0089(t) = Pm(t)\u00E2\u0088\u0092 1X \u00E2\u0080\u00B2d,1 +Xl +X\u00E2\u0080\u00B2d,2\u00CE\u00B4(t), (2.31)15the structure and form of which are reminiscent of the single-machine swing equation.Thus, indeed, the choice of c = [1,\u00E2\u0088\u00921]T in (2.24)\u00E2\u0080\u0093(2.25) results in the well-known reducedinter-area dynamical model in [6]. \u00042.2.2 Common-frequency ModelAssume that the electrical distances between geographically different parts of the network arenegligible, so that all generator frequencies follow the same transient behaviour [54]. As adirect consequence, in (2.8), \u00E2\u0088\u0086\u00CF\u0089g = \u00E2\u0088\u0086\u00CF\u0089, for all g \u00E2\u0088\u0088 G. Then, the aggregate-frequency dynamicsin (2.24) becomecT1G\u00E2\u0088\u0086\u00CF\u0089\u00CB\u0099(t) = cTdiag(MG)\u00E2\u0088\u00921(PmG (t)\u00E2\u0088\u0092DG\u00E2\u0088\u0086\u00CF\u0089(t) +H\u00CE\u00B4G(t) +WP (t)). (2.32)For simplicity, and without loss of generalization, we impose a normalization condition on theweighting factor c, so that cT1G = 1. Also define scalar variable \u00CE\u00B4 := cTdiag(MG)\u00E2\u0088\u00921H\u00CE\u00B4G .Then we obtain the following reduced (2 +G)th-order common-frequency model:\u00CE\u00B4\u00CB\u0099(t) = Heff\u00E2\u0088\u0086\u00CF\u0089(t), (2.33)\u00E2\u0088\u0086\u00CF\u0089\u00CB\u0099(t) = cTdiag(MG)\u00E2\u0088\u00921PmG (t)\u00E2\u0088\u0092Deff\u00E2\u0088\u0086\u00CF\u0089(t) + \u00CE\u00B4(t) +WeffP (t), (2.34)P\u00CB\u0099mG (t) = diag(\u00CF\u0084G)\u00E2\u0088\u00921 (P rG \u00E2\u0088\u0092 PmG (t)\u00E2\u0088\u0092RG\u00E2\u0088\u0086\u00CF\u0089(t)) , (2.35)where Heff \u00E2\u0088\u0088 R, Deff \u00E2\u0088\u0088 R, and Weff \u00E2\u0088\u0088 R1\u00C3\u0097N are expressed as, respectively,Heff := cTdiag(MG)\u00E2\u0088\u00921H1G, (2.36)Deff := cTdiag(MG)\u00E2\u0088\u00921DG , (2.37)Weff := cTdiag(MG)\u00E2\u0088\u00921W. (2.38)In this reduced-order model, Weff captures the effective locational effects of load changes acrossthe network on aggregate-frequency dynamics, and Deff represents the effective damping con-stant. Furthermore, and of particular note, Heff = 0 because H1G = 0G; we offer a proof of thisin Appendix A.2. Assuming that the network topology does not vary, \u00CE\u00B4\u00CB\u0099(t) \u00E2\u0089\u00A1 0, so \u00CE\u00B4(t) \u00E2\u0089\u00A1 \u00CE\u00B4(0),16for t > 0, and thus we can remove the dynamic state \u00CE\u00B4(t) from the system described by (2.33)\u00E2\u0080\u0093(2.35).2.2.3 Aggregate-governor ModelTo further reduce the model in (2.34)\u00E2\u0080\u0093(2.35) to a second-order system, we aggregate the governorstates by defining scalar variable Pm := cTdiag(MG)\u00E2\u0088\u00921PmG and, using (2.35), we getP\u00CB\u0099m(t) = cTdiag(MG)\u00E2\u0088\u00921P\u00CB\u0099mG (t)= cTdiag(MG)\u00E2\u0088\u00921diag(\u00CF\u0084G)\u00E2\u0088\u00921 \u00C2\u00B7(P rG \u00E2\u0088\u0092 PmG (t)\u00E2\u0088\u0092RG\u00E2\u0088\u0086\u00CF\u0089(t)). (2.39)In practice, the turbine-governor time constants are similar in value, i.e., \u00CF\u0084g \u00E2\u0089\u0088 \u00CF\u0084 , for all g \u00E2\u0088\u0088G [55]. Various options have been proposed for the aggregate parameter \u00CF\u0084 in the literature.The average of all entries in \u00CF\u0084G is utilized in [14, 56]. More recently, other choices have emergedbased on minimizing the error between the reduced- and higher-order models [15]. Assumingthat \u00CF\u0084g \u00E2\u0089\u0088 \u00CF\u0084 , for all g \u00E2\u0088\u0088 G, (2.39) can be simplified, and we arrive at the following reducedsecond-order model:\u00E2\u0088\u0086\u00CF\u0089\u00CB\u0099(t) = Pm(t)\u00E2\u0088\u0092Deff\u00E2\u0088\u0086\u00CF\u0089(t) + \u00CE\u00B4(0) +WeffP (t),P\u00CB\u0099m(t) = \u00CF\u0084\u00E2\u0088\u00921(P r \u00E2\u0088\u0092 Pm(t)\u00E2\u0088\u0092Reff\u00E2\u0088\u0086\u00CF\u0089(t)),(2.40)where the aggregate reference power and effective inverse droop constant are, respectively,P r := cTdiag(MG)\u00E2\u0088\u00921P rG , (2.41)Reff := cTdiag(MG)\u00E2\u0088\u00921RG . (2.42)The reduced-order model in (2.40) features the following key attributes: (i) it combines allsynchronous-generator frequencies into one aggregate frequency via the weighting factor c, and(ii) it is network cognizant as it distinguishes the effect of different load changes on systemaggregate-frequency dynamics via the vector Weff .172.3 Case StudiesIn this section, we apply the model-order reduction methodology described in this chapter totwo test systems: the Western Electric Coordinating Council (WECC) 3-machine 9-bus andNew England 39-bus test systems. The reduced second-order model is verified with time-domainsimulations of a nonlinear differential-algebraic model that includes the two-axis synchronous-generator, turbine-governor, and exciter models performed using PSAT [57]. Hence, in theremainder of this chapter, we refer to simulation results conducted using the detailed systemmodel as \u00E2\u0080\u009CPSAT\u00E2\u0080\u009D.Figure 2.2: Network topology for WECC 9-bus test system.2.3.1 WECC 9-Bus Test SystemVarying Weighting Factor cThe typology for WECC 9-bus test system is shown in Fig. 2.2. Consider the dynamic responsefollowing a load increase at bus 5, i.e., \u00E2\u0088\u0086PL,5(t) := PL,5(t)\u00E2\u0088\u0092 PL,5(0) = \u00E2\u0088\u00920.2 p.u., t > 0. Usingvarious values for the weighting factor c, we compare the mismatch in frequency and mechanical-power dynamics resulting from the full-order model versus the proposed reduced-order onein (2.40). As shown in Fig. 2.3 by the blue traces, the na\u00C2\u00A8\u00C4\u00B1ve choice of c = 13 [1, 1, 1]T (i.e., theaverage) results in large mismatch between the reduced- and full-order models. With c = MG ,on the other hand, the mismatch greatly reduces. Finally, as shown by green traces in Fig. 2.3,the choice c = c? = 125 [19, 4, 2]T (tuned via trial and error) attains nearly perfect match between18the dynamic response resulting from PSAT simulations and reduced-order models. Note that,in Fig. 2.3a, comparisons are made against the frequency deviations of generator 1, which isrepresentative of the other generators.0 1 2 3 4 5-0.7-0.6-0.5-0.4-0.3-0.2-0.10(a)0 1 2 3 4 5-101234 10-3(b)Figure 2.3: WECC 9-bus test system: comparison of trajectories resulting from PSATsimulations and reduced-order models with various values of weighting factors c andload disturbance at bus 5, \u00E2\u0088\u0086PL,5 = \u00E2\u0088\u00920.2 p.u. (a) Aggregate-frequency deviations.(b) Aggregate-turbine-mechanical-power deviations.Varying Location of Load DisturbanceSuppose generators in the WECC system are responding to a generation-load mismatch causedby a load increase located at bus 5 or 8, i.e., \u00E2\u0088\u0086PL,5(t) = \u00E2\u0088\u00920.2 p.u. or \u00E2\u0088\u0086PL,8(t) = \u00E2\u0088\u00920.2 p.u.,t > 0. As shown by the solid traces in Fig. 2.4, the two load disturbances lead to differentdynamics in generator frequencies and turbine mechanical powers when the detailed model isused to simulate the system. These locational effects are captured in the reduced-order modelwith c = c?, as demonstrated by dashed and dash-dot traces that nearly overlap with the solidones in Fig. 2.4. In contrast, the effect of disturbance location is not evident with c = Mg,because both load disturbances result in the same dynamic response, as shown by the reddotted traces in Fig. 2.4.190 1 2 3 4 5-0.7-0.6-0.5-0.4-0.3-0.2-0.10(a)0 1 2 3 4 500.20.40.60.811.21.41.61.8 10-3(b)Figure 2.4: WECC 9-bus test system: comparison of trajectories resulting from PSATsimulations and reduced-order models with load disturbances at bus 5 or 8, whichare distinguishable with c = c?, but not with c = MG . (a) Aggregate-frequencydeviations. (b) Aggregate-turbine-mechanical-power deviations.2.3.2 New England 39-bus Test SystemCenter-of-inertia QuantitiesIn Chapter 2.3.1, we compare aggregate-frequency dynamics with those of one generator, as itis representative of the other two generators in the system. However, this is not the case forthe New England 39-bus test system, which has 10 generators. Thus, we compare aggregate-frequency dynamics resulting from the proposed reduced-order model with the so-called center-of-inertia (COI) frequency, which is defined as\u00E2\u0088\u0086\u00CF\u0089COI :=\u00EF\u00A3\u00AB\u00EF\u00A3\u00AD\u00E2\u0088\u0091g\u00E2\u0088\u0088GMg\u00EF\u00A3\u00B6\u00EF\u00A3\u00B8\u00E2\u0088\u00921\u00E2\u0088\u0091g\u00E2\u0088\u0088GMg\u00E2\u0088\u0086\u00CF\u0089g. (2.43)We will compute \u00E2\u0088\u0086\u00CF\u0089COI by sampling electrical frequency deviations of all generators from thePSAT simulation, and this will serve as the benchmark against which we assess the applicabilityof the proposed reduced-order model.20Figure 2.5: Network topology for New Ebgland 39-bus test system.Varying Weighting Factor cThe typology for New England 39-bus test system is shown in Fig. 2.5. Similar to Chapter 2.3.1,consider the dynamic response following a load increase at bus 20, i.e., \u00E2\u0088\u0086PL,20(t) := PL,20(t)\u00E2\u0088\u0092PL,20(0) = \u00E2\u0088\u00921 p.u., t > 0. Using various values for the weighting factor c, we compare themismatch in frequency and mechanical-power dynamics resulting from the full-order modelversus the proposed reduced-order one in (2.40). As shown in Fig. 2.6 by the blue traces,choosing c = c\u00C2\u00AF (i.e., the average) results in some mismatch between the PSAT simulation andreduced-order models. With c = MG , on the other hand, the mismatch is reduced. Finally,as shown by green traces in Fig. 2.6, the choice c = c? (tuned via trial and error) attainsnearly perfect match between the dynamic response resulting from the PSAT simulation andreduced-order models. The value of c\u00C2\u00AF and c? are reported in Table. 2.1.Table 2.1: Values of weighting factor c used in case studies involving the New England39-bus test system.Generator g 1 2 3 4 5 6 7 8 9 10c\u00C2\u00AF 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100 0.100c? 0.096 0.109 0.109 0.109 0.109 0.109 0.109 0.109 0.071 0.071210 1 2 3 4 5 6-2.5-2-1.5-1-0.50(a)0 1 2 3 4 5 600.0020.0040.0060.0080.010.0120.0140.016(b)Figure 2.6: New England 39-bus test system: comparison of trajectories resulting fromPSAT simulations and reduced-order models with various values of weighting fac-tors c and load disturbance at bus 20, \u00E2\u0088\u0086PL,20 = \u00E2\u0088\u00921 p.u. (a) Aggregate-frequencydeviations. (b) Aggregate-turbine-mechanical-power deviations.0 1 2 3 4 5 6-2.5-2-1.5-1-0.50(a)0 1 2 3 4 5 600.0020.0040.0060.0080.010.0120.0140.016(b)Figure 2.7: New England 39-bus test system: comparison of trajectories resulting fromPSAT simulations and reduced-order models with load disturbances at bus 15 or20, which are distinguishable with c = c?, but not with c = MG . (a) Aggregate-frequency deviations. (b) Aggregate-turbine-mechanical-power deviations.Varying Location of Load DisturbanceSuppose generators in the New England 39-bus system are responding to a generation-loadmismatch caused by a load increase located at bus 15 or 20, i.e., \u00E2\u0088\u0086PL,15(t) = \u00E2\u0088\u00921 p.u. or22\u00E2\u0088\u0086PL,20(t) = \u00E2\u0088\u00921 p.u., t > 0. As shown by the solid traces in Fig. 2.7, the two load disturbanceslead to different dynamics in generator frequencies and turbine mechanical powers when thefull-order model is used to simulate the system. These locational effects are captured in thereduced-order model with c = c?, as demonstrated by dashed and dash-dot traces that nearlyoverlap with the solid ones in Fig. 2.7. In contrast, the effect of disturbance location is notevident with c = Mg, because both load disturbances result in the same dynamic response, asshown by the red dotted traces in Fig. 2.7.2.4 Concluding RemarksIn this chapter, we discussed a reduced-order model to estimate the system dynamic responsewith limited information. The introduction of weighting factor c enables the flexibility of tuningthe reduced-order model to match the dynamic behaviour resulting from PSAT simulations ofthe detailed model. Another crucial parameter Weff helps to differentiate dynamics arising fromload changes at different buses. Based on these two parameters, the proposed second-ordermodel is analytically tractable, and closed-form expressions reveal well-defined relationshipsbetween parameters of original and reduced systems. We apply the proposed method to theWECC 9-bus and New England 39-bus test systems, from which comparisons are made betweenresults obtained via the reduced-order and detailed models from PSAT. The results show thatthe proposed method can account for locational effects of load disturbances. In the next chapter,we will illustrate the application of this method for power system with DERs in order to designinertia and damping constants.23Chapter 3DER Controller DesignIn this chapter, we propose a method to design DER controllers that emulate synchronousgenerators. Specifically, the proposed method determines the synthetic-inertia and dampingconstants by leveraging a modified version of the reduced second-order power system dynamicalmodel from Chapter 2. The main advantage of the proposed method is that it involves a setof algebraic equations that can be solved simultaneously to obtain DER synthetic-inertia anddamping constants that achieve both desired steady-state and dynamic frequency response atthe same time.3.1 Modified Reduced-order Model to Incorporate DERsIn this section, we incorporate DERs into the power system described in Chapter 2.1. Then,using a similar derivation as in Chapter 2.2, we obtain a reduced second-order model for thepower system with DERs. To this end, consider the same system from Chapter 2 with Nbuses collected in the set N = {1, . . . , N}, and the set of G buses connected to conventionalturbine-based generators is denoted by G = {1, . . . , G}. Furthermore, denote the set of D busesconnected to frequency-responsive DERs (or their aggregates) by D = {G + 1, . . . , G + D}.Also denote the set of S generator and DER connected buses (collectively power sources) asS = D \u00E2\u0088\u00AA G, and to contain notational burden, we assume that D \u00E2\u0088\u00A9 G = \u00E2\u0088\u0085. Finally, the set ofL load buses are denoted by L = {G+D+ 1, . . . , N}. Other system descriptions and notationremain as in Chapter 2.243.1.1 Frequency-responsive DER ModelInspired by the model for synchronous generator in (2.8), assume the following model for eachDER d \u00E2\u0088\u0088 D:\u00CE\u00B4\u00CB\u0099d(t) = \u00E2\u0088\u0086\u00CF\u0089d(t),Md\u00E2\u0088\u0086\u00CF\u0089\u00CB\u0099d(t) = \u00E2\u0088\u0092Dd\u00E2\u0088\u0086\u00CF\u0089d(t)\u00E2\u0088\u0092\u00E2\u0088\u0091`\u00E2\u0088\u0088NdPd`(t) + PL,d(t), (3.1)where Md and Dd denote, respectively, the synthetic-inertia and damping constants for DER d,PL,d is the non-frequency-sensitive active power injected at bus d, and Pd` is the active powerflow from bus d to `. It is important to note that, unlike the inertia and damping constantsfor a synchronous generator, the Md and Dd are not constrained by the physics of a rotatingmass. Thus, we are free to assign values to tuneable parameters Md and Dd to achieve desiredpost-disturbance dynamic and steady-state behaviour.The model in (3.1) represents the dynamic behaviour of the DER (or an aggregation ofDERs, e.g., a large PV system [1]) connected to bus d. Taken together with the synchronous-generator and network models in Chapter 2.1, we develop the state-space model for a powersystem incorporating DERs next.3.1.2 Full State-space ModelWe combine conventional generators and frequency-responsive DERs into the same dynamicequations. To further reduce the notational burden, define \u00CE\u00B8 := [\u00CE\u00B8TG , \u00CE\u00B8TD, \u00CE\u00B8TL ]T. For the setof synchronous generators G and DERs D, collect their mechanical power states into PmS =[PmG , 0D]T = [Pm1 , . . . , PmG , 0D]T. Similarly, collect non-frequency-responsive injections at thesebuses into PS = [PL,G , PL,D]T = [PL,1, . . . , PL,G, PL,G+1, . . . , PL,G+D]T. Furthermore, collecttheir angular positions and frequency variations into vector \u00CE\u00B4S = [\u00CE\u00B4g, \u00CE\u00B4d]T = [\u00CE\u00B41, . . . , \u00CE\u00B4G, \u00CE\u00B4G+1,. . . , \u00CE\u00B4G+D]T and \u00E2\u0088\u0086\u00CF\u0089S = [\u00E2\u0088\u0086\u00CF\u0089G ,\u00E2\u0088\u0086\u00CF\u0089D]T = [\u00CF\u00891, . . . , \u00CF\u0089G, \u00CF\u0089G+1, . . . , \u00CF\u0089G+D]T, respectively. Then, wecan combine generator dynamics in (2.8) and DER dynamics in (3.1) to formulate the following25state-space model:\u00CE\u00B4\u00CB\u0099S(t) = \u00E2\u0088\u0086\u00CF\u0089S(t), (3.2)diag(MS)\u00E2\u0088\u0086\u00CF\u0089\u00CB\u0099S(t) = PmS (t)\u00E2\u0088\u0092 diag(DS)\u00E2\u0088\u0086\u00CF\u0089S(t) +K\u00CE\u00B8(t) + PS(t), (3.3)diag(\u00CF\u0084G)P\u00CB\u0099mG (t) = PrG \u00E2\u0088\u0092 PmG (t)\u00E2\u0088\u0092 diag(RG)\u00E2\u0088\u0086\u00CF\u0089G(t), (3.4)where MS = [MG ,MD]T = [M1, . . . ,MG,MG+1, . . . ,MG+D]T and DS = [DG , DD]T = [D1,. . . , DG, DG+1, . . . , DG+D]T. In (3.3), similar to (2.16), the entries of K \u00E2\u0088\u0088 RS\u00C3\u0097N are formed byappropriately evaluating (2.6) based on the network topology asK = \u00E2\u0088\u0092[diag(1SX\u00E2\u0080\u00B2d)+BSS BSL], (3.5)where the X \u00E2\u0080\u00B2d = [X\u00E2\u0080\u00B2d,1, . . . , X\u00E2\u0080\u00B2d,G, X\u00E2\u0080\u00B2d,G+1, . . . , X\u00E2\u0080\u00B2d,G+D]T. Here, we abuse notation slightly andkeep the same variable name X \u00E2\u0080\u00B2d for d-axis transient reactance. As DERs do not have tran-sient reactance, so we assume the values of X \u00E2\u0080\u00B2d for DERs are infinitely large (hence 1/X\u00E2\u0080\u00B2d isinfinitely small). Matrices BSS , BSL, BLS , and BLL are analogous to the matrices defined inAppendix A.1 and can be expressed asBSS = BSS \u00E2\u0088\u0092 diag([BSS BSL]1N +1SX \u00E2\u0080\u00B2d), (3.6)BSL = BSL, (3.7)BLS = BLS , (3.8)BLL = BLL \u00E2\u0088\u0092 diag([BLS BLL]1N), (3.9)where submatrices BSS , BSL, BLS , and BLL are extracted from the standard (and appropriatelyreordered) network admittance matrix Y for a lossless system, i.e.,Y = j\u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0BSS BSLBLS BLL\u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB . (3.10)26Define P := [PTS , PTL ]T. Then, following a similar derivation as in (2.9)\u00E2\u0080\u0093(2.12), we obtain thefollowing expression:\u00CE\u00B8(t) = B\u00E2\u0088\u00921 (D\u00CE\u00B4S(t) + P (t)) , (3.11)where1B =\u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0BSS BSLBLS BLL\u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB , D =\u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0diag(1SX\u00E2\u0080\u00B2d)0L\u00C3\u0097S\u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB . (3.12)Finally, by substituting (3.11) into (3.3), we obtain the full-order2 power-system dynamicalmodel as follows:\u00CE\u00B4\u00CB\u0099S(t) = \u00E2\u0088\u0086\u00CF\u0089S(t), (3.13)diag(MS)\u00E2\u0088\u0086\u00CF\u0089\u00CB\u0099S(t) = PmS (t)\u00E2\u0088\u0092 diag(DS)\u00E2\u0088\u0086\u00CF\u0089S(t) +H\u00CE\u00B4S(t) +WP (t), (3.14)diag(\u00CF\u0084G)P\u00CB\u0099mG (t) = PrG \u00E2\u0088\u0092 PmG (t)\u00E2\u0088\u0092 diag(RG)\u00E2\u0088\u0086\u00CF\u0089G(t), (3.15)where matrices H \u00E2\u0088\u0088 RS\u00C3\u0097S and W \u00E2\u0088\u0088 RS\u00C3\u0097N capture network effects and are expressed as,respectively,H = KB\u00E2\u0088\u00921D, (3.16)W = KB\u00E2\u0088\u00921 +[diag(1S) 0S\u00C3\u0097L]. (3.17)The model in (3.13)\u00E2\u0080\u0093(3.15) captures frequency dynamics for a power system that containssynchronous generators and DERs with controllers that mimic synchronous generators.3.1.3 Model-order ReductionAssume that all generators and DERs follow the same frequency dynamics, so that \u00E2\u0088\u0086\u00CF\u0089s = \u00E2\u0088\u0086\u00CF\u0089,for all s \u00E2\u0088\u0088 S. Furthermore, define scalar variable Pm := cTdiag(MS)\u00E2\u0088\u00921PmS , where the c isthe weighting factor introduced in Chapter 2. Assuming that \u00CF\u0084g \u00E2\u0089\u0088 \u00CF\u0084 , for all g \u00E2\u0088\u0088 G, and thenfollowing the same model-order reduction method as described in Chapter 2.2, we arrive at the1We use the same notation for matrices B and D as Chapter 2 to retain the equations in similar form asbefore. However, the entries in these two matrices are modified to incorporate DER buses.2We refer to the model in (3.13)\u00E2\u0080\u0093(3.15) as the \u00E2\u0080\u009Cfull-order\u00E2\u0080\u009D model as it will serve as the basis for subsequentmodel-order reduction. However, in the simulation case studies, more detailed machine models are utilized.27following second-order model:\u00E2\u0088\u0086\u00CF\u0089\u00CB\u0099(t) = Pm(t)\u00E2\u0088\u0092Deff\u00E2\u0088\u0086\u00CF\u0089(t) + \u00CE\u00B4(0) +WeffP (t),P\u00CB\u0099m(t) = \u00CF\u0084\u00E2\u0088\u00921(P r \u00E2\u0088\u0092 Pm(t)\u00E2\u0088\u0092Reff\u00E2\u0088\u0086\u00CF\u0089(t)),(3.18)where the Deff \u00E2\u0088\u0088 R, Weff \u00E2\u0088\u0088 R1\u00C3\u0097N , aggregate reference power and effective inverse dampingconstant are expressed as, respectively,Deff := cTdiag(MS)\u00E2\u0088\u00921DS , (3.19)Weff := cTdiag(MS)\u00E2\u0088\u00921W, (3.20)P r := cTdiag(MS)\u00E2\u0088\u00921[P rG , 0D]T, (3.21)Reff := cTdiag(MS)\u00E2\u0088\u00921[RG ,0D]T. (3.22)Furthermore, we rewrite (3.18) in state space, asx\u00CB\u0099red = Aredxred +Bredured. (3.23)The state vector and input, xred \u00E2\u0088\u0088 R2, ured \u00E2\u0088\u0088 RN+2, system matrices, Ared \u00E2\u0088\u0088 R2\u00C3\u00972, Bred \u00E2\u0088\u0088R2\u00C3\u0097(N+2) are given byxred = [\u00E2\u0088\u0086\u00CF\u0089(t), Pm(t)]T , ured = [P (t), Pr, \u00CE\u00B4(0)]T ,Ared =\u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0 \u00E2\u0088\u0092Deff 1\u00E2\u0088\u0092Reff\u00CF\u0084\u00E2\u0088\u00921 \u00E2\u0088\u0092\u00CF\u0084\u00E2\u0088\u00921\u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB , Bred =\u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0 Weff 0 10TN \u00CF\u0084\u00E2\u0088\u00921 0\u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB . (3.24)This model contains the aggregated frequency response for both synchronous generators andDERs, furthermore it also retains aggregated governor dynamics for generators. The term \u00CE\u00B4(0)is constant, and it does not affect the relationship between \u00E2\u0088\u0086\u00CF\u0089(t) and P (t), as we will shownext.283.2 Transfer Function from Load Changes to AggregateFrequencyWe focus on the input-output relationship from the load disturbance to the aggregate frequencyin order to design DER synthetic-inertia and damping constants to meet requirements on fre-quency response. As the state-space model in (3.23) is a linear system, we can rewrite it asa transfer function. Particularly, s-domain transfer function from the load P (t) to aggregate-frequency deviation \u00E2\u0088\u0086\u00CF\u0089(t) can be obtained from (3.23) as\u00E2\u0088\u0086\u00CF\u0089(s) =(s+ \u00CF\u0084\u00E2\u0088\u00921)Weffs2 + (\u00CF\u0084\u00E2\u0088\u00921 +Deff)s+ (Reff +Deff)\u00CF\u0084\u00E2\u0088\u00921P (s)=(s+ \u00CE\u00B6)ks2 + 2\u00CE\u00BE\u00CF\u0089ns+ \u00CF\u00892nP (s), (3.25)where the parameters k, \u00CE\u00B6, \u00CE\u00BE, and \u00CF\u0089n are given byk := Weff , \u00CE\u00B6 := \u00CF\u0084\u00E2\u0088\u00921,\u00CF\u0089n :=\u00E2\u0088\u009A\u00CF\u0084\u00E2\u0088\u00921(Deff +Reff), \u00CE\u00BE :=12Deff + \u00CF\u0084\u00E2\u0088\u00921\u00E2\u0088\u009A\u00CF\u0084\u00E2\u0088\u00921(Reff +Deff). (3.26)The transfer function (3.25) represents the net aggregate-frequency deviations due to simulta-neous load changes at multiple buses. To see this, we emphasize that P (s) \u00E2\u0088\u0088 RN is a vectorcontains all buses instead of a scalar value. Moreover, the locational effects of any individualload change at bus i \u00E2\u0088\u0088 N is captured by the vector quantity k \u00E2\u0088\u0088 R1\u00C3\u0097N .3.3 Frequency Regulation Using DERsIn order to motivate ideas presented in this section, we refer to Fig. 3.1 in the following dis-cussion. This figure depicts the frequency response for a power system that only includessynchronous generators. The original system synchronous frequency is \u00CF\u0089s, then the introduc-tion of a load disturbance causes the system frequency to oscillate. At time tnadir, the systemfrequency reaches a minimum value, so-called frequency nadir \u00E2\u0088\u0086\u00CF\u0089nadir. Then generators invokeprimary frequency response, which leads the system to a new stable frequency after some time.Compared to the original system frequency, the steady-state frequency deviation is \u00E2\u0088\u0086\u00CF\u0089ss. Ourgoal is to achieve the desired dynamic and steady-state system frequency response using DERs.29Figure 3.1: System frequency response [1]Specifically, by leveraging the transfer-function model in (3.25), we solve for DER synthetic-inertia and damping constants that result in predefined requirements on \u00E2\u0088\u0086\u00CF\u0089ss and \u00E2\u0088\u0086\u00CF\u0089nadir.Suppose we introduce load step change \u00E2\u0088\u0086P at time t = 0, which includes an individual loadchange at bus i \u00E2\u0088\u0088 N . By substituting P (s) = 1s\u00E2\u0088\u0086P into (3.25), we get\u00E2\u0088\u0086\u00CF\u0089(s) =(s+ \u00CE\u00B6)k\u00E2\u0088\u0086Ps (s2 + 2\u00CE\u00BE\u00CF\u0089ns+ \u00CF\u00892n). (3.27)Assume that the system is underdamped, inverse Laplace transform of (3.27) yields\u00E2\u0088\u0086\u00CF\u0089(t) = \u00E2\u0088\u0086\u00CF\u0089ss(1\u00E2\u0088\u0092 e\u00E2\u0088\u0092\u00CE\u00BE\u00CF\u0089nt\u00E2\u0088\u009A1\u00E2\u0088\u0092 \u00CE\u00BE2(sin (\u00CF\u0089dt+ \u00CF\u0095)\u00E2\u0088\u0092 \u00CF\u0089n\u00CE\u00B6sin (\u00CF\u0089dt))), (3.28)where parameters \u00CF\u0089d and \u00CF\u0095, as well as steady-state frequency offset \u00E2\u0088\u0086\u00CF\u0089ss are expressed as,respectively,\u00CF\u0089d = \u00CF\u0089n\u00E2\u0088\u009A1\u00E2\u0088\u0092 \u00CE\u00B62, (3.29)\u00CF\u0095 = tan\u00E2\u0088\u00921(\u00CE\u00B6\u00E2\u0088\u00921\u00E2\u0088\u009A1\u00E2\u0088\u0092 \u00CE\u00B62), (3.30)\u00E2\u0088\u0086\u00CF\u0089ss =\u00CE\u00B6k\u00E2\u0088\u0086P\u00CF\u00892n. (3.31)303.3.1 Steady-state ResponseWe can further simplify (3.31) by substituting (3.26) to get\u00E2\u0088\u0086\u00CF\u0089ss =Weff\u00E2\u0088\u0086PReff +Deff. (3.32)Notice from (3.32) that the steady-state frequency offset \u00E2\u0088\u0086\u00CF\u0089ss is related to Deff , Weff , and Reff .Of the three parameters, Weff as expressed in (3.20) is a fixed quantity based on the networktopology. On the other hand, according to (3.19) and (3.22), Deff and Reff are functions of DSand MS , which include the inertia and damping constants for both generators and DERs, i.e.,MS = [MG ,MD]T and DS = [DG , DD]T. Thus, we can suitably assign values to Md and Dd,d \u00E2\u0088\u0088 D, which result in values of Meff and Deff that achieve desired steady-state frequency offset\u00E2\u0088\u0086\u00CF\u0089ss. Since there are two tuneable parameters, there are conceivably infinitely many ways toobtain the desired steady-state frequency offset. Below, we constrain the feasible solution spaceby invoking an additional requirement on the frequency nadir.3.3.2 Dynamic ResponseAn important attribute of the dynamic frequency response is the frequency nadir. As shown inFig. 3.1, the system frequency reaches \u00E2\u0088\u0086\u00CF\u0089nadir whend\u00E2\u0088\u0086\u00CF\u0089dt\u00E2\u0088\u00A3\u00E2\u0088\u00A3t=tnadir= 0 for the first time followingthe load disturbance. Indeed, we can solve for tnadir by taking the derivative of (3.28) settingthe resultant to zero, to obtain the following:0 = e\u00E2\u0088\u0092\u00CE\u00BE\u00CF\u0089ntnadir (\u00CF\u0089n\u00CE\u00BE sin (\u00CF\u0089dtnadir + \u00CF\u0095)\u00E2\u0088\u0092 \u00CF\u0089d cos (\u00CF\u0089dtnadir + \u00CF\u0095)+ \u00CE\u00B6\u00E2\u0088\u00921\u00CF\u00892n\u00CE\u00BE sin (\u00CF\u0089dtnadir)\u00E2\u0088\u0092 \u00CE\u00B6\u00E2\u0088\u00921\u00CF\u0089n\u00CF\u0089d cos (\u00CF\u0089dtnadir)).(3.33)Elementary trigonometric manipulations help to simplify the above to0 = \u00CF\u0089n\u00CE\u00BE (sin(\u00CF\u0089dtnadir) cos\u00CF\u0095+ cos(\u00CF\u0089dtnadir) sin\u00CF\u0095)\u00E2\u0088\u0092 \u00CF\u0089d (cos(\u00CF\u0089dtnadir) cos\u00CF\u0095\u00E2\u0088\u0092 sin (\u00CF\u0089dtnadir) sin\u00CF\u0095)+ \u00CE\u00B6\u00E2\u0088\u00921\u00CF\u00892n\u00CE\u00BE sin (\u00CF\u0089dtnadir)\u00E2\u0088\u0092 \u00CE\u00B6\u00E2\u0088\u00921\u00CF\u0089n\u00CF\u0089d cos (\u00CF\u0089dtnadir) .(3.34)31Then, dividing both sides by cos(\u00CF\u0089dtnadir) yields0 = \u00CF\u0089n\u00CE\u00BE (tan(\u00CF\u0089dtnadir) cos\u00CF\u0095+ sin\u00CF\u0095)\u00E2\u0088\u0092 \u00CF\u0089d (cos\u00CF\u0095\u00E2\u0088\u0092 tan(\u00CF\u0089dtnadir) sin\u00CF\u0095)+ \u00CE\u00B6\u00E2\u0088\u00921\u00CF\u00892n\u00CE\u00BE tan (\u00CF\u0089dtnadir)\u00E2\u0088\u0092 \u00CE\u00B6\u00E2\u0088\u00921\u00CF\u0089n\u00CF\u0089d.(3.35)Furthermore, according to (3.30), we have thatsin\u00CF\u0095 =\u00E2\u0088\u009A1\u00E2\u0088\u0092 \u00CE\u00B62, (3.36)cos\u00CF\u0095 = \u00CE\u00B6, (3.37)which can be substituted into (3.35) to yield0 = \u00CF\u0089n\u00CE\u00BE(tan(\u00CF\u0089dtnadir)\u00CE\u00B6 +\u00E2\u0088\u009A1\u00E2\u0088\u0092 \u00CE\u00B62)\u00E2\u0088\u0092 \u00CF\u0089d(\u00CE\u00B6 \u00E2\u0088\u0092 tan(\u00CF\u0089dtnadir)\u00E2\u0088\u009A1\u00E2\u0088\u0092 \u00CE\u00B62)+ \u00CE\u00B6\u00E2\u0088\u00921\u00CF\u00892n\u00CE\u00BE tan (\u00CF\u0089dtnadir)\u00E2\u0088\u0092 \u00CE\u00B6\u00E2\u0088\u00921\u00CF\u0089n\u00CF\u0089d.(3.38)Finally, we can simplify the above astan (\u00CF\u0089dtnadir) =\u00CF\u0089d\u00CE\u00BE\u00CF\u0089n \u00E2\u0088\u0092 \u00CE\u00B6 , (3.39)from which we can solve for tnadir astnadir = \u00CF\u0089\u00E2\u0088\u00921d tan\u00E2\u0088\u00921(\u00CF\u0089d\u00CE\u00BE\u00CF\u0089n \u00E2\u0088\u0092 \u00CE\u00B6). (3.40)In order to emphasize the relationship between tnadir and tuneable parameters Deff and Reff ,we substitute (3.26) and (3.29) into (3.40) to gettnadir = \u00CF\u0081\u00E2\u0088\u00921 tan\u00E2\u0088\u00921(2\u00CF\u0081Deff \u00E2\u0088\u0092 \u00CF\u0084\u00E2\u0088\u00921), (3.41)where\u00CF\u0081 =\u00E2\u0088\u009A(\u00CF\u0084\u00E2\u0088\u00921 \u00E2\u0088\u0092 \u00CF\u0084\u00E2\u0088\u00923)(Reff +Deff). (3.42)32Substituting \u00E2\u0088\u0086\u00CF\u0089ss and tnadir from (3.32) and (3.41), respectively, into (3.28), we obtain\u00E2\u0088\u0086\u00CF\u0089nadir = \u00E2\u0088\u0086\u00CF\u0089 (tnadir )=Weff\u00E2\u0088\u0086PReff +Deff(1\u00E2\u0088\u0092 e\u00E2\u0088\u0092\u00CE\u00BE\u00CF\u0089ntnadir\u00E2\u0088\u009A1\u00E2\u0088\u0092 \u00CE\u00BE2(sin (\u00CF\u0089dtnadir + \u00CF\u0095)\u00E2\u0088\u0092 \u00CF\u0089n\u00CE\u00B6sin (\u00CF\u0089dtnadir))). (3.43)When we design system dynamical response by (3.41) and (3.43), it involves k, \u00CE\u00BE, \u00CE\u00B6, and \u00CF\u0089n.From (3.26), these parameters are relate to Deff , Weff and Reff , which are then related to thetuneable DER synthetic-inertia and damping constants.3.3.3 Solution StrategyBy giving specific requirements on both \u00E2\u0088\u0086\u00CF\u0089ss and \u00E2\u0088\u0086\u00CF\u0089nadir, we utilize the steady-state re-sponse (3.32) and dynamic response (3.43) together to design the tuneable parameters Md andDd, d \u00E2\u0088\u0088 D. Considering the power system operational constraints for steady-state frequencyand frequency nadir, we obtain the following system of equations:\u00EF\u00A3\u00B1\u00EF\u00A3\u00B4\u00EF\u00A3\u00B2\u00EF\u00A3\u00B4\u00EF\u00A3\u00B3\u00E2\u0088\u0086\u00CF\u0089ss =Weff\u00E2\u0088\u0086PReff+Deff,\u00E2\u0088\u0086\u00CF\u0089nadir =Weff\u00E2\u0088\u0086PReff+Deff(1\u00E2\u0088\u0092 e\u00E2\u0088\u0092\u00CE\u00BE\u00CF\u0089ntnadir\u00E2\u0088\u009A1\u00E2\u0088\u0092\u00CE\u00BE2(sin (\u00CF\u0089dtnadir + \u00CF\u0095)\u00E2\u0088\u0092 \u00CF\u0089n\u00CE\u00B6 sin (\u00CF\u0089dtnadir))),(3.44)s.t. 0 \u00E2\u0089\u00A4 \u00E2\u0088\u0086\u00CF\u0089ss \u00E2\u0089\u00A4 \u00E2\u0088\u0086\u00CF\u0089max, (3.45)0 \u00E2\u0089\u00A4 \u00E2\u0088\u0086\u00CF\u0089nadir \u00E2\u0089\u00A4 \u00E2\u0088\u0086\u00CF\u0089max, (3.46)\u00E2\u0088\u0086\u00CF\u0089ss \u00E2\u0089\u00A4 \u00E2\u0088\u0086\u00CF\u0089nadir, (3.47)where the \u00E2\u0088\u0086\u00CF\u0089max is the maximum safety frequency drop for the power system. Since (3.44)represents a nonlinear set of equations, it may have multiple solutions, so the constraints (3.45)\u00E2\u0080\u0093(3.47) are needed to establish feasible operating points. The frequency deviation range is usuallyheld within \u00C2\u00B11%, so feasible frequency values range from 59.4 Hz to 60.6 Hz for the 60 Hz-gridin north America [58, 59].The method proposed above computes the DER synthetic-inertia and damping constants atthe same time to achieve the desired system frequency performance. The equations in (3.44)\u00E2\u0080\u0093(3.47) are nonlinear, so we can not express Md and Dd in closed-form expression. We can obtainthe numerical solution by using nonlinear equation solver in Matlab.333.4 Case StudiesIn this section, we verify the transfer function relationship between load disturbance and ag-gregate frequency by a modified WECC 9-bus test system. Furthermore, we demonstrate theproposed DER-controller design method via New England 39-bus test system.3.4.1 Verifying Transfer FunctionIn this section, we verify that the transfer function obtained in (3.25) matches the PSAT time-domain simulations conducted with detailed machine models. We modify the WECC 9-bus testsystem and replace the generator at bus 3 with a DER, as shown in Fig. 3.2. We embed themodel in (3.1) into the DER. Furthermore, we set the weighting factor c = c? = 125 [19, 4, 2]T asdetailed in Chapter. 2.3.1. Pertinent model parameters and ratings are listed in Appendix B.For convenience, power and impedance values are in per unit [pu] with base 100 [MVA], unlessotherwise specified.The time-domain simulation are performed using PSAT [57] for the modified WECC 9-bustest system shown in Fig. 3.2. At time t = 0, the load at bus 5 undergoes a step increase of\u00E2\u0088\u0086P = 0.2 p.u. We evaluate the transfer function and pertinent parameters from (3.25)\u00E2\u0080\u0093(3.32)using two values of weighting factor c: (i) c = c? as obtained in Chapter 2.3.1, and (ii) the naivechoice of setting c = MS . The transfer function result with c are obtained from (3.25)\u00E2\u0080\u0093(3.32).Consider the trajectories in Fig. 3.3, the blue solid trace corresponds to simulations ofthe transfer function c = c?, while the red dashed trace corresponds to those carried out withFigure 3.2: Modified WECC 9-bus test system with DER (marked as square).340 1 2 3 4 5-0.8-0.6-0.4-0.20Figure 3.3: Comparison of trajectories resulting from transfer function and PSAT.transfer function c = MS . The green dotted line is the dynamic response from PSAT simulationand serves as the benchmark. We observe that trajectories from the transfer function with c = c?are very similar to the benchmark result. In steady-state, the trajectory from transfer functionwith c = MS does not match the PSAT simulation. As the weighting factor c give the modeltuning flexibility and Weff have location information, the transfer function with c = c? matchesthe benchmark result very well.3.4.2 Designing Frequency ResponseIn this section, we design tuneable DER synthetic-inertia and damping constants to achievedesired frequency response using the method outlined in Chapter. 3.3. To this end, we modifythe New England 39-bus test system and replace generators 8\u00E2\u0080\u009310 with DERs as shown inFig. 3.4. The inertia and damping constants of DERs are unknown variables in this setting, andwe determine suitable values for them in order to meet the desired system frequency response.Pertinent model parameters and ratings are listed in Appendix B. For convenience, power andimpedance values are in per unit [pu] with base 100 [MVA], unless otherwise specified.We design system steady-state response \u00E2\u0088\u0086\u00CF\u0089ss and dynamic frequency response \u00E2\u0088\u0086\u00CF\u0089nadir atthe same time based on the modified New England 39-bus test system. At time t = 0, theload at bus 20 undergoes a step increase of \u00E2\u0088\u0086P = 1 p.u. We solve the system of equationsin (3.44)\u00E2\u0080\u0093(3.47) to compute the effective inertia and damping constants, which then help toinform the DER synthetic-inertia and damping constants. We use the optimal weighting factor35Figure 3.4: Modified New England 39-bus test system with DERs (marked as squares).0 1 2 3 4 5 6-2.5-2-1.5-1-0.50Figure 3.5: \u00E2\u0088\u0086\u00CF\u0089nadir and \u00E2\u0088\u0086\u00CF\u0089ss design result.c = c? obtained in Chapter 2.3.2.We investigate three settings with different frequency control aims: (i) \u00E2\u0088\u0086\u00CF\u0089nadir = \u00E2\u0088\u00921.80 rad/sand \u00E2\u0088\u0086\u00CF\u0089ss = \u00E2\u0088\u00920.80 rad/s, (ii) \u00E2\u0088\u0086\u00CF\u0089nadir = \u00E2\u0088\u00922.00 rad/s and \u00E2\u0088\u0086\u00CF\u0089ss = \u00E2\u0088\u00921.00 rad/s, (iii) \u00E2\u0088\u0086\u00CF\u0089nadir =\u00E2\u0088\u00922.20 rad/s and \u00E2\u0088\u0086\u00CF\u0089ss = \u00E2\u0088\u00921.20 rad/s. We plot resulting frequency trajectories obtained in eachof the three cases in Fig. 3.5. We observe that, indeed, the DER synthetic-inertia and damp-ing constants obtained via the proposed method are able to achieve prescribed dynamic and36steady-state frequency objectives. The time at which frequency nadir occurs, i.e., tnadir, varieswith different \u00E2\u0088\u0086\u00CF\u0089nadir values, because tnadir and \u00E2\u0088\u0086\u00CF\u0089nadir are related by (3.43). Furthermore, wereport the accuracy of the proposed method in Table 3.1. The comparison shows that both thesteady-state and dynamic frequency response meet the design requirements with good accu-racy. The error between desired and actual values can be accounted for by the approximationstaken to obtain the reduced-order system model, which is used to determine DER-controllerparameters. For example, as the weighting factor c obtained by trial and error, it introduceserror into the resulting model. Effective parameters Deff and Reff values are shown in Table 3.2.Tuneable parameters DD and MD are shown in Table 3.3.Table 3.1: Comparison of \u00E2\u0088\u0086\u00CF\u0089nadir [rad/s] and \u00E2\u0088\u0086\u00CF\u0089ss [rad/s].Case 1 2 3Desired \u00E2\u0088\u0086\u00CF\u0089nadir -1.80 -2.00 -2.20Actual \u00E2\u0088\u0086\u00CF\u0089nadir -1.8290 -1.997 -2.261Desired \u00E2\u0088\u0086\u00CF\u0089ss -0.80 -1.00 -1.20Actual \u00E2\u0088\u0086\u00CF\u0089ss -0.789 -1.006 -1.202Table 3.2: Deff and Reff values.Case 1 2 3Deff 3.4496 2.4043 1.8050Reff 8.6377 5.0304 3.8431Table 3.3: DD and MD values for designing \u00E2\u0088\u0086\u00CF\u0089ss and \u00E2\u0088\u0086\u00CF\u0089nadir.D8 D9 D10 M8 M9 M10\u00E2\u0088\u0086\u00CF\u0089nadir = \u00E2\u0088\u00921.80, \u00E2\u0088\u0086\u00CF\u0089ss = \u00E2\u0088\u00920.80 9.84 10.80 10.78 2.63 5.30 5.19\u00E2\u0088\u0086\u00CF\u0089nadir = \u00E2\u0088\u00922.00, \u00E2\u0088\u0086\u00CF\u0089ss = \u00E2\u0088\u00921.00 18.60 12.42 11.96 2.62 7.76 11.67\u00E2\u0088\u0086\u00CF\u0089nadir = \u00E2\u0088\u00922.20, \u00E2\u0088\u0086\u00CF\u0089ss = \u00E2\u0088\u00921.20 23.33 24.59 28.88 23.61 26.22 25.37373.5 Concluding RemarksIn this chapter, we proposed a method to regulate system frequency using DERs along with amodified version of reduced second-order power system dynamical model. The solution strat-egy proposed computes the synthetic-inertia and damping constants for DERs according toprescribed steady-state and dynamic frequency response. As the transfer function is the foun-dation of this method, we verify it via numerical case studies involving the WECC 9-bus testsystem. We apply the proposed design method to the New England 39-bus test system. Resultsshow that we indeed obtain synthetic-inertia and damping constants that lead to the desiredfrequency response.38Chapter 4Load-change Detection LeveragingFrequency MeasurementsIn this chapter, we develop a method to detect and identify load changes in the power systemusing only synchronous generator frequency measurements in near-real time. Central to theproposed detection method is the use of power-system model developed in Chapter. 2.4.1 Model FormulationConsider the same power system full-order model from Chapter 2 with nodes collected in theset N = G \u00E2\u0088\u00AA L, where G and L denotes the sets of synchronous generator and load buses,respectively. Transmission lines are collected in the set E \u00E2\u008A\u0082 N \u00C3\u0097N .In Chapter 2.1.3, the power-system dynamical model was introduced in (2.19)\u00E2\u0080\u0093(2.21). Thedetection method focus on the load-change compared to the previous time point, so we aim todevelop a perturbative model around an initial operating point, we will find the definition ofthe following small-signal variables for generator g useful: variations in rotor angular position\u00E2\u0088\u0086\u00CE\u00B4g(t) := \u00CE\u00B4g(t) \u00E2\u0088\u0092 \u00CE\u00B4g(0), variations in turbine mechanical power \u00E2\u0088\u0086Pmg (t) := Pmg (t) \u00E2\u0088\u0092 Pmg (0),and variations in loads \u00E2\u0088\u0086P (t) = P (t)\u00E2\u0088\u0092 P (0). For the set of synchronous generators G, collectvariations in rotor angular positions into vectors \u00E2\u0088\u0086\u00CE\u00B4G = [\u00E2\u0088\u0086\u00CE\u00B41, . . . ,\u00E2\u0088\u0086\u00CE\u00B4G]T and variations inturbine mechanical power states into \u00E2\u0088\u0086PmG = [\u00E2\u0088\u0086Pm1 , . . . ,\u00E2\u0088\u0086PmG ]T. With these definitions inplace, the large-signal model in (2.19)\u00E2\u0080\u0093(2.21) can be cast as a perturbative model that helps to39identify nonzero entries in \u00E2\u0088\u0086P (t), as follows:\u00E2\u0088\u0086\u00CE\u00B4\u00CB\u0099G(t) = \u00E2\u0088\u0086\u00CF\u0089G(t), (4.1)diag(MG)\u00E2\u0088\u0086\u00CF\u0089\u00CB\u0099G(t) = \u00E2\u0088\u0086PmG (t)\u00E2\u0088\u0092 diag(DG)\u00E2\u0088\u0086\u00CF\u0089G(t) +H\u00E2\u0088\u0086\u00CE\u00B4G(t) +W\u00E2\u0088\u0086P (t), (4.2)diag(\u00CF\u0084G)\u00E2\u0088\u0086P\u00CB\u0099mG (t) = \u00E2\u0088\u0086PrG \u00E2\u0088\u0092\u00E2\u0088\u0086PmG (t)\u00E2\u0088\u0092 diag(RG)\u00E2\u0088\u0086\u00CF\u0089G(t), (4.3)where \u00E2\u0088\u0086P rG represents variations in reference power inputs, and matrices H \u00E2\u0088\u0088 RG\u00C3\u0097G and W \u00E2\u0088\u0088RG\u00C3\u0097N are expressed as, respectively,H = KB\u00E2\u0088\u00921D = KA, (4.4)W = KB\u00E2\u0088\u00921 +[diag(1G) 0G\u00C3\u0097L]. (4.5)With a view that \u00E2\u0088\u0086P (t) in (4.2) represents a disturbance input to the power system, entries ofW indicate the effects of load-variation disturbances at different buses on generator frequencydynamics. Next, we outline a method to estimate load changes in the network using generatorfrequency measurements in conjunction with (4.2).4.2 Load-change Disturbance Estimation and IdentificationIn this section, we outline a method to estimate load changes and identify their locations usingmeasurements of generator frequencies. We further leverage the sparsity structure of the load-variation vector to estimate load disturbances using fewer measurements.4.2.1 Problem FormulationSuppose we have measurements of variations in synchronous-generator frequencies sampled att = k\u00E2\u0088\u0086t, k \u00E2\u0088\u0088 Z, where \u00E2\u0088\u0086t > 0 is the time interval between consecutive samples. Denotethe measurement at time t = k\u00E2\u0088\u0086t (or time step k) by \u00E2\u0088\u0086\u00CF\u0089G [k]. Then, at time step k, we canrearrange terms in (4.2) to get\u00E2\u0088\u0086y[k] = W\u00E2\u0088\u0086P [k], (4.6)40where\u00E2\u0088\u0086y[k] = diag(MG)\u00E2\u0088\u0086\u00CF\u0089\u00CB\u0099G [k]\u00E2\u0088\u0092\u00E2\u0088\u0086PmG [k] + diag(DG)\u00E2\u0088\u0086\u00CF\u0089G [k]\u00E2\u0088\u0092H\u00E2\u0088\u0086\u00CE\u00B4G [k]. (4.7)Our goal is to estimate the value of \u00E2\u0088\u0086P [k] that satisfies (4.6) by evaluating \u00E2\u0088\u0086y[k] using onlinemeasurements of \u00E2\u0088\u0086\u00CF\u0089G [k]. In addition to measurements of \u00E2\u0088\u0086\u00CF\u0089G [k], (4.7) requires samples of\u00E2\u0088\u0086\u00CF\u0089\u00CB\u0099G [k], \u00E2\u0088\u0086\u00CE\u00B4G [k], and \u00E2\u0088\u0086PmG [k], which are easily computed using measurements of \u00E2\u0088\u0086\u00CF\u0089G [k], asshown below.Computing \u00E2\u0088\u0086\u00CF\u0089\u00CB\u0099G [k]We approximate \u00E2\u0088\u0086\u00CF\u0089\u00CB\u0099G [k] as the slope between two consecutive measurements: \u00E2\u0088\u0086\u00CF\u0089G [k \u00E2\u0088\u0092 1] and\u00E2\u0088\u0086\u00CF\u0089G [k]. Specifically,\u00E2\u0088\u0086\u00CF\u0089\u00CB\u0099G [k] \u00E2\u0089\u0088 \u00E2\u0088\u0086t\u00E2\u0088\u00921 (\u00E2\u0088\u0086\u00CF\u0089G [k]\u00E2\u0088\u0092\u00E2\u0088\u0086\u00CF\u0089G [k \u00E2\u0088\u0092 1]) . (4.8)Computing \u00E2\u0088\u0086\u00CE\u00B4G [k]Discretization of the differential equation in (4.1) results in the following expression:\u00E2\u0088\u0086\u00CE\u00B4G [k + 1] = \u00E2\u0088\u0086\u00CE\u00B4G [k] + \u00E2\u0088\u0086t\u00E2\u0088\u0086\u00CF\u0089G [k]. (4.9)Given the initial value \u00E2\u0088\u0086\u00CE\u00B4G [0] = 0G, we can compute \u00E2\u0088\u0086\u00CE\u00B4G [k] at all time steps k > 0.Computing \u00E2\u0088\u0086PmG [k]We first note that variations in reference power inputs, denoted by \u00E2\u0088\u0086P rG(t), in (4.3) are theoutput of automatic generation control (AGC), which is used to restore system frequency backto synchronous value. In the timescales that we consider, we assume that the reference powerinputs do not change significantly, i.e., \u00E2\u0088\u0086P rG(t) \u00E2\u0089\u0088 0G, t \u00E2\u0089\u00A5 0, in (4.3). Then, discretizationof (4.3) yields\u00E2\u0088\u0086PmG [k + 1] = A\u00E2\u0088\u0086PmG [k] +B\u00E2\u0088\u0086\u00CF\u0089G [k], (4.10)41where matrices A,B \u00E2\u0088\u0088 RG\u00C3\u0097G are given byA = exp(\u00E2\u0088\u0092diag(\u00CF\u0084G)\u00E2\u0088\u00921\u00E2\u0088\u0086t),B = diag(RG)A\u00E2\u0088\u0092 diag(RG).Then, given the initial value \u00E2\u0088\u0086PmG [0] = 0G, we can compute \u00E2\u0088\u0086PmG [k], k > 0, using the recurrencerelation in (4.10).We close this discussion with a few remarks on the estimation of \u00E2\u0088\u0086P [k] from (4.6). Bysuitably shifting time indices in (4.8), (4.9), and (4.10), and substituting them into (4.7), weget\u00E2\u0088\u0086y[k] =(\u00E2\u0088\u0086t\u00E2\u0088\u00921diag(MG) + diag(DG))\u00E2\u0088\u0086\u00CF\u0089G [k]\u00E2\u0088\u0092 (\u00E2\u0088\u0086t\u00E2\u0088\u00921diag(MG) +Bd + \u00E2\u0088\u0086tH)\u00E2\u0088\u0086\u00CF\u0089G [k \u00E2\u0088\u0092 1]\u00E2\u0088\u0092Ad\u00E2\u0088\u0086PmG [k \u00E2\u0088\u0092 1]\u00E2\u0088\u0092H\u00E2\u0088\u0086\u00CE\u00B4G [k \u00E2\u0088\u0092 1]. (4.11)The observation in (4.11) indeed requires only measurements of synchronous-generator frequen-cies, and information at time step k \u00E2\u0088\u0092 1 is known from previous frequency samples. However,the estimation problem in (4.6) cannot be solved by na\u00C2\u00A8\u00C4\u00B1vely inverting W . This is because,in general, not all buses are connected to a generator, i.e., G < N , hence W is not a squarematrix. In fact, the estimation problem in (4.6) is under-determined, and it cannot be solvedvia conventional least-squares estimation. One way to bypass this issue is to assume that theload variations remain constant and stack up observations over multiple time steps and formu-late an over-determined estimation problem. However, in an effort to minimize the number ofmeasurements needed, we next outline a method to solve the estimation problem by leveragingthe expectedly sparse nature of load changes in the system.4.2.2 Solution AlgorithmAt any given point in time, we expect the vast majority of loads in the system to remain closeto their steady-state operating point, i.e., for a given time step k \u00E2\u0088\u0088 Z, \u00E2\u0088\u0086P [k] is a sparse vector.Thus, we cast the load-change detection and identification problem as an optimization program42that minimizes the number of nonzero elements in \u00E2\u0088\u0086P [k]:min\u00E2\u0088\u0086P [k]\u00E2\u0080\u0096\u00E2\u0088\u0086P [k]\u00E2\u0080\u00960s.t. \u00E2\u0088\u0086y[k] = W\u00E2\u0088\u0086P [k]. (4.12)The problem in (4.12) is NP-hard due to the unavoidable combinatorial search [60]. A commonsolution approach is to relax the nonconvex and discontinuous `0-norm cost function with theconvex and continuous `1-norm, yielding the following convex optimization problem:min\u00E2\u0088\u0086P [k]\u00E2\u0080\u0096\u00E2\u0088\u0086P [k]\u00E2\u0080\u00961s.t. \u00E2\u0080\u0096\u00E2\u0088\u0086y[k]\u00E2\u0088\u0092W\u00E2\u0088\u0086P [k]\u00E2\u0080\u009622 \u00E2\u0089\u00A4 \u000F, (4.13)where \u000F \u00E2\u0089\u00A5 0 is a small tolerance introduced to facilitate computations. The optimization prob-lem in (4.13) is known as the basis pursuit denoising (BPDN) problem and is widely usedin compressed sensing for signal reconstruction, where noise is inherent in signal measure-ments [61]. There are several algorithms that can be utilized to solve the BPDN problemin (4.13) such as the alternating direction method of multipliers (ADMM) [62] and the spectralprojected gradient for `1-norm minimization (SPGL1) algorithm [63, 64]. For \u00E2\u0088\u0086y[k] = W\u00E2\u0088\u0086P [k]to have at least one solution, the matrix W must be full rank. We utilize the SPGL1 algorithmin [63] to solve the BPDN problem in (4.13), which relies only on matrix-vector operations andis computationally tractable for large-scale estimation problems [64].4.3 Case StudiesWe demonstrate the proposed method via numerical case studies involving the Western Elec-tricity Coordinating Council (WECC) 3-machine 9-bus test system; the one-line diagram forthis system can be found in Fig. 2.2. Synchronous generators are connected to buses collectedin G = {1, 2, 3}, loads are connected to buses in L = {5, 6, 8}.434.3.1 Simulation SetupAlthough our analytical development is grounded in a simplified synchronous-generator modeland leverages several approximations (e.g., lossless, small voltage-angle differences, etc.), we ver-ify the proposed method with time-domain simulations for a detailed, lossy, and nonlinear DAEmodel of the power network that includes generator dynamics from a two-axis synchronous-machine model, a turbine-governor model, and an exciter model. The simulations are per-formed using PSAT [57]. Synthetic measurements are sampled from the PSAT simulation atdiscrete intervals of \u00E2\u0088\u0086t = 0.0333 s, which is within the capability of current measurement tech-nology [65]. In our setup, at time t = 0, the load at bus 6 increases by \u00E2\u0088\u0086P6 = 0.07 p.u. as astep change, which is followed by a ramp increase at t = 2 s to reach a net steady-state changeof \u00E2\u0088\u0086P6 = 0.16 p.u. at t = 3 s. Other loads remain constant.4.3.2 Load-change EstimationWe restrict measurements sampled from the time-domain simulation to only \u00E2\u0088\u0086\u00CF\u0089G [k] and com-pute \u00E2\u0088\u0086\u00CF\u0089\u00CB\u0099G [k], \u00E2\u0088\u0086\u00CE\u00B4G [k], and \u00E2\u0088\u0086PmG [k] using (4.8), (4.9), and (4.10), respectively. As shown inFig. 4.1a and Fig. 4.1b, the computed dynamic state trajectories match the actual ones rea-0 1 2 3 400.010.020.030.040.050.060.07(a)0 1 2 3 4-1.4-1.2-1-0.8-0.6-0.4-0.20(b)Figure 4.1: WECC test system: load change detection, estimation, and identificationalgorithm performance evaluation. (a)(b) Demonstrating that dynamic state tra-jectories reconstructed from \u00E2\u0088\u0086\u00CF\u0089G [k] measurements (dashed traces) are accurate ascompared to actual trajectories (solid traces).440 1 2 3 400.050.10.150.2 Bus 5Bus 6Bus 8Figure 4.2: Estimated load changes due to a ramp-followed-by-step disturbance in loadat bus 6 using measurements \u00E2\u0088\u0086\u00CF\u0089G [k].sonably well. With the measured \u00E2\u0088\u0086\u00CF\u0089G [k], \u00E2\u0088\u0086y[k] is evaluated using (4.11). Then, at each timestep, we solve the optimization problem in (4.13) using the SPGL1 algorithm [63, 64]. Theestimated load changes are plotted in Fig. 4.2. The load change is detected and its locationis identified at k = 1. The load changes estimated from measurements obtained during thetransient period match the actual values very closely. In steady state, load change at bus 6 isestimated to be \u00E2\u0088\u0086P\u00CC\u00826 = 0.1608 p.u., while others remain at zero, as desired. Indeed, the proposedload-change detection and identification method using only generator frequency measurementsis very effective.4.4 Concluding RemarksIn this chapter, we proposed a method to estimate load changes using generator frequencymeasurements along with a reduced-order power system dynamical model. The model accountsfor locational effects of load disturbances on generator frequency dynamics. The utility of theproposed method was demonstrated via numerical case studies involving the WECC test system.Results show that we can accurately detect and identify load disturbances in the network.45Chapter 5Conclusions5.1 Contributions of the ThesisIn this thesis, we propose two applications based on reduced-order power system dynamicmodel. The first part is that using DERs to design system after load change performance. Thesecond part is that load change detection and identification.In Chapter 2, we present a reduced second-order power system dynamical model, which ac-counts for locational effects of load disturbances on aggregate-frequency dynamics correspondingto the inertial- and primary-frequency response time scales. The utility of the proposed modelwas demonstrated via numerical case studies involving the WECC 9-bus and New England39-bus test system. There are two aspects of improvement compared to the most relevantmodel developed in prior work. First, by choosing the contribution of each generator speedto the aggregate-frequency response, we have the flexibility of tuning the reduced-order modelso that its dynamical response more closely matches that of the original full-order system.Furthermore, by embedding the power-flow constraints within the reduced-order model, it candistinguish between dynamics arising from load disturbances that occur at different locationsin the network.Chapter 3 outlines an approach to design DER synthetic-inertial and damping constantsbased on the previous second-order model. We can meet a specific requirement for systemsteady-state and dynamic response after load disturbance by adjusting the DER synthetic-inertia and damping constants. In order to solve for the DER constants numerically, we obtain46the power system transfer function described the system dynamics after load change. Theweighting factor provides us a way to disaggregate the synthetic value into individual DERsparameters. The effectiveness of this method was verified via numerical case studies.In Chapter 4, we proposed a method to estimate load changes using generator frequencymeasurements along with a reduced-order power system dynamical model. The model accountsfor locational effects of load disturbances on generator frequency dynamics. We obtain thesolution using a convex relaxation of the sparse-signal recovery problem. The utility of theproposed method was demonstrated via numerical case studies involving the WECC 9-bus testsystem. Results show that we can accurately detect and identify load disturbances in thenetwork.5.2 Future WorkThere are several potential areas for extension of this thesis work:\u00E2\u0080\u00A2 The weighting factor introduced in Chapter 2 is shown to be effective for second-orderpower system dynamical model. However, the value is obtained by trial and error. Thetuning process would be more efficient with an automated program, perhaps by castingthe problem as an optimization method. The general idea is to solve for the optimalweighting factor that minimizes the error between the full system model and the proposedsecond-order model.\u00E2\u0080\u00A2 The model proposed in Chapter 2 ignores the network losses in order to simplify thesystem. In future work, we can incorporate losses into the model to improve modelaccuracy.\u00E2\u0080\u00A2 The change detection framework in Chapter 4 focuses on load-change disturbances only.However, generator and transmission-line outages are also severe problems in the system.Hence, a compelling direction for future work, is to detect and identify other contingenciesusing a similar method.47Bibliography[1] S. S. Guggilam, C. Zhao, E. Dall\u00E2\u0080\u0099Anese, Y. C. Chen, and S. V. Dhople, \u00E2\u0080\u009COptimizing derparticipation in inertial and primary-frequency response,\u00E2\u0080\u009D IEEE Transactions on PowerSystems, vol. 33, no. 5, pp. 5194\u00E2\u0080\u00935205, Sep. 2018. \u00E2\u0086\u0092 pages xi, 4, 25, 30[2] A. Ipakchi and F. Albuyeh, \u00E2\u0080\u009CGrid of the future,\u00E2\u0080\u009D IEEE power and energy magazine,vol. 7, no. 2, pp. 52\u00E2\u0080\u009362, Mar 2009. \u00E2\u0086\u0092 page 1[3] A. Woyte, V. Van Thong, R. Belmans, and J. Nijs, \u00E2\u0080\u009CVoltage fluctuations on distributionlevel introduced by photovoltaic systems,\u00E2\u0080\u009D IEEE Transactions on energy conversion,vol. 21, no. 1, pp. 202\u00E2\u0080\u0093209, Mar 2006. \u00E2\u0086\u0092 page 1[4] A. Kulmala, S. Repo, and P. Ja\u00C2\u00A8rventausta, \u00E2\u0080\u009CCoordinated voltage control in distributionnetworks including several distributed energy resources,\u00E2\u0080\u009D IEEE Transactions on SmartGrid, vol. 5, no. 4, pp. 2010\u00E2\u0080\u00932020, Jul 2014. \u00E2\u0086\u0092 page 1[5] H. Farhangi, \u00E2\u0080\u009CThe path of the smart grid,\u00E2\u0080\u009D IEEE power and energy magazine, vol. 8,no. 1, Jan 2010. \u00E2\u0086\u0092 page 1[6] A. Chakrabortty and J. H. Chow, Measurement-Based Methods for Model Reduction ofPower Systems Using Synchrophasors. New York, NY: Springer New York, 2013, pp.159\u00E2\u0080\u0093197. \u00E2\u0086\u0092 pages 2, 15, 16[7] B. Stott, J. Jardim, and O. Alsac, \u00E2\u0080\u009CDC power flow revisited,\u00E2\u0080\u009D IEEE Transactions onPower Systems, vol. 24, no. 3, pp. 1290 \u00E2\u0080\u0093 1300, 2009. \u00E2\u0086\u0092 page 2[8] A. J. Germond and R. Podmore, \u00E2\u0080\u009CDynamic aggregation of generating unit models,\u00E2\u0080\u009DIEEE Transactions on Power Apparatus and Systems, vol. PAS-97, no. 4, pp. 1060\u00E2\u0080\u00931069,Jul 1978. \u00E2\u0086\u0092 pages 2, 3[9] J. H. Chow, J. R. Winkelman, M. A. Pai, and P. W. Sauer, \u00E2\u0080\u009CSingular perturbationanalysis of large-scale power systems,\u00E2\u0080\u009D International Journal of Electrical Power &Energy Systems, vol. 12, no. 2, pp. 117\u00E2\u0080\u0093126, Apr 1990. \u00E2\u0086\u0092 pages 2, 3[10] J. Zaborszky, K.-W. Whang, G. Huang, L.-J. Chiang, and S.-Y. Lin, \u00E2\u0080\u009CA clustereddynamic model for a class of linear autonomous systems unsing simple enumerativesorting,\u00E2\u0080\u009D IEEE Transactions on Circuits and Systems, vol. 29, no. 11, pp. 747\u00E2\u0080\u0093758, Nov1982. \u00E2\u0086\u0092 page 3[11] J. M. Undrill and A. E. Turner, \u00E2\u0080\u009CConstruction of power system electromechanicalequivalents by modal analysis,\u00E2\u0080\u009D IEEE Transactions on Power Apparatus and Systems,vol. PAS-90, pp. 2049\u00E2\u0080\u00932059, 1971. \u00E2\u0086\u0092 page 348[12] I. J. Perez-arriaga, G. C. Verghese, and F. C. Schweppe, \u00E2\u0080\u009CSelective modal analysis withapplications to electric power systems, part I: Heuristic introduction,\u00E2\u0080\u009D IEEETransactions on Power Apparatus and Systems, vol. PAS-101, no. 9, pp. 3117\u00E2\u0080\u00933125, Sep1982. \u00E2\u0086\u0092 page 3[13] D. Chaniotis and M. A. Pai, \u00E2\u0080\u009CModel reduction in power systems using Krylov subspacemethods,\u00E2\u0080\u009D IEEE Transactions on Power Systems, vol. 20, no. 2, pp. 888\u00E2\u0080\u0093894, May 2005.\u00E2\u0086\u0092 page 3[14] D. Apostolopoulou, P. W. Sauer, and A. D. Dom\u00C4\u00B1\u00C2\u00B4nguez-Garc\u00C2\u00B4\u00C4\u00B1a, \u00E2\u0080\u009CBalancing authorityarea model and its application to the design of adaptive AGC systems,\u00E2\u0080\u009D IEEETransactions on Power Systems, vol. 31, no. 5, pp. 3756\u00E2\u0080\u00933764, Sep 2016. \u00E2\u0086\u0092 pages 3, 17[15] S. S. Guggilam, C. Zhao, E. Dall\u00E2\u0080\u0099Anese, Y. C. Chen, and S. V. Dhople, \u00E2\u0080\u009CEngineeringinertial and primary-frequency response for distributed energy resources,\u00E2\u0080\u009D in 2017 IEEE56th Annual Conference on Decision and Control (CDC), Dec 2017, pp. 5112\u00E2\u0080\u00935118. \u00E2\u0086\u0092pages 3, 17[16] A. Al-Digs, S. V. Dhople, and Y. C. Chen, \u00E2\u0080\u009CTime-varying injection shift factors topredict post-contingency dynamic line flows,\u00E2\u0080\u009D in Allerton Conference on Communication,Control, and Computing, Oct 2017, pp. 302\u00E2\u0080\u0093306. \u00E2\u0086\u0092 page 3[17] H. Alatrash, A. Mensah, E. Mark, G. Haddad, and J. Enslin, \u00E2\u0080\u009CGenerator emulationcontrols for photovoltaic inverters,\u00E2\u0080\u009D IEEE Transactions on Smart Grid, vol. 3, no. 2, pp.996\u00E2\u0080\u00931011, June 2012. \u00E2\u0086\u0092 page 3[18] A. Anzalchi, M. M. Pour, and A. Sarwat, \u00E2\u0080\u009CA combinatorial approach for addressingintermittency and providing inertial response in a grid-connected photovoltaic system,\u00E2\u0080\u009Din 2016 IEEE Power and Energy Society General Meeting (PESGM), July 2016, pp. 1\u00E2\u0080\u00935.\u00E2\u0086\u0092 page 3[19] S. Ghosh, S. Kamalasadan, N. Senroy, and J. Enslin, \u00E2\u0080\u009CDoubly fed induction generator(dfig)-based wind farm control framework for primary frequency and inertial responseapplication,\u00E2\u0080\u009D IEEE Transactions on Power Systems, vol. 31, no. 3, pp. 1861\u00E2\u0080\u00931871, May2016. \u00E2\u0086\u0092 page 3[20] J. Zhao, X. Lyu, Y. Fu, X. Hu, and F. Li, \u00E2\u0080\u009CCoordinated microgrid frequency regulationbased on dfig variable coefficient using virtual inertia and primary frequency control,\u00E2\u0080\u009DIEEE Transactions on Energy Conversion, vol. 31, no. 3, pp. 833\u00E2\u0080\u0093845, Sep. 2016.[21] J. Van de Vyver, J. D. M. De Kooning, B. Meersman, L. Vandevelde, and T. L.Vandoorn, \u00E2\u0080\u009CDroop control as an alternative inertial response strategy for the syntheticinertia on wind turbines,\u00E2\u0080\u009D IEEE Transactions on Power Systems, vol. 31, no. 2, pp.1129\u00E2\u0080\u00931138, March 2016. \u00E2\u0086\u0092 page 3[22] S. Dong, Y. Chi, and Y. Li, \u00E2\u0080\u009CActive voltage feedback control for hybrid multiterminalhvdc system adopting improved synchronverters,\u00E2\u0080\u009D IEEE Trans. Power Del., vol. 31,no. 2, pp. 445\u00E2\u0080\u0093455, Apr. 2016. \u00E2\u0086\u0092 page 3[23] S. Dong and Y. C. Chen, \u00E2\u0080\u009CAdjusting synchronverter dynamic response speed via dampingcorrection loop,\u00E2\u0080\u009D IEEE Trans. Energy Convers., vol. 32, no. 2, pp. 608\u00E2\u0080\u0093619, Jun. 2017.49[24] S. Dong and Y. C. Chen, \u00E2\u0080\u009CA method to directly compute synchronverter parameters fordesired dynamic response,\u00E2\u0080\u009D IEEE Trans. Energy Convers., vol. 33, no. 2, pp. 814\u00E2\u0080\u0093825,Jun. 2018. \u00E2\u0086\u0092 page 3[25] H. Cha\u00C2\u00B4vez, R. Baldick, and S. Sharma, \u00E2\u0080\u009CGovernor rate-constrained opf for primaryfrequency control adequacy,\u00E2\u0080\u009D IEEE Transactions on Power Systems, vol. 29, no. 3, pp.1473\u00E2\u0080\u00931480, May 2014. \u00E2\u0086\u0092 page 4[26] F. Teng, M. Aunedi, D. Pudjianto, and G. Strbac, \u00E2\u0080\u009CBenefits of demand-side response inproviding frequency response service in the future GB power system,\u00E2\u0080\u009D Frontiers inEnergy Research, vol. 3, p. 36, 2015.[27] F. Teng, V. Trovato, and G. Strbac, \u00E2\u0080\u009CStochastic scheduling with inertia-dependent fastfrequency response requirements,\u00E2\u0080\u009D IEEE Transactions on Power Systems, vol. 31, no. 2,pp. 1557\u00E2\u0080\u00931566, March 2016. \u00E2\u0086\u0092 page 4[28] F. Katiraei, M. Iravani, and P. Lehn, \u00E2\u0080\u009CSmall-signal dynamic model of a micro-gridincluding conventional and electronically interfaced distributed resources,\u00E2\u0080\u009D IETgeneration, transmission & distribution, vol. 1, no. 3, pp. 369\u00E2\u0080\u0093378, 2007. \u00E2\u0086\u0092 page 4[29] T. S. Borsche, T. Liu, and D. J. Hill, \u00E2\u0080\u009CEffects of rotational inertia on power systemdamping and frequency transients,\u00E2\u0080\u009D in 2015 54th IEEE Conference on Decision andControl (CDC), Dec 2015, pp. 5940\u00E2\u0080\u00935946. \u00E2\u0086\u0092 page 4[30] S. S. Guggilam, C. Zhao, E. Dall\u00E2\u0080\u0099Anese, Y. C. Chen, and S. V. Dhople, \u00E2\u0080\u009CEngineeringinertial and primary-frequency response for distributed energy resources,\u00E2\u0080\u009D in 2017 IEEE56th Annual Conference on Decision and Control (CDC), Dec 2017, pp. 5112\u00E2\u0080\u00935118. \u00E2\u0086\u0092page 4[31] Y. C. Chen, T. Banerjee, A. D. Dom\u00C4\u00B1\u00C2\u00B4nguez-Garc\u00C2\u00B4\u00C4\u00B1a, and V. V. Veeravalli, \u00E2\u0080\u009CQuickest lineoutage detection and identification,\u00E2\u0080\u009D IEEE Transactions on Power Systems, vol. 31,no. 1, pp. 749\u00E2\u0080\u0093758, Jan 2016. \u00E2\u0086\u0092 page 4[32] R. Emami and A. Abur, \u00E2\u0080\u009CExternal system line outage identification using phasormeasurement units,\u00E2\u0080\u009D IEEE Transactions on Power Systems, vol. 28, no. 2, pp.1035\u00E2\u0080\u00931040, 2013.[33] H. Zhu and G. B. Giannakis, \u00E2\u0080\u009CSparse overcomplete representations for efficientidentification of power line outages,\u00E2\u0080\u009D IEEE Transactions on Power Systems, vol. 27,no. 4, pp. 2215 \u00E2\u0080\u0093 2224, 2012.[34] J. E. Tate and T. J. Overbye, \u00E2\u0080\u009CLine outage detection using phasor angle measurements,\u00E2\u0080\u009DIEEE Transactions on Power Systems, vol. 23, no. 4, pp. 1644 \u00E2\u0080\u0093 1652, 2008. \u00E2\u0086\u0092 page 4[35] Y. Shin, E. J. Powers, M. Grady, and A. Arapostathis, \u00E2\u0080\u009CPower quality indices fortransient disturbances,\u00E2\u0080\u009D IEEE Transactions on Power Delivery, vol. 21, no. 1, pp.253\u00E2\u0080\u0093261, Jan 2006. \u00E2\u0086\u0092 page 4[36] W. G. Morsi and M. E. El-Hawary, \u00E2\u0080\u009CFuzzy-wavelet-based electric power qualityassessment of distribution systems under stationary and nonstationary disturbances,\u00E2\u0080\u009DIEEE Transactions on Power Delivery, vol. 24, no. 4, pp. 2099\u00E2\u0080\u00932106, Oct 2009.50[37] R. H. G. Tan and V. K. Ramachandaramurthy, \u00E2\u0080\u009CVoltage sag acceptability assessmentusing multiple magnitude-duration function,\u00E2\u0080\u009D IEEE Transactions on Power Delivery,vol. 27, no. 4, pp. 1984\u00E2\u0080\u00931990, Oct 2012. \u00E2\u0086\u0092 page 4[38] L. Xie, Y. Chen, and P. R. Kumar, \u00E2\u0080\u009CDimensionality reduction of synchrophasor data forearly event detection: linearized analysis,\u00E2\u0080\u009D IEEE Transactions on Power Systems, vol. 29,no. 6, pp. 2784\u00E2\u0080\u00932794, Nov 2014. \u00E2\u0086\u0092 page 4[39] I. Kamwa, A. K. Pradhan, and G. Joos, \u00E2\u0080\u009CRobust detection and analysis of power systemoscillations using the Teager-Kaiser energy operator,\u00E2\u0080\u009D IEEE Transactions on PowerSystems, vol. 26, no. 1, pp. 323\u00E2\u0080\u0093333, Feb 2011.[40] M. J. Smith and K. Wedeward, \u00E2\u0080\u009CEvent detection and location in electric power systemsusing constrained optimization,\u00E2\u0080\u009D in 2009 IEEE Power Energy Society General Meeting,Jul 2009, pp. 1\u00E2\u0080\u00936.[41] G. Rigatos, P. Siano, and A. Piccolo, \u00E2\u0080\u009CNeural network-based approach for early detectionof cascading events in electric power systems,\u00E2\u0080\u009D IET Generation, TransmissionDistribution, vol. 3, no. 7, pp. 650\u00E2\u0080\u0093665, Jul 2009. \u00E2\u0086\u0092 page 4[42] G. W. Hart, \u00E2\u0080\u009CNonintrusive appliance load monitoring,\u00E2\u0080\u009D Proceedings of the IEEE, vol. 80,no. 12, pp. 1870\u00E2\u0080\u00931891, Dec 1992. \u00E2\u0086\u0092 page 4[43] Y. Su, K. Lian, and H. Chang, \u00E2\u0080\u009CFeature selection of non-intrusive load monitoring systemusing stft and wavelet transform,\u00E2\u0080\u009D in 2011 IEEE 8th International Conference one-Business Engineering, Oct 2011, pp. 293\u00E2\u0080\u0093298.[44] A. F. Ebrahim and O. Mohammed, \u00E2\u0080\u009CHousehold load forecasting based on apre-processing non-intrusive load monitoring techniques,\u00E2\u0080\u009D in 2018 IEEE GreenTechnologies Conference, Apr 2018, pp. 107\u00E2\u0080\u0093114. \u00E2\u0086\u0092 page 4[45] A. G. Phadke and B. Kasztenny, \u00E2\u0080\u009CSynchronized phasor and frequency measurementunder transient conditions,\u00E2\u0080\u009D IEEE Transactions on Power Delivery, vol. 24, no. 1, pp.89\u00E2\u0080\u009395, Jan 2009. \u00E2\u0086\u0092 page 4[46] F. Hu, K. Sun, A. Del Rosso, E. Farantatos, and N. Bhatt, \u00E2\u0080\u009CMeasurement-basedreal-time voltage stability monitoring for load areas,\u00E2\u0080\u009D IEEE Transactions on PowerSystems, vol. 31, no. 4, pp. 2787\u00E2\u0080\u00932798, July 2016. \u00E2\u0086\u0092 page 4[47] K. Vu, M. M. Begovic, D. Novosel, and M. M. Saha, \u00E2\u0080\u009CUse of local measurements toestimate voltage-stability margin,\u00E2\u0080\u009D IEEE Transactions on Power Systems, vol. 14, no. 3,pp. 1029\u00E2\u0080\u00931035, Aug 1999. \u00E2\u0086\u0092 page 5[48] I. Smon, G. Verbic, and F. Gubina, \u00E2\u0080\u009CLocal voltage-stability index using tellegen\u00E2\u0080\u0099stheorem,\u00E2\u0080\u009D IEEE Transactions on Power Systems, vol. 21, no. 3, pp. 1267\u00E2\u0080\u00931275, Aug2006. \u00E2\u0086\u0092 page 5[49] D.-J. Lee and L. Wang, \u00E2\u0080\u009CSmall-signal stability analysis of an autonomous hybridrenewable energy power generation/energy storage system part i: Time-domainsimulations,\u00E2\u0080\u009D IEEE Transactions on Energy Conversion, vol. 23, no. 1, pp. 311\u00E2\u0080\u0093320,2008. \u00E2\u0086\u0092 page 551[50] P. W. Sauer and M. A. Pai, Power System Dynamics and Stability. Upper Saddle River,NJ: Prentice-Hall, Inc., 1998. \u00E2\u0086\u0092 pages 6, 8, 10[51] B. Chen, A. Al-Digs, and Y. C. Chen, \u00E2\u0080\u009CA network-cognizant aggregate-frequencyreduced-order power system dynamical model,\u00E2\u0080\u009D in 2018 North American PowerSymposium (NAPS), Sep. 2018, pp. 1\u00E2\u0080\u00936. \u00E2\u0086\u0092 page 6[52] P. W. Sauer and B. Hoveida, \u00E2\u0080\u009CConstrained stochastic power flow analysis,\u00E2\u0080\u009D ElectricPower Systems Research, vol. 5, no. 2, pp. 87 \u00E2\u0080\u0093 95, 1982.[53] R. A. Horn and C. R. Johnson, Matrix Analysis. New York, NY: Cambridge UniversityPress, 2013. \u00E2\u0086\u0092 page 11[54] M. D. Ilic\u00C2\u00B4 and Q. Liu, Toward Sensing, Communications and Control Architectures forFrequency Regulation in Systems with Highly Variable Resources. New York, NY:Springer New York, 2012, pp. 3\u00E2\u0080\u009333. \u00E2\u0086\u0092 page 16[55] R. Ramanujam, Power System Dynamics: Analysis and Simulation. PHI Learning Pvt.Ltd., 2009. \u00E2\u0086\u0092 page 17[56] P. M. Anderson and M. Mirheydar, \u00E2\u0080\u009CA low-order system frequency response model,\u00E2\u0080\u009DIEEE Transactions on Power Systems, vol. 5, no. 3, pp. 720\u00E2\u0080\u0093729, 1990. \u00E2\u0086\u0092 page 17[57] F. Milano, \u00E2\u0080\u009CAn open source power system analysis toolbox,\u00E2\u0080\u009D IEEE Trans. on PowerSystems, vol. 20, no. 3, pp. 1199\u00E2\u0080\u00931206, Aug 2005. \u00E2\u0086\u0092 pages 18, 34, 44[58] N. P. C. Council, \u00E2\u0080\u009CNPCC regional reliability reference directory# 12 under-frequencyload shedding program requirements,\u00E2\u0080\u009D 2009. \u00E2\u0086\u0092 page 33[59] H.-Q. TransE\u00C2\u00B4nergie, \u00E2\u0080\u009CTransmission provider technical requirements for the connection ofpower plants to the hydro-que\u00C2\u00B4bec transmission system,\u00E2\u0080\u009D 2009. \u00E2\u0086\u0092 page 33[60] M. Hyder and K. Mahata, \u00E2\u0080\u009CAn approximate L0 norm minimization algorithm forcompressed sensing,\u00E2\u0080\u009D in 2009 IEEE International Conference on Acoustics, Speech andSignal Processing, Apr 2009, pp. 3365\u00E2\u0080\u00933368. \u00E2\u0086\u0092 page 43[61] J. F. C. Mota, J. M. F. Xavier, P. M. Q. Aguiar, and M. Puschel, \u00E2\u0080\u009CDistributed basispursuit,\u00E2\u0080\u009D IEEE Transactions on Signal Processing, vol. 60, no. 4, pp. 1942\u00E2\u0080\u00931956, Apr2012. \u00E2\u0086\u0092 page 43[62] J. F. C. Mota, J. M. F. Xavier, P. M. Q. Aguiar, and M. Pu\u00C2\u00A8schel, \u00E2\u0080\u009CD-ADMM: Acommunication-efficient distributed algorithm for separable optimization,\u00E2\u0080\u009D IEEETransactions on Signal Processing, vol. 61, no. 10, pp. 2718\u00E2\u0080\u00932723, May 2013. \u00E2\u0086\u0092 page 43[63] E. van den Berg and M. P. Friedlander, \u00E2\u0080\u009CSPGL1: A solver for large-scale sparsereconstruction,\u00E2\u0080\u009D Jun 2007, http://www.cs.ubc.ca/labs/scl/spgl1. \u00E2\u0086\u0092 pages 43, 45[64] E. van den Berg and M. P. Friedlander, \u00E2\u0080\u009CProbing the pareto frontier for basis pursuitsolutions,\u00E2\u0080\u009D SIAM Journal on Scientific Computing, vol. 31, no. 2, pp. 890\u00E2\u0080\u0093912, 2008. \u00E2\u0086\u0092pages 43, 4552[65] US DOE & FERC. (2006, Feb) Steps to establish a real-time transmission monitoringsystem for transmission owners and operators within the eastern and westerninterconnections. [Online]. Available:http://energy.gov/sites/prod/files/oeprod/DocumentsandMedia/\final 1839.pdf \u00E2\u0086\u0092 page4453Appendix AModel DetailsA.1 Detailing BGG, BGL, BLG, BLL MatricesBased on the structures of (2.9) and (2.10), matrices BGG , BGL, BLG , and BLL can be expressedasBGG = BGG \u00E2\u0088\u0092 diag([BGG BGL]1N +1GX \u00E2\u0080\u00B2d), (A.1)BGL = BGL, (A.2)BLG = BLG, (A.3)BLL = BLL \u00E2\u0088\u0092 diag([BLG BLL]1N), (A.4)where submatrices BGG, BGL, BLG, and BLL are extracted from the standard (and appropri-ately reordered) network admittance matrix Y for a lossless system, i.e.,Y = j\u00EF\u00A3\u00AE\u00EF\u00A3\u00AF\u00EF\u00A3\u00B0BGG BGLBLG BLL\u00EF\u00A3\u00B9\u00EF\u00A3\u00BA\u00EF\u00A3\u00BB . (A.5)54A.2 Showing H1G = 0GSubstitute (2.14) and (2.18) into (2.22) to getH = \u00E2\u0088\u0092(diag(1GX \u00E2\u0080\u00B2d)+ J)J\u00E2\u0088\u00921diag(1GX \u00E2\u0080\u00B2d)= \u00E2\u0088\u0092diag(1GX \u00E2\u0080\u00B2d)(J\u00E2\u0088\u00921diag(1GX \u00E2\u0080\u00B2d)+ diag(1G)), (A.6)where J = BGG \u00E2\u0088\u0092BGLB\u00E2\u0088\u00921LLBLG . Below, we proceed to show that H1G = 0G by proving that0G =(J\u00E2\u0088\u00921diag(1GX \u00E2\u0080\u00B2d)+ diag(1G))1G, (A.7)which is equivalent to showing that0G =(diag(1GX \u00E2\u0080\u00B2d)+ J)1G=1GX \u00E2\u0080\u00B2d+BGG1G \u00E2\u0088\u0092BGLB\u00E2\u0088\u00921LLBLG1G. (A.8)We proceed to prove (A.8) as follows. First, using (A.1), we getBGG1G = BGG1G \u00E2\u0088\u0092 diag([BGG BGL]1N +1GX \u00E2\u0080\u00B2d)1G= BGG1G \u00E2\u0088\u0092[BGG BGL]1N \u00E2\u0088\u0092 1GX \u00E2\u0080\u00B2d= \u00E2\u0088\u0092BGL1L \u00E2\u0088\u0092 1GX \u00E2\u0080\u00B2d= \u00E2\u0088\u0092BGL1L \u00E2\u0088\u0092 1GX \u00E2\u0080\u00B2d. (A.9)where the last equality above follows by substituting (A.2). Next, using (A.4), we have thatBLL1L = BLL1L \u00E2\u0088\u0092 diag([BLG BLL]1N)1G= BLL1L \u00E2\u0088\u0092[BLG BLL]1N= \u00E2\u0088\u0092BLG1G = \u00E2\u0088\u0092BLG1G, (A.10)55where the last equality above follows by substituting (A.3). Finally, we substitute (A.9)and (A.10) into the right-hand side of (A.8) to obtain1GX \u00E2\u0080\u00B2d\u00E2\u0088\u0092BGL1L \u00E2\u0088\u0092 1GX \u00E2\u0080\u00B2d+BGLB\u00E2\u0088\u00921LLBLL1L = 0G, (A.11)and H1G = 0G, as desired.56Appendix BParameters for Numerical CaseStudiesSynchronous frequency, \u00CF\u0089s = 2pi60[rad sec\u00E2\u0088\u00921].Parameters modified from the standard IEEE WECC 9-bus test system: Generator dampingconstants D1 = D2 = 5, droop constants R1 = R2 = 0.04, inertia M1 = 13.64, M2 = 6.4[sec];DER damping constant D3 = 5, inertia M3 = 3.01[sec].Parameters modified from the standard IEEE New England 39-bus test system: Generatordamping constants D1 = \u00C2\u00B7 \u00C2\u00B7 \u00C2\u00B7 = D7 = 2, droop constants R1 = \u00C2\u00B7 \u00C2\u00B7 \u00C2\u00B7 = R7 = 0.04, inertiaM1 = 4.6372, M2 = \u00C2\u00B7 \u00C2\u00B7 \u00C2\u00B7 = M7 = 5.2838[sec].57"@en .
"Thesis/Dissertation"@en .
"2019-09"@en .
"10.14288/1.0380452"@en .
"eng"@en .
"Electrical and Computer Engineering"@en .
"Vancouver : University of British Columbia Library"@en .
"University of British Columbia"@en .
"Attribution-NonCommercial-NoDerivatives 4.0 International"@* .
"http://creativecommons.org/licenses/by-nc-nd/4.0/"@* .
"Graduate"@en .
"A reduced-order power-system dynamical model and applications in design and monitoring of power systems"@en .
"Text"@en .
"http://hdl.handle.net/2429/71307"@en .