Dynamic modeling and stability analysis of power networks using dq 0 transformations with a uniﬁed reference frame

. The dq 0 reference frame has become popular for modeling and control of traditional electric machines and small power sources. However, its widespread use for modeling and analysis of large-scale, general power systems is still a pending issue. One problem that arises when considering dq 0 models is that they are typically based on local reference frames, and therefore linking different models is not straightforward. In this paper we propose to approach this problem by modeling the network and its components using a dq 0 transformation that is based on a uniﬁed reference frame. We demonstrate this idea on the basis of synchronous machines and photovoltaic generators, and we also establish a dq 0-based dynamic model of a transmission network. The resulting models all use a uniﬁed reference frame, and therefore can be directly linked to each other in simulation and analytically. The paper is accompanied by a free software package (Levron, Y. and Belikov, J. Toolbox for Modeling and Analysis of Power Networks in the DQ0 Reference Frame. 2016. www.mathworks.com/matlabcentral/ﬁleexchange/58702) that constructs the proposed dynamic models and provides tools for dynamic simulations and stability studies based on dq 0 quantities.


INTRODUCTION
In recent years, the increasing penetration of small distributed generators and fast power electronics based devices has given rise to new challenges in modeling the dynamics of power systems.This becomes evident when considering large-scale systems for which it is especially important to maintain the balance between accuracy and complexity.This challenge has led to new dynamic models of power systems that are based on direct-quadraturezero (dq0) quantities [1].Such models are not as general as three-phase (abc)-based (EMTP-like) models and are advantageous mainly when the transmission network and units are balanced.However, dq0 models combine two properties of interest: similar to transient models, dq0 models are derived from physical representations, and are therefore accurate at high frequencies, and as a result the assumption of slowly varying signals is not re-quired.In addition, dq0 models are time-invariant, which allows defining an operating point, and enables smallsignal analysis.A detailed comparison of simulation techniques based on the abc and dq0 reference frames can be found in [2,3].
The dq0 transformation has been traditionally used in transient analysis of electric machines [4] and is increasingly used today for modeling distributed sources, complex loads, renewable generators, and on devices based on power electronics [5][6][7].However, widespread use of this transformation for modeling large-scale power systems is still a pending issue.Toward this end, one problem that arises when considering dq0 models is that units are typically based on local reference frames, and therefore linking different models is not straightforward.The idea of modeling a transmission network and loads in a common dq reference frame was presented in [1,8], assuming that the network consists only of resistors and inductors.These works consider a common reference frame, meaning that all the units and the network share the same reference, see [9].Work [10] presents dq0based models of three-phase networks with RL elements.Work [11] develops a small-signal dq0 model of a microgrid that includes synchronous machines and electronically interfaced distributed generators.In addition, an attempt to derive a model of a general network is presented in [12], where it is shown how to construct a state-space model of the network that is nonminimal.A further extension of this paper is presented in [13].
In this paper we continue these ideas and propose a method for modeling the network and its components using a dq0 transformation that is based on a unified reference frame, where the main objective is to create dq0based models of large-scale power systems that are both accurate and time-invariant.The resulting models all use a unified reference frame, and therefore can be directly linked to each other in simulation and analytically.We open the paper by recalling the basic transformation from one reference frame to another and demonstrate this process using the classic example of a synchronous machine connected to an infinite bus.The paper then extends this example and presents a dq0-based dynamic model of a general transmission network, using ideas presented in [10,14] and based on the network topology defined in MATPOWER [15].We also develop models of several typical units.This approach is demonstrated numerically in a dedicated software tool available in [16].Several numeric examples are described in detail: a 4-bus system showing comparison between quasi-static, abc-, and dq0-models; a 14-bus system illustrating the joint dynamics of synchronous machines, photovoltaic generators, and the transmission network; and a large 200-bus system demonstrating small-signal analysis.

