An Analytical Model for Small Signal Stability Analysis in Unbalanced Electrical Power Systems

: This paper presents a general model to carry out a small signal stability analysis in electrical power systems (EPSs) that operate in unbalanced conditions. The classic traditional approach is based on a phase representation where it is supposed that the triphasic electric grid does not suffer any variations during its operation. With the presence of unbalances in transmission lines and loads, it is necessary to develop a general model that allows answering the needs and challenges with which modern electric systems must deal. The present work ﬁrstly intends to address the three-phase representation of an EPS, including the controllers. The proposed model is applied to a classical test system, a three machine-nine bus system, considering all the dynamic and algebraic variations associated with angular stability analysis. The proposed approach to small signal stability analysis shows results that differ from the classical analysis. The results are backed up with time domain simulations, and therefore, these results could be used in the calculation of the controllers that operate in unbalanced multimachine systems.


Introduction
The study of stability of electrical power systems (EPSs) is relevant for their operation and planning, as it thoroughly describes their behavior when facing events such as small and large perturbations [1]. The classical analysis assumes that the transmission grid is balanced and does not present large variations [2]. This allows making phase modeling suppositions and implementing simplifications and equivalent models.
However, with the dependency of the world on electricity and the wide variety of users and loads in the distribution networks, it is not possible to assume that the power transmission through the lines is symmetric and that the charges will remain balanced. Furthermore, with the integration of new generation technologies with the distribution grid, the studies of the implementation and operation need to be more precise. This is why the model of the consumption loads can no longer be represented with the classical balanced method. The challenge for the new studies in electrical power system stability is to describe and model the systems considering these new scenarios in order create better analysis.
Studies on this kind of analysis reported in the international literature consist of the proposition and practical verification of new models that accurately describe all the parameters of the system. Some works study transients [3][4][5][6] and therefore need to solve the non-linear equations for power flow. In particular, concerning small signal stability, there is an additional problem not always considered, and it is the balance point that defines the linear model that is not fixed and varies with time [7].
Accordingly, the solutions proposed in the literature describe several types of approaches that solve this inconvenience [8]. Reference [9] was based on simulations and the use of estimation methods complemented with simulations so as to obtain the eigenvalues of the system. Additionally, mathematical modeling tools such as dynamic phasors are used to directly describe the system.
Estimation methods just allow knowing the types of oscillation modes, oscillation frequencies, and the damping coefficient. However, having a mathematical model that analytically describes the behavior of the system enables carrying out more specific analyses such as the calculation of the participation factors.
In References [10,11], results were obtained by using the estimation model of a reduced system, which consists of a distributed generator under different load conditions. The results indicate that the more unbalanced the load of the system, the greater the possibility of losing the stability of the permanent regime is.
Furthermore, multiple articles address small signal stability applied to power electronic interfaces in renewable generation systems [12][13][14]. In general, a dynamic phasors model is proposed based on "dq" coordinates, which are eventually validated through simulations. However, the techniques proposed are not applied to large-scale systems.
Reference [2] proposed a model based on the phasors used in ZPNtransformation to solve the problem of the variant balance point. However, the model is only applied to a single machine infinite-bus system. Then, the case in [15] uses a derivation of this model applied to a small system with the inclusion of induction motors and photovoltaic modules. Reference [16] applied this concept to microgrid analysis.
Despite the favorable results presented, there is no work that compares the classic mathematical models with those proposed in previous reports under real modeling conditions. This kind of study can verify the applicability of the models and can be possibly used to improve the studies of controllers applied to unbalanced system. In this article, we present a novel solution to the problem of small signal stability analysis for the multimachine systems under unbalanced conditions and compare it to the classical analysis in the search of a validation.
From the proposal in [2], where the model was only applied to a one machine test case, the present work models multimachine EPS in unbalanced conditions, then generalizes the small signal stability analysis, and finally, compares the results obtained with a monophasic model used in the classical literature. The simulation tools are DigSilent and Matlab to obtain the unbalanced three-phase power flow, as well as the linear equivalent model.
To the best of our knowledge, this kind of study cannot be performed in any available commercial software for large-scale systems. To name a very popular software solution, DigSilent PowerFactory is based on simulation and is not prepared for an eigenvalue analysis with unbalanced conditions [17]. In this work, we propose and validate an analytic derivation for the generalization of a small signal stability study. Thus, the solution of algebraic differential equations is not required. The methodology we introduce in this article covers an area that has not been deeply developed and aims at analyzing unbalanced systems with a large number of machines and buses, considering an analytical model of the entire system. Furthermore, there is no analytically-oriented model for controller tuning in unbalanced systems, and the proposed approach to small signal stability analysis opens the way for the adjustment of the controller parameters under unbalanced conditions. This work is at the cutting-edge of research in small signal stability analysis.
The present report is structured as follows: Section 2 presents the description of the problem and how it is approached. Section 3 describes the modeling and linearization results. Section 4 explains the case study and results obtained. Finally, Section 5 presents the conclusions and proposed future work.

