Dynamic Modeling, Stability, and Control of Power Systems with Distributed Energy Resources
I Introduction
Significant infrastructural changes are currently being implemented on power system networks around the world by maximizing the penetration of renewable energy, by installing new transmission lines, by adding flexible loads, by promoting independence in power production by disintegrating the grid into microgrids, and so on [1]. The shift of energy supply from large central generating stations to smaller producers such as wind farms, solar photovoltaic (PV) farms, rooftop PVs, and energy storage systems, collectively known as distributed energy resources (DERs) or inverterinterfaced resources (IBRs), is accelerating at a very rapid pace. Hundreds of power electronic devices are being added, creating hundreds of new control points in the grid. This addition is complimented by an equal progress in sensing technology, whereby highprecision, high samplingrate, GPSsynchronized dynamic measurements of voltages and currents are now available from sensors such as Phasor Measurement Units (PMUs) [2]. With all these transformational changes happening in the grid, operators are inclining to explore new control methods that go far beyond how the grid is controlled today.
In the current stateofart, power system controllers, especially the ones that are responsible for transient stability and power oscillation damping, are all operated in a decentralized and uncoordinated fashion using local output feedback only. A survey of these controllers is given in Appendix A. With rapid modernization of the grid, these local controllers, however, will not be tenable over the longterm. Instead, systemwide coordinated controllers will become essential. Such controllers where signals measured at one part of the network are communicated to other remote parts for feedback are called widearea controllers [3].
Widearea control alone, however, will not be enough either. It may improve the damping performance of the legacy system, but will not be able to keep up with the unpredictable rate at which DERs are being added to the grid. Every time a new DER is added it will be almost impossible for an operator to retune all the widearea control gains to accommodate the change in dynamics. DERs have high variability and intermittency, and need to be operated in a plugandplay fashion. Accordingly their controllers need to be local, decentralized, and modular in both design and implementation. In other words, neither should the design of one DER controller depend on that of another nor should either of these two controllers need to be updated when a third DER is added in the future. The overall control architecture for the future grid needs to be a combination of these decentralized plugandplay DER controllers and distributed widearea controllers.
The objective of this article is twofold. The first objective is to present a suite of new control methods for developing this combined control architecture. For brevity, the design will be limited to only one particular application  namely, adding damping to the oscillations in power flows after both small and large disturbances, also called power oscillation damping (POD) in short [4]. POD is one of the most critical realtime control problems in today’s power grid, and its importance is only going to increase with DER integration. The applicability of the design, however, go far beyond just POD to many other grid control problems such as frequency control, voltage control, and congestion relief. Efficacy of the methods for enhancing transient stability as a bonus application will be illustrated via simulations. While many papers on decentralized DER control and distributed widearea control exist in the literature (surveys on these two control methods are given in Appendix B and Appendix C, respectively), very few have studied the simultaneous use of both. Moreover, most DER controllers reported so far lack the modularity and plugandplay characteristic explained above. The control methods presented in this article address all of these challenges. The second objective is to present a comprehensive list of mathematical models of the various components of a power grid ranging from synchronous generators, their internal controllers, loads, wind and solar farms, batteries to the power electronic device interfaces and associated control mechanisms for each of these components. While these models may be individually wellcited in the literature, very few references so far have collected all of them together to understand the holistic dynamic behavior of an entire grid.
The rest of the article is organized as follows. Section II describes the dynamic model of a power system, integrated with different types of DERs. A general framework for modeling is provided first, followed by details of each individual component model. Section III demonstrates the impacts of DERs on power system dynamics through numerical simulations. Motivated by these simulation results, Section IV develops new decentralized DER control laws using the idea of retrofit control [5], as well as distributed widearea controllers for damping of low frequency oscillations using sparse optimal control. The effectiveness of this combined control strategy is demonstrated on the IEEE 68bus power system with wind and solar farms. The article concludes with a list of open research problems.
Ii Power System Models
Symbol  Numerical value  Description 