PRELIMINARIES AND A MOTIVATING EXAMPLE
Consider a reference frame rotating with an angle of θ (t).For instance, in a synchronous machine, θ (t) is typically selected to be the rotor electrical angle.Let ζ represent the quantity to be transformed (current, voltage, or flux), and use the compact notation with Further in this paper we aim to explore an alternative reference frame that rotates with the unified angle ω s t, where ω s is the steady-state system frequency or the frequency of an infinite bus, see [17] for more details.We will show that such selection provides several advantages in the connectivity of a large variety of system components and also can be used to derive compact and easyto-use small-signal models describing the power system dynamics.The dq0 transformation can be redefined with respect to a new reference frame rotating with an angle ω s t by direct substitution of θ = ω s t in (1) as A formula that allows conversion of signals from the standard reference frame (defined with respect to θ ) to the new frame (defined with respect to ω s t) is given as where δ (t) = θ (t) − ω s t + π/2, the variables ζd , ζq are defined with respect to θ , and ζ d , ζ q are defined with respect to ω s t.
To explain the idea of a unified reference frame, we start by recalling a classical model of a synchronous machine connected to an infinite bus, and present both models using a unified reference frame.Assume a simplified synchronous machine represented as an ideal voltage source behind a synchronous inductance Ld .The machine is connected to an infinite bus, as shown in Fig. 1.Both the infinite bus and the internal synchronous machine are represented as voltage sources.The synchronous machine voltage is ṽd = 0, ṽq = V e , v 0 = 0, with a reference angle of θ .In this example θ is the electric angle of the rotor in respect to the stator.The infinite bus has a constant frequency of ω s , so its voltage is given by v d = V g , v q = 0, v 0 = 0, with a reference angle of ω s t.Now the goal is to construct a dynamic model of the complete system based on dq0 signals.However, a potential problem is that the two voltage sources are defined with respect to two different reference frames (θ and ω s t).To solve this, we choose ω s t as a unified reference frame for both the infinite bus and synchronous machine, and construct a model of the synchronous machine using signals in this reference frame.The synchronous machine voltage is now obtained by substituting ṽd = 0, ṽq = V e , v 0 = 0 in (3), which yields In addition, the dynamic behavior of the angle δ is described by which is the classic swing equation.The term J is the rotor constant of inertia, poles is the number of machine poles (must be even), P m is the mechanical power, and K d is the damping constant.The three-phase power can be computed as P 3φ = 3 2 (v d i d + v q i q + 2v 0 i 0 ).Let δ = φ 1 , then the combination of ( 4) and ( 5) results in a statespace model given as v d =V e cos(φ 1 ), v q = V e sin(φ 1 ), v 0 = 0.
Note that the model represents only the machine's voltage source and does not include the synchronous inductance, which is modeled separately in the next section.

MODELING POWER SYSTEM COMPONENTS USING A UNIFIED REFERENCE FRAME
In this section we show that this idea can be extended and that complex networks and various power system units can be connected to each other by using signals that are defined with respect to the unified reference frame.This allows us to obtain a complete model of the system that is based entirely on dq0 quantities, and thus provides a constructive method to analyze the system's dynamic behavior.

State-space models of linear passive transmission networks
We open this section by developing models of single passive components [18], using the unified reference frame.The discussion that follows shows how to combine these models to describe the dynamics of general transmission networks.The dynamic model of an inductor is given in the abc reference frame as and can be converted to the dq0 frame as follows.Observe that the differentiation of T ω s results in and the derivative of T ω s can be expressed as Substitute ( 7) and ( 9) into ( 8), and use relations (2) to get This equation describes a state-space model of the threephase inductor.Analogously, the model of a capacitor C is given as And for a resistor R the model is given by simple static relations where I 3 denotes the 3 × 3 identity matrix.Similarly, by combining elementary passive components, any balanced transmission network can be modeled based on the unified reference frame.Define the dq0 signals v d,n , v q,n , v 0,n to be the voltages on bus n; i d,n , i q,n , i 0,n to be the injected currents to bus n; and where N is the number of buses.Assume a network with the standard branch [15] as in Fig. 2.
This network can be represented in the minimal state-space form using dq0 signals as [13] where In addition, often disconnected buses should be eliminated.The need for this arises in several situations: • Sometimes certain buses are not connected to either a generator or a load.In this case the current injected into the bus is zero.• Frequently loads are modeled as shunt elements and are integrated into the network model.In this case load buses appear as disconnected buses with zero current.• In many scenarios there is a need to analyze the dynamics or stability of a certain subset of units in the network, typically only the generators.In such cases elimination of the disconnected buses results in a simpler dynamic model in which the inputs and outputs relate only to the required subset of buses.This is usually done after integrating the loads into the network model as shunt elements.
The dynamic model after disconnected buses have been eliminated is where ũ and ỹ represent the input/output following the reduction procedure.Bus elimination is achieved by controlling the inputs of the buses being eliminated, so that the corresponding outputs are zeroed.Assume for instance that the ith input and ith output are eliminated.This is done using the dynamic model output equation y = C ξ ξ + D ξ u.In case the matrix D ξ of the original dq0 model is diagonal, the output y i may be zeroed by controlling the ith input so that u i = −D −1 i,i C i ξ , where C i is the ith row in matrix C ξ .Thus, the input u i and the output y i are eliminated and will not appear in the reduced model.In more complex systems, there is a need to transform the state vector and to compute a new dynamic model such that the corresponding rows in C ξ are also zero.This may be achieved by a LU decomposition of the relevant rows in C ξ .More details regarding this procedure may be found in [16].