Description of the Stability Problem for Three-Phase Unbalanced Systems
Reference [18] showed how an EPS can be described through a set of differential and algebraic non-linear equations. This equation system contains dynamic and algebraic variables of the system: where: • x ∈ R n represents the state variables of the machines; • y ∈ R m represents the the algebraic variables of the system; • u ∈ R d represents the input variables, in this case the parameter associated with the control systems of the generators; • f (x, y, u) represents the differential equations of the generators; • g(x, y) represents the algebraic equations associated with the generators and traditional load flow equations; • n, m represent the number of state variables and algebraic equations.
The calculation needed to directly and simultaneously solve the set of algebraic and differential equations, defined in Equations (1) and (2), is numerically very complex because of the non-linearity associated with the equations of the system [18]. As an alternative path, the procedure used linearizes the equation system around an operation point defined through the power flow. Thus, the following linear equivalent is obtained: where: • ∆x represents linearized state variables; • ∆y represents linearized algebraic variables.
In classical small signal stability studies, it is considered that the input vector ∆u = 0, and therefore, the eigenvalues of matrix A sys are evaluated. In this context, the present article proposes a model that allows generalizing the structure of the A, B, C, and D matrices for the cases in which the system operates under unbalanced conditions. The approach used by the authors is based on what was proposed in [2]. A sequence component methodology is used, which consists of representing an unbalanced three-phase system, through three three-phase balanced systems denominated as: positive sequence, negative sequence, and zero sequence. The description of each system considers the following:

•
Positive sequence: this is the real system considering that it operates in balanced conditions. • Negative sequence: this has the same topology as the positive sequence, but the equivalent voltage sources of the generators are not represented. • Zero sequence: this has the same topology as the negative system, but only here, the grounding connection and the existence of neutrals are modeled.
The use of sequence components allows representing in the bases A, B, and C the variables associated with the algebraic equations. In general terms, matrix A does no suffer any modifications in the quantity of variables that define it. However, matrices B, C, and D are almost completely modified since each of the variables that define them are tripled in quantity.
Once the new A, B, C, and D matrices are obtained, the equivalent matrix A sys (p) is redefined, and a small signal stability analysis is carried out. This study intends to determine the eigenvalues of the system (3), which are obtained by solving the following equation: where: • I is the identity matrix; • λ represents the eigenvalue vector of the system (3).
From the classical approach of permanent regime stability, the system of Equation (3) will be stable, only if all the real parts of the eigenvalues obtained from Equation (4) have a negative real part [18]. A commonly used indicator for stability is the damping ratio, which indicates how much each eigenvalue contributes to the stability of the system. It can be calculated with Equation (5). where: • ζ the damping ratio; • σ the real part of the eigenvalue λ; • ω the real imaginary part of the eigenvalue λ.
When carrying out a study related to small signal stability in unbalanced systems, the main issue observed is the increased number of algebraic variables. As the classical analysis uses a phase modeling reduction, when the three phases of the system are analyzed, the number of algebraic variables studied is three times as much. This does not happen with the number of dynamic variables as they do not need new variables to describe the dynamics of the machine. This implies that the possible unbalance in the EPS is represented by the values of the matrices B, C, and D that include the algebraic variables.