base angular speed for a 60Hz power system. Unit of is rad/sec.  
number of buses  
th bus voltage  
,  active and reactive power injected from the th component  
state of the th component  
control input of the th component  
model parameter depending on operating point  
admittance matrix  
, , , ,  index set of the buses connecting to generators, loads, wind farms, solar farms, energy storages, and that of nonunit buses. These sets are disjoint, and 
First, the dynamic models of the four core components of a power system are developed  namely, generation, transmission, load, and energy storage. The generating units are classified into conventional power plants and DERs such as wind generators and PV generators. Each model follows from firstprinciples of physics. Note that in reality a generation facility, whether that be conventional generation or wind/solar generation, and energy storage facilities contain many generating units and storage devices inside them. In the following, the terms generator, wind farm, solar farm, and energy storage system are used to refer to an aggregate of those individual units representing the overall facilities. Similarly, the term load is used to refer to an aggregate of all consumers inside the associated demand area. Each aggregated unit comes with its own individual bus, such as a generator bus or a load bus. The buses are connected to each other through a network of power transmission lines. The power system may also contain buses where no generator, wind/solar farm, load, or energy storage system is connected. These buses are called nonunit buses. The term component is used to refer to either a unit with its bus or the nonunit bus. An example of these connected components is shown in Fig. 1.
As will be shown in the following, a general form for the dynamic model of the th component of a power system, whether that component be a generator, load, storage, wind farm, or solar farm, can be written as
(1) 
for , where the nomenclature of this model is summarized in Table I. Details of the two functions and for each component of the grid will be described shortly. Throughout the article, complex variables will be written in bold fonts, (for example ). All symbols with superscript will denote setpoints.
The components are interconnected by a transmission network. Let denote the admittance matrix of the network (for details of the construction of this matrix, please see “Construction of Admittance Matrices”). The power balance across the transmission lines follows from Kirchoff’s laws as
(2) 
where is the elementwise multiplication, is the elementwise complex conjugate operator, , , and are the stacked representations of , and for . From (2), is determined for a given and . The overall dynamics of a power system can be described by the combination of (1) and (2). A signalflow diagram of this model is shown in Fig. 2.
The power system model (1)(2) is operated around its equilibrium. This is determined as follows. The steadystate value of , , , and , and parameter in (1) must satisfy
(3) 
and (2). The steadystate value of is assumed to be zero without loss of generality. A standard procedure of finding the steadystate values consists of two steps: first, find for satisfying (2), and then, given the triple find the pair satisfying (3) for each . The first step is called power flow calculation, details of which are given in “Brief Tutorial on Power Flow Calculation”. The second step, often called an initialization step, will be described later in this section.
In the following, the statespace models of generators, loads, energy storage systems, wind farms, solar farms, and nonunit buses, conforming to the structure in (1), are derived. For easier understanding, each subsection starts with a qualitative description of the respective component model followed by its statespace representation. Some parts may refer to equations that appear later in the text. To simplify the notation, the subscript is omitted unless otherwise stated.
Iia Generators
A generator consists of a synchronous machine, an energy supply system (or a primemover), and an excitation system [4]. The excitation system induces currents in the excitation winding, and thereby magnetizes the rotor. The primemover generates mechanical power to rotate the rotor in this magnetic field. The synchronous machine converts the mechanical power to electrical power, which is transmitted to the rest of the grid. The dynamics of the primemover is usually ignored because of its slow timescale.
IiA1 Synchronous Machine
While various types of synchronous machine models are available in the literature (for example see [4]), in this article a wellknown model called the oneaxis model or fluxdecay model is used. This model consists of the electromechanical swing dynamics (4) and the electromagnetic voltage dynamics (5). For simplicity, the mechanical power in (4) is assumed to be constant.
IiA2 Excitation System
Typically, the excitation system consists of an exciter, an Automatic Voltage Regulator (AVR) that regulates the generator voltage magnitude to its setpoint value, and a Power System Stabilizer (PSS) that ensures smallsignal stability. The exciter with AVR is modeled as (6), where is a control input representing an additional voltage reference signal to the AVR. The PSS is taken as a typical speedfeedback type controller which consists of two stage leadlag compensators and one highpass washout filter [7].
The statespace representation of the overall generator model can be written as follows (definitions of , , , are given in Table I while those of the other symbols are provided in Table II):
Symbol  Numerical value  Description 