Physical synchronous machine
For completeness, consider a more sophisticated (physical) model of a synchronous machine [4].The model presented herein captures the interaction of the direct-axis magnetic field with the quadrature-axis mmf and of the quadrature-axis magnetic field with the direct-axis mmf, as well as the effects of resistances, transformer voltages, field winding dynamics, and salient poles.The parameters are explained in Table 1.
By following [4] and omitting laborious algebraic manipulations, the resulting state-space model of a synchronous machine in the dq0 reference frame (with respect to ω s t) is given by with outputs defined as where . In this model, the state variables are φ 1 = λ d , φ 2 = λ q , φ 3 = λ 0 , φ 4 = λ f , φ 5 = ω, δ = φ 6 ; the inputs are v d , v q , v 0 , v f , T m ; and the outputs are i d , i q , i 0 , i f , ω.Unlike the simplified model (6), this model includes the inductance terms L d , L q , L 0 .
A convenient property of the model presented above is that its inputs and outputs are defined with respect to the unified reference frame rotating with ω s t, and therefore it can be directly connected to the network.For instance, connecting the synchronous machine model to Mechanical torque an infinite bus is immediate.An infinite bus with a frequency of ω s is defined in this reference frame by v d = √ 2V g , v q = v 0 = 0, where V g is the infinite bus RMS voltage.Therefore, a synchronous machine connected to an infinite bus can be modeled by direct substitution of these voltages into (3.2).

Photovoltaic generator
Similarly to the synchronous machine, photovoltaic inverters can also be represented in the unified reference frame that rotates with ω s t.Several models of photovoltaic generators are available in the recent literature, as reviewed in [10].Regardless of the model used, dq0 signals defined in respect to a local reference frame can be converted to a unified reference frame.Here we demonstrate this process using a photovoltaic inverter model based on [19].The inverter is modeled by where P pv denotes the DC power at the output of the photovoltaic array, C b is the bus capacitor (connected at the input of the switching devices), v b is the bus capacitor voltage, and v re f is the reference voltage for the bus capacitor control circuitry.The term 3 2 (v d îd + v q îq ) is the inverter output power in terms of dq0 signals.Based on [19], a simple PI controller is used here to control the bus voltage v b .Hence, K p and K i denote the coefficients for the proportional and integral terms, respectively.The second part describes the inverter output capacitor as and where C sh is the inverter output capacitance.

NUMERIC TEST CASES
This section shows several test cases that illustrate the main idea of modeling and analysis in the unified reference frame.Three networks with 4, 14, and 200 buses, whose parameters are taken from [15], are considered.
The first example presents a brief comparison between the quasi-static, abc, and dq0 models.The second and third examples are devoted to the small-signal analysis of the 14-and 200-bus networks.

Comparison of quasi-static, abc, and dq0 models
Figure 3 illustrates comparison between conventional abc and dq0 models presented in terms of sparsity and total number of nonzero elements for various test-case networks taken from the MATPOWER database [15].The sparsity indices are computed as the ratio between zero elements in all system matrices and the total number of elements.Observe that both models have approximately the same number of nonzero elements; however, the abc model is slightly sparser.Consider now a simple 4-bus network.Its single-line diagram is shown in Fig. 4.
The quasi-static and dq0 models were constructed using the software package available in [16].The abc model was implemented using MATLAB toolbox Simscape Power Systems.The transient behavior is analyzed under changing operating conditions as depicted in Fig. 5.The input voltage is subsequently changed from the initial value as follows.The ramp signal with the slope of 200 and upper limit of 10 kV is used to change the d component of the input voltage v d,1 at time t = 0.02 s.Then the input is stepped from 10 to 20 kV at t = 0.12 s.The output current dq components are measured on bus 2. Observe that all models coincide in the steady state as expected.However, the abc and dq0 models are able to more accurately describe the transient behavior as they capture high-frequency phenomena.