Modeling the Dynamics of the Machines
In order to show the development proposed, the generalization of the oscillation equation of the electric machines is firstly contemplated. The modeling of the machines uses as a base Figure 1. The model used is called the "two-axis model", as it is a reduction of the complete model of the synchronous machine and does not include the subtransient terms and the flux linkage terms. For the detailed origin of this model, see [18]. The main modification that the oscillation equation suffers is the introduction of the associated losses to the negative sequence. The following is considered for the oscillation equation in a synchronous machine of the current algebraic variable associated with the negative sequence, given by the term 1 where: • ω i the angular speed of the machine (rad/s); − R s i ) the resistance of core losses due to unbalance (pu).

Internal Voltage And control Systems' Differential Equations
The synchronous machine model has its origin in the application of the coordinates in the direct and quadrature axes ("dq"), as was described in [18]. The internal voltage shown in Figure 1 is decomposed intoĖ d q i andĖ q q i . The differential equations that define them are presented below: where: • T do i and T qo i the direct transient time constants and of quadrature (s); • X d i and X q i the direct reactance and of quadrature (pu); • X d i and X q i the transient direct reactance and of quadrature (pu); • I d i and I q i the direct axis current and of quadrature (pu); • E f d i the equivalent field voltage (pu).
Synchronous machines can adjust their internal voltage through what is known as field control or excitation current/power. We will use the standard IEEE Type 1 model to represent the voltage control monitoring dynamics: where: • V re f i the reference voltage of the field control (pu); In the differential equation forV R i appears the concept V Ti , which corresponds to the voltage in the terminals/clampings of the machine and consists of an algebraic variable. Thus, in this case, we will consider that V Ti = [V a i + V b i + V c i ]/3 in the same way as in [2]. A method to reduce the oscillations that may occur in the EPS operations consists of introducing the power system stabilizer (PSS). This controller takes a speed signal and transforms it into a voltage signal, which powers/sustains/fuels the exciter system. This device is composed of three blocks: gain, washout, and phase compensators. The dynamic equations that characterize it are: where: • X 1i , X 2i , X 3i the internal variables of the PSS controller; • K PSS i the gain of the PSS controller;

Stator Algebraic Equations
In order to write the algebraic equations of the model, the equivalent circuits associated with the negative and zero positive sequence networks are analyzed. The following denomination is used to make reference to the voltage and current of the system. For any given machine i, it is considered that: the real part of the voltage β phase sequence ρ (pu); the imaginary part of the voltage β phase sequence ρ (pu); • I ρ β(Im) i the imaginary part of the current β phase sequence ρ (pu); • θ β i the current gap of a β phase (rad).

Positive Sequence
The positive sequence corresponds to the model displayed in Figure 1: considering: Then, the equation is: When separated into the real and imaginary parts, it is: Additionally, it is necessary to add the equations of the I d and I q currents:

Negative Sequence
From Figure 2, the following equivalent circuit is obtained, where R neg i and X neg i represent the negative sequence impedance parameters of the machine: When solving the Kirchoff law: By separating into the real and imaginary parts, the resulting equations of the machine are: where: • R

Zero Sequence
From Figure 3 and considering that R cer i and X cer i represent the zero sequence impedance parameters of the machine, the resulting equations to be used are the following: When posing Kirchoff's voltage law: When separating the real and imaginary parts, the resulting equations of the machine are: where: • R cer i the zero sequence resistance (pu); • X cer i the zero sequence reactance (pu); • R cer n i the neutral resistance (pu); • X cer n i the neutral reactance (pu); • V cer aRe i the real part of the zero sequence voltage (pu); • V cer aIm i the imaginary part of the zero sequence voltage (pu);