rotor angle relative to the frame rotating at . Unit of is (rad)  
frequency deviation, that is, rotor angular velocity relative to  
qaxis voltage behind transient reactance  
field voltage  
PSS state  
exciter field voltage  
PSS output  
additional voltage reference to AVR  
30  inertia constant (sec)  
0.1  damping coefficient  
0.1  daxis transient opencircuit time constant (sec)  
,  1.8  d and qaxis synchronous reactance 
0.3  daxis transient reactance  
8.0  steadystate mechanical power  
0.05  time constant of exciter (sec)  
20  AVR gain  
150  PSS gain  
1.0  setpoints for the field voltage  
1.0  setpoints for the bus voltage magnitude  
10  washout filter time constant (sec)  
,  0.02, 0.07  leadlag time constants of the first stage of PSS (sec) 
,  0.02, 0.07  leadlag time constants of the second stage of PSS (sec) 
Synchronous machine:
Electromechanical swing dynamics:
(4) 
Electromagnetic dynamics:
(5) 
Excitation System:
Exciter with AVR:
(6) 
where represents the setpoint of .
PSS:
(7) 
where
(8) 
IiB Nonunit Buses
Nonunit buses are simply modeled by the Kirchoff’s power balance law, namely for ,
(10) 
In reference to (1) this means , , and are empty vectors.
IiC Loads
Loads are commonly modeled by algebraic power balance equations, although extensive literature also exists for dynamic loads (for example, see [8, 9]). The wellknown static load models are:
(11)  
(12)  
(13) 
where is the complex conjugate operator, , , and are constant. Therefore, for , in reference to (1) a load at the th bus can be represented by and being empty vectors, being either , , or , and the output equation being either (11), (12) or (13). For the simulations later in this article, constant impedance loads will be used. For a given triple , the load impedance will be calculated such that .
IiD Wind Farms
A wind farm model typically consists of a wind turbine, a doublyfed induction generator (DFIG), and a backtoback (B2B) converter with associated controllers. A battery with DC/DC converter can be added to the B2B converter if needed. Fig. 4 shows the physical architecture of a wind farm with its bus while Fig. 5 shows a signalflow diagram of the model. When the battery is not connected, the current in Fig. 4 is regarded as zero. The symbols for the wind farm model are listed in Table III.
IiD1 Wind Turbine
The wind turbine, as shown in Fig. 4, converts aerodynamic power coming from the wind to mechanical power that is transmitted to the DFIG. The turbine is typically modeled as a oneinertia or twoinertia model (the latter is followed in this article) consisting of a lowspeed shaft, a highspeed shaft, and a gearbox [10]. For simplicity, the aerodynamic power in the twoinertia dynamics (15) is assumed to be constant as wind speeds usually change slowly.
IiD2 Dfig
The DFIG shown in Fig. 4 converts the mechanical power from the turbine into electrical power. The DFIG consists of a threephase rotor and a threephase stator. The stator is connected to the wind bus to transmit the electrical power into the grid while the rotor is connected to the B2B converter with associated controllers that control the rotor winding voltage. The stator and rotor are coupled electromagnetically, which is reflected in the dynamics of the stator and rotor currents expressed in a rotating dq reference frame [11].
IiD3 B2B Converter with its Controllers
The B2B converter is used for regulating the DFIG rotor voltages , as well as the reactive power flowing from the stator to the converter. The B2B converter consists of two threephase voltage source converters, namely, the rotorside converter (RSC) and the gridside converter (GSC), linked via a common DC line [12]. Each of the converters is equipped with a controller. The explanation of the models of GSC, RSC, and their controllers is as follows.
Following [12], the GSC dynamics is expressed as the variation of the ACside current in dq reference frame.
The GSC controller consists of an innerloop controller and an outerloop controller [12]. The objective of the outerloop controller is to generate a reference signal of the GSC currents for regulating both of the DC link voltage and the reactive power flowing into the GSC to their respective setpoints. The outerloop controller is designed as a PI controller as (19). The innerloop controller aims at regulating to the generated reference signals , by the control of the duty cycles , . Following [12], in this article the controller is designed such that the transfer function from (or ) to (or ) is a desired firstorder system when the duty cycles are not saturated. The controller is implemented as (20).
The RSC model is described as
(14) 
where and are the DFIG rotor currents, and are the duty cycles of the RSC, and is the DC link voltage. In this article, the RSC resistance and inductance are considered to be negligible, that is, . This assumption is always satisfied by incorporating the two into the DFIG rotor circuit. Thus, the RSC model used in this article is described as (21).
The RSC is equipped with an innerloop controller and an outerloop controller. The outerloop controller generates reference signals for the DFIG rotor currents and for regulating the stator voltage magnitude and the highspeed shaft speed to their setpoints while the innerloop controller aims at regulating the RSC currents [13]. This control action is actuated through the control of the duty cycles of the B2B converter.
The RSC and GSC are connected by a DC link equipped with a capacitor whose dynamics is derived from the power balance through the B2B converter [12].
IiD4 Battery and DC/DC converter
A battery is used for charging or discharging of electricity whenever needed. The battery comes with a DC/DC converter that steps up/down the battery terminal voltage. Both devices are also sometimes used for suppressing the fluctuations in the output power by controlling the DC/DC converter. The model for each is described as follows.
The DC/DC converter is modeled by buck (stepdown) and boost (stepup) models. These models are widely available in the literature [14]. When the converter dynamics are sufficiently fast, simpler models where the output voltage and current are explicit functions of the duty ratio can be derived. In this article this simple model is used.
IiD5 Interconnection to Grid
The net active and reactive power injected by the wind farm to the grid are determined as the sum of the power leaving from the stator and that consumed by the B2B converter.
Symbol  Numerical value  Description 
Wind turbine  
,  angular velocity of lowspeed shaft and highspeed shaft  
torsion angle (rad)  
aerodynamic power input depending on wind speed  
,  inertia coefficients of the lowspeed and highspeed shafts (sec)  
,  friction coefficients of the lowspeed and highspeed shafts  
torsional stiffness (1/rad)  
damping coefficient of turbine  
90  gear ratio  
mechanical synchronous frequency (rad/sec)  
DFIG  
,  d and qaxis rotor currents  
,  d and qaxis stator currents  
,  d and qaxis rotor voltages  
electromechanical torque converted by DFIG  
power flowing from DFIG to bus  
number of wind generators inside the farm  
magnetizing reactance  
,  stator and rotor leakage reactance  
,  stator and rotor resistance  
GSC and its controller  
,  d and qaxis currents flowing from AC side to DC side  
,  d and qaxis duty cycles  
power flowing from bus to GSC  
steadystate value of  
steadystate DC link voltage  
,  innerloop controller state  
,  outer controller state  
,  reference signal of and  
,  additional control input signals  
,  inductance and resistance of GSC  
,  P gains of outer controller  
,  I gains of outer controller  
GSC current dynamics’ time constant to be designed  
RSC and its controller  
,  d and qaxis duty cycles  
,  innerloop controller state  
,  reference signal of and  
,  additional control input signals  
,  P gains of outer controller  
,  P gains of the innerloop controller  
,  I gains of the innerloop controller  
DC link  
DCside voltage  
DC link capacitance  
conductance representing switching loss of the B2B converter  
BuckandBoost DC/DC converter  
current injecting from DC/DC converter into DC link  
voltage at battery side  
step down/up gain  
Battery  
battery voltage  
current injected from the battery  
battery capacity  
battery conductance  
,  resistance and inductance of battery circuit 
InnerLoop controller of GSC:
(20) 
where and are defined in (18), and in (19), in (24), and is a saturation function whose output is restricted within the range of .
DC link:
(24) 
where and are defined in (18), and in (21), and in (16), and in (25). When the battery and DC/DC are not connected, .
In reference to (1), the wind farm model with the battery and DC/DC converter can be summarized as:

(28) 
and and in (1) follow from (15)(27) for . The steadystate is determined as follows. Note that, given the total generated power , there exists a degree of freedom for determining and satisfying (27) in steady state. Thus, not only the triple but also the pair needs to be known. In this setting, the pair satisfying (3) is uniquely determined.
IiE Solar Farms
Symbol  Numerical value  Description 
DC/AC converter and its controller  
,  d and qaxis currents flowing from ACside to DCside  
,  d and qaxis duty cycles  
power injecting from solar bus  
,  innerloop controller state  
,  outerloop controller state  
,  reference signal of and  
,  additional control input signals on duty cycles  
number of PV generators inside farm  
,  inductance and resistance of DC/AC converter  
steadystate power injecting from solar bus  
,  PI gains of daxis outerloop controller  
,  PI gains of qaxis outerloop controller  
0.7  design parameter representing time constant of converter current dynamics  
DC link  
DC link voltage  
DC link capacitance  
conductance representing switching loss of DC/AC converter  
BuckandBoost DC/DC converter  
current flowing from DC/DC converter to DC link  
voltage at PV array side  
step down/up gain so that solar farm is operated at MPP  
PV array  
current flowing from the PV array to the DC link  
series resistance inside PV array model  
voltage of constant voltage source inside PV array 
A solar farm model consists of a PV array, a buckandboost DC/DC converter, a DC/AC converter with a controller, and a DC link [15], as shown in Fig. 6. The signalflow diagram for the system is shown in Fig. 7. The dynamics of the DC/AC converter, its controller, and DC link are similar to those in the wind farm model, given in (31)(34). The models of the PV array and DC/DC converter are described as follows.
IiE1 PV array
The PV array is a parallel interconnection of circuits, each of which contains seriesconnected PV cells, as shown in Fig. 8 (a). Each PV cell is assumed to be identical. Typically, a PV cell has nonlinear IV characteristics, as shown by the blue line in Fig. 8 (b) [15]. Assuming that the PV cell is operated around the socalled maximum power point (MPP) where the cell output power is maximized, the IV curve around this point can be approximated by a linear function as shown by the red line. In that case, the PV array can be modeled as a series connection of a constant voltage source with value and a resistance whose value is . This PV array model is described as (29).
IiE2 DC/DC converter
A buckandboost DC/DC converter is used to ensure PV cell operation around MPP. This can be done by determining the converter step down/up gain such that the steadystate PV array output voltage and current coincide with the maximum point on the IV curve. While the gain can be dynamically regulated by controllers, for simplicity, the gain is supposed to be constant. The DC/DC converter model is described as (30).
InnerLoop controller of DC/AC converter:
(33) 
IiF Energy Storage Systems
The energy storage system consists of a battery, a buckandboost DC/DC converter, a DC/AC converter, and a controller, as shown in Fig. 9. The basic functions of these four components are to charge/discharge electricity, to step down/up the battery terminal voltage, to rectify the threephase current to a DC current, and to regulate the DC voltage in between the converters, respectively. When the energy storage system is connected to DC line, the DC/AC converter is not needed. The dynamics of the DC/AC converter, its controller, DC/DC converter, and DC link are similar to those described in equations (31), (32)(33), (30), and (34), respectively.
Iii Impact of DERs on Power System Dynamics
Given a DERintegrated power system model (1)(2), the question is  how does the penetration of DERs and their controllers dictate the stability and dynamic performance of the grid? This section demonstrates these impacts using numerical simulation of the IEEE 68bus power system model [16]. The network diagram is shown in Fig. 10. The DER bus, denoted as Bus 69, connects to Bus 22. The reactance between Bus 22 and 69 represents the transformer for stepping down the grid voltage to the DER voltage. Its value is taken as .
First, consider the DER to be a solar farm as in (29)(34) with . The other bus indices , , are shown in Fig. 10. The model of this PVintegrated power system is the combination of (1)(2) where for , , and are defined as (4)(8), (11), (10), and (29)(34). Note that is the number of PV generators inside the farm. A question here is  how does affect smallsignal stability of the grid model? Fig. 11 shows the 13 dominant eigenvalues of the linearized power system model at a desired equilibrium. The eigenvalues around for start moving to the right as the value of is increased, and finally cross the imaginary axis when , resulting in an unstable system. Each PV generator is rated at 2 MW; therefore means that the net steadystate power output of the solar farm is MW, which is 3.85% of the total generated power of the system. This may look like a small percentage, but in terms of the stability limit the amount of solar penetration is quite close to critical. This pole shift happens due to the fact that the equilibrium changes with . When a fault (modeled as an impulse function causing the initial conditions of and in (31) to move from their equilibrium values) is induced oscillations in the transient response of the states can easily be seen. Fig. 12 shows the frequency deviation of all 16 synchronous generators for the cases where . The results indicate that as increases the PVintegrated power system, without any DER control, becomes oscillatory with poor damping.
Next, the solar farm at Bus 69 is replaced by a wind farm without a battery or a DC/DC converter. Fig. 13 shows the first 14 dominant eigenvalues of the linearized windintegrated power system at a desired equilibrium. The eigenvalues around for start moving to the right as the value of is increased, and finally cross the imaginary axis when , resulting in an unstable system. Each wind generator is rated at 2 MW; therefore, means that the net steadystate power output of the farm is MW. Thus, compared to the PV penetration, the wind penetration in this case poses a greater threat to smallsignal stability. To investigate the difference between the two, the singular value plot of the frequency response of each model from the d and qaxis bus voltages to the injected power is shown in Fig. 14 (a) and (b), respectively. The figure shows that the wind farm model has a resonance peak at 0.157 Hz, and the amplitude of the peak increases as is increased. This is an interesting observation since 0.157 Hz lies in the range of frequencies for the lowfrequency (0.1 Hz to 2 Hz) oscillations of the synchronous generators, commonly called interarea oscillations [17]. The wind injection at Bus 69, thus, stimulates an interarea mode in this case. The resonance mode actually stems from the internal characteristics of the DFIG dynamic model. Details of this phenomenon can be found in [18]. The PV model, on the other hand, does not show any such resonance peak.
One potential way to combat the poorly damped oscillation would be to tune the PI gains of the converter controller (19)(20) and (22)(23). However, such tuning must be done extremely carefully with full knowledge of the entire grid model, since high values of these gains can jeopardize closedloop stability. This is shown in Fig. 15, where the first ten dominant eigenvalues of the linearized closedloop model for are shown. The integrator gains that are more than 13.5 end up destabilizing the power system. This is because highgain controllers stimulate the negative coupling effect between the DER and the rest of the grid. These observations show how many of the power system damping controllers used in today’s grid can easily become invalid tomorrow. Much more systematic control mechanisms need to be built for the future grid to accommodate deep DER penetration while increasing flexibility and robustness.
Iv New Approaches for Control of DERintegrated Power Systems
Iva Local control of DERs
To counteract the destabilizing effects that may be caused by deep penetration of DERs in a power grid, as shown in the previous section, local control mechanism for each individual DER needs to be built. A brief survey of local controllers used in today’s grid, both with and without DERs, is summarized in Appendix B. One drawback of existing approaches, however, is that although controller implementation is decentralized, their design is not necessarily so. This means that the controllers are designed jointly based on full knowledge of the entire power system model. As DER penetration is growing at an unpredictably high rate, grid operators must make sure that if more and more DERs are installed in the future the existing controllers would not need to be retuned or redesigned from scratch. In other words, the DER controllers to be designed must have plugandplay capability.
A control method called retrofit control recently proposed in [5, 18] can fulfill this objective. A brief summary of this approach is presented as follows. The dynamical system , where is input, is said to be stable if the autonomous system under is asymptotically stable. Consider a power system integrated with solar and wind farms. For , let the dynamic model of the DER connected to the th bus be rewritten as