A 14-bus network
Consider now the IEEE 14-bus system.This network was modified by replacing synchronous machines connected to buses 6 and 8 by photovoltaic generators.Simulation parameters of synchronous machines are given in Table 2.Note that mechanical input powers (P m ≈ T m ω s ) are used as new external inputs.Parameters of the infinite bus are v d,1 = 1.5 • 10 5 V and v q,1 = v 0,1 = 0 V.Moreover, ω s = 2π50 rad/s is used to define the frequency of the unified reference frame.Simulation parameters for photovoltaic generators are given in Table 3.
Further we perform a small-signal stability analysis of an interconnected system, which is done in several steps: first the system's operating point is computed by solving the system power flow equations and converting the resulting complex voltages and currents to dq0 quantities.Next, unit models are linearized in the neighborhood of this operating point, and the overall system is described using state equations.We start with the nominal case and use the initial values provided in Tables 2  and 3.An array of Bode plots representing the frequency responses of the linear model is illustrated in Fig. 8 for P m,2 → P 2 and P m,3 → P 2 input-output pairs.These figures demonstrate that synchronous machines are weakly coupled at low frequencies but nonetheless affect each other at a certain resonance frequency of about 30 rad/s.In addition, as expected, high frequencies in the mechanical powers are filtered by the network inductances and synchronous machines inertia.The transient process is shown in Fig. 9 for active powers.As expected, the infinite bus reacts to the changes in power generation.Table 4 provides numeric information for the largest eigenvalues.These results are then visualized in Fig. 10(a), which shows the root locus of dominant poles for 3 intervals of changing operating conditions.In addition, steps are depicted in Fig. 10(b).
Next, the dynamic behavior is analyzed under the changing operating conditions as presented in Fig. 10(b).The nonlinear models of units are then re-linearized in the neighborhood of new operating points when the transient process is complete and the system is in the steady state.Time instances are selected as t ∈ {1.9, 3.9, 5.9, 8}.Signals are subsequently changed from their initial values as follows.The mechanical powers P m,2 and P m,3 of SG1 and SG2 at buses 2 and 3 are increased from 120 to 180 MW at t = 2 s and from 60 to 90 MW at t = 6 s, respectively.The DC power P pv,6 at the output of the photovoltaic array PV1 connected to bus 6 is decreased from 30 to 15 MW at t = 4 s, while P pv,8 is kept constant (30 MW).From Fig. 9 it can be seen that the provided changes in input powers directly affect the steady-state conditions and cause a change of the operating point.The presented analysis may assist in understanding the dynamic behavior of the system.For example, from Fig. 10(a) it follows that poles move toward the origin, which means that the response of the system becomes slower with the increase of external powers.In addition, observe that a change in the DC power P pv,6 of the photovoltaic array connected to bus 6 has a minor effect (zoomed area in Fig. 10(a)) on the overall dynamics of the system, because poles are slightly moved from their positions.This is due to the fact that the photovoltaic inverter connected to the grid is less powerful than the respective synchronous machines.Clearly, if the amount or the overall capacity of photovoltaic energy sources is increased, the effect will be more significant.This, in particular, means that the proposed approach may help to estimate the amount of renewable energy produced by PVs within the stability limits.