Network Equations
In order to obtain the network equations, which consist of those algebraic equations associated with the power flow, it is necessary to make a generalization of Ohm's law. For only one element, the relation between voltage/power and current is given by the following expression: where: •Ż the impedance element (pu); •Ẏ the admittance element (pu). AsẎ = (Ż) −1 is fulfilled for the multimachine case, an admittance element matrix is used as follows: When considering the electric system in balanced conditions, the admittance matrix has a range given by the number of bars that the network has.
In the case of an EPS under unbalanced conditions, it is not possible to use a phase equivalent model, and it is necessary to model the elements in its three components a, b, and c. In this case, each element Y ij of the admittance matrixẎ system corresponds to the other 3 × 3 matrix: where: • Y αβ ij the admittance between bus i and bus j between phase α and phase β (pu).
In this case, the three-phase representation of the lines and transformers connected to the different bars of the EPS is considered. Thus, the network equations are used for each phase β ∈ (a, b, c) and for each i bar.
Therefore, for each bar, there are three active and reactive power equations. Overall, there will be six algebraic equations for each bar. To represent the foregoing, the factor γ is used for all the necessary equations of the linear model.

•
Algebraic equations for a generation bus: • Algebraic equations for a load bus:

Linearization
The equivalent model for the small signal stability analysis is obtained by linearizing the differential and algebraic equations based on the operation points obtained from the power flow.
For an n bar and m machine system (m < n), the following system matrices are defined. The variables of the general system are displayed in Equation (39) where "x i " contains the state variables for an i bar, with a total of ten state variables. Simultaneously, "y i " contains fourteen algebraic variables. Below, the linearized expressions are shown: The resulting equation of the linearization can be founded in Appendix A.

Summary of the Matrices' Structures
To form the corresponding matrices of all the system, diagonal matrices are formed: For the stator equations: Because of space limitations, the following equalities were determined: The matrices of all the systems are formed in (43). The network equations for the generation bars are shown below: (1) . . . where: Following the same logic, the network equations for the charge bars correspond to: 0 =     D 6 (m+1,1) · · · D 6 (m+1,m) . . . . . . . . .
where the A sys of Equation (3) has the following form: considering:

Testing System and Modeling Methodology
A three generator-nine bar system is used, widely found in the international literature, as a testing system to apply the proposed model. The topology of the system is shown in Figure 4, and the information of the elements is included in Appendices A and B. For the permanent regime symmetric operation condition, Table 1 shows the considered load distribution. In order to implement the unbalanced model described in the sections above, two load conditions are presented: the nominal values of the loads and the case previous to instability in the classical model. The two scenarios analyzed are presented below: • Scenario 1: balanced nominal charges; • Scenario 2: charge increase by 45%.
The model proposed in this paper allows analyzing an unbalanced system in its topology, such as unbalanced lines or operation points where the loads do not have the consumed power distributed in a symmetric form in its three phases. As a first approach to the unbalanced study, unbalanced load is considered, modeled as constant power and distributed in a proportional form among its three phases The three-phase unbalanced model with initial balanced conditions is used as a control case because it can directly be compared to the classical analysis. This case is not real since distribution systems are naturally unbalanced. Two more asymmetrical distributions of loads were studied, named Forms A and B. These are used to analyze the effect of different kinds of unbalance on the system stability. The way in which the asymmetrical load is distributed in both forms is shown in Tables 2 and 3.  Table 4 displays all the cases considered in order to evaluate the model proposed in this paper. The methodology for simulation used in this work (see Figure 5) starts first (1) with gathering the information of the system, especially focused on the parameters used in the unbalanced representation. The second step (2) is to create the model of the system in the software DigSilent, including all the unbalanced parameters. The third step (3) is to modify the loads of the system according to the study cases. In the fourth step (4), the load flow is calculated for every single study case, and the results are obtained. The fifth step (5) is to use the system parameters, load flow results, and initial conditions as the inputs for a Matlab script that calculates the resulting matrices A, B, C and D for each study case. In Step (6), the matrices are used for the eigenvalue analysis of the system, and the damping ratios for these eigenvalues are also calculated. Finally, (7) the time domain simulation is carried out in DigSilent to contrast the results.
The main criteria for deciding about the stability of the eigenvalues is based on eigenvalue analysis. If an eigenvalue has a real part that is positive, it is considered unstable [18]. The system must not have any unstable eigenvalue to be classified as stable.
To check the validity of the results of the proposed model, a time-based simulation of the system is constructed. The simulation includes a small perturbation and contrasts the results.