A 200-bus network
In this section, we show that the proposed approach can be applied for the analysis of larger systems.In particular, consider a 200-bus system.A detailed description and numeric values can be found in [15] based on [20].Models of the transmission network, generators, and loads are constructed using [21].Generators are described using the swing equation ( 6), and loads are represented based on the balanced series RL impedances.The small-signal model is obtained similarly to the case of the 14-bus network.The complete (with attached units and loads) dq0 model has 1470 states and 48 external input signals (mechanical powers).All the components are modeled in the unified reference frame rotating with ω s = 2π50 rad/s.The sparsity pattern of the matrix A (sparsity index is 0.49%) and the distribution of the dominant eigenvalues are shown in Fig. 11.The dynamic behavior is further analyzed under changing operating conditions to demonstrate similarities and differences between quasi-static and dq0 models.Scenario (similarities): In the first test case the system is simulated under the nominal conditions, and then the overall active power demand of all loads is increased by 15% and 35% from 2229 to 2563 and 3009 MW, respectively.Table 5 provides numeric information, in which the first three eigenvalues with the largest real parts are shown for both quasi-static (qs) and dq0 models.Observe that both quasi-static and dq models provide almost identical dominant poles.Scenario (differences): Consider first the case when the moment of inertia (J) of the synchronous machine connected to bus 196 is changed.The moment of inertia is first increased and then decreased by 4 times from the nominal value.The dominant poles of the quasi-static and dq models are similar for higher values but differ for a smaller moment of inertia.Table 6 provides numeric information, in which the first several eigenvalues with the largest real part are shown.Observe that the dq model predicts that the system is unstable, but the quasistatic model fails to predict this instability.
Consider now the case when several line resistances are changed.In particular, we change the branch resistances for the lines connecting buses 135 → 133, 136 → 133, and 189 → 187.First, branch resistances are increased by 1.5 times and then decreased by 5 times from their nominal values.Table 7 provides the numeric information, in which the first three eigenvalues with the This example shows that the traditional small-signal stability analysis based on quasi-static approximation of the network has to be used carefully as in some cases it may give misleading results.

DISCUSSION AND CONCLUSIONS
A developing solution to model the dynamics of largescale power systems is to use dq0 quantities, assuming that the transmission network is symmetrically configured.This approach combines two properties of interest: similar to transient models, dq0-based models are derived from physical models and are therefore accurate at high frequencies.In addition, similarly to timevarying phasor models, dq0 models are time invariant.This property allows one to define an operating point and enables small-signal analysis.As a result, dq0-based analysis is expected to be beneficial when considering the stability of large-scale power networks that include a variety of distributed and renewable power sources.
A current challenge is to merge various dq0-based models appearing in the literature to obtain a complete model of a large power system.In this paper we approach this problem by representing various dq0 models in a reference frame rotating with a unified angle.This enables direct connections between units of different types and the network and provides a means to perform a smallsignal stability analysis of large-scale systems.This approach is demonstrated on the basis of two standard units (synchronous machine and photovoltaic generator), and in addition, a model of the transmission network based on the unified reference frame is provided.The paper is accompanied by a free software package available in [16] that constructs the proposed dynamic models and provides tools for dynamic simulations and stability studies based on dq0 quantities.Three particular examples are presented in this paper.

Fig. 3 . 4 Fig. 4 .
Fig. 3. Comparison between abc and dq0 models.The top plot presents the sparsity index, and the bottom plot illustrates the total number of nonzero elements.

Fig- ure 6
shows a single-line diagram of the system, which is composed of two 300 and 150 MW synchronous generators and two 20 MW photovoltaic energy sources.In this figure, the values of P in MW and Q in MVar represent the operating point, which is obtained by solving the power flow equations.Bus 1 represents an infinite bus modeled by a voltage source.The complete system is constructed in the unified dq0 reference frame.A sketch of the signal flow diagram of the complete system is shown in Fig.7.

2 SG1P 2 Fig. 6 .
Fig.6.Single-line diagram of a 14-bus system with two synchronous machines (SG) and two renewable energy sources (PV).Values of P and Q denote the operating point (power flow solution).

6 Fig. 10 .
Fig. 10.Eigenanalysis: root locus of the dominant eigenvalues when external input signals are changed.Dotted arrows indicate the direction of poles change.

Fig. 11 .
Fig. 11.Sparsity pattern of the matrix A and the distribution of the dominant eigenvalues for the dq0 model of the feedbackconnected 200-bus system.

Table 4 .
Change of the dominant eigenvalues

Table 5 .
Dominant poles: increase in active power consumption

Table 6 .
Dominant poles: change in the moment of inertia

Table 7 .
Dominant poles: change in branch resistances largest real part are shown.Observe that although the quasi-static model provides quite similar results for the 'increase' case, it fails to predict instability when resistances are decreased.