Results
The resulting eigenvalues obtained for Scenario 1 (nominal loads) were 30; the system has 30 dynamic variables and 140 algebraic variables. In Figure 6 is shown the eigenvalues close to the origin considering Cases 1, 2, 3, and 4. The eigenvalues are determined from the Matlab script. The classic model is considered for Case 1, and Cases 2, 3, and 4 use the unbalanced model. It is possible to observe that the systems are stable since none of the eigenvalues has a positive real part. However, unbalanced models tend to be less stable because some real parts tend to move closer to the right of the complex plane.
Analyzing the eigenvalues of only the unbalanced Cases 3 and 4, they have a further reduction of the stability due to the considered load distributions. To verify this situation, Table 5 presents the results of the damping percentages obtained for all the cases of Scenario 1. With respect to the damping percentages, when comparing all the modeled cases, there is a tendency in Cases 2, 3, and 4 to be less stable/steady. This is verified as the damping values are lower in the balanced case (Case 1). To check these results, the time domain simulation is carried out. A three-phase short circuit in bus 8 is activated at Second 10 of the simulation for one millisecond of duration. Figures 7 and 8 show the resulting angle of the rotor and power delivery from Generator 2.  It is clearly visible that Cases 1 and 2 with the balanced initial condition are both stable. Cases 3 and 4 are stable, but the steady state performance has sustained variations.
The results obtained from Scenario 2 are presented in Figure 9, which shows the eigenvalues calculated with the Matlab software.
As was expected, the system tends to instability in all the cases due to the load increase (see Table 6). This indeed occurs with the classical model. However, with such representation, it is not possible to model an asymmetrical load distribution in the three phases. Scenario 2 is a load situation that stresses the operation of the system. Case 5 is stable, and the analysis of the eigenvalues shows that the system is on the verge of instability. The cases that used the proposed approach, Cases 6, 7, and 8, are unstable because all of them have real positive eigenvalues. This clearly means that the proposed model has opposite results to the classical analysis. The different load conditions of Cases 7 and 8 generate different results and show the unstable eigenvalue to be closer to the origin compared to Case 6.  Although Case 5 is stable despite the operation condition, it has very low damping percentages and is close to instability with the analytical model. The real time domain simulation of Scenario 2 showed that all of the systems were unstable in all conditions. This reinforces the idea that the classical model does not properly represents unbalanced conditions and that the proposed approach does a better performance.
To verify the instability in the time domain simulations, Figure 10 shows the frequency response for Cases 2 and 6. It is evident that Case 6 has an unstable condition as the frequency drops even without the short-circuit at Second 10. This behavior is identical in Cases 7 and 8.  Figure 10. Angles and frequency of Machine 2 in Case 2 versus Case 6.

Conclusions
The present research shows the necessary analytical methodology to carry out small signal stability studies in electric systems that operate in the unbalanced condition. In this context, a traditional testing system widely used in the international literature is used.
It is verified that the load variations and imbalances are increasingly more common in modern electric systems and cannot be represented by a classical model, which considers only one representation per phase because the analytical system does not show the imbalance present in the time domain simulations.
These results are particularly important in real systems where the imbalance can be present, such as distribution networks under different load conditions. For future work, it is intended to apply this model to real systems with a large number of variables and to consider the presence of self-sustained perturbations in time. The new approach presented in this paper to small signal stability analysis opens the way for the adjustment of the controller parameters in unbalanced multimachine systems.
Additionally, the research team is currently working on a generalized model, which allows including non-conventional renewable energies and aspects related to distributed generation so as to represent all the new elements influencing the stability of electrical power systems.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The linearization of the general equations of the system is given below.
Appendix A.1. Differential and Mechanic Equations where:

. Stator Algebraic Equations
Positive sequence: Negative sequence: Zero sequence:

. Network Equations
By linearizing the active and reactive expressions, we obtain the linear equations (A15) and (A16) for the generation bars.
By linearizing the active and reactive expressions, we obtain the linear equations (A17) and (A18) for the charging bars.

Appendix B
The details of the testing system are presented below.