Study on the Method for Analyzing Electric Network Resonance Stability

: With the increasing utilization of power electronic equipment in power systems, there has been an increase in the occurrence of oscillatory behavior from unknown sources in recent years. This paper puts forward the concept of electric network resonance stability (ENRS) analysis and tries to classify the above-mentioned oscillations into the category of ENRS. With this method, many complex power system oscillations can be analyzed with the linear network theory, which is mathematically mature. The objective of this paper is to establish a systematic approach to analyze ENRS. By introducing the s-domain nodal admittance matrix (NAM) of the electric network, this paper transforms the judgment of ENRS into the zero-point solution of the determinant of the s-domain NAM. First, the zero-points of the determinant of the s-domain NAM are proved to correspond to the eigenvalues of the system. Then, a systematic approach is proposed to analyze ENRS, which includes the identiﬁcation of the dominant resonance region and the determination of the key components related to resonance modes. The effectiveness of the proposed approach for analyzing ENRS is illustrated through case studies.


Introduction
Oscillations are common phenomena in power system operations [1]. Generally, power system oscillations can be divided into three categories: (1) oscillations of the generator shaft (torsional interaction); (2) oscillations among the generator rotors; and (3) inherent resonance in the electric network.
With the increasing utilization of power electronic equipment in power systems, the third category-namely the inherent resonance in the electric network-becomes more severe and complex. For some electronic equipment, a negative resistance effect appears within a particular frequency range, which may lead to resonance instability. In the Hebei Province of China, since 2011 the interaction of doubly-fed induction generator (DFIG) wind turbines with nearby series capacitors has often caused subsynchronous oscillations with a typical frequency of 3-10 Hz. This leads to abnormal vibrations of transformers and wind generator tripping [2,3]. Moreover, unattenuated subsynchronous oscillations once occurred in the dc side of the Xiamen VSC-HVDC project in the Fujian Province, with typical frequency of about 25 Hz [4]. Furthermore, torsional oscillations once appeared in the Xinjiang Region, which is regarded as an interaction between generator shafts and the local power grid upon the integration of permanent magnetic synchronous generator (PMSG) wind farms [5]. Oscillations in the scenario of wind farm integration through a modular multilevel converter-based HVDC (MMC-HVDC) have also been reported [6,7]. Therefore, the inherent resonance structure and its relevant stability in power systems should be fully researched.
The inherent resonance structure of an electric network contains information on four aspects of the network: (1) the number of resonance modes within a target frequency range, and for each of the

Definition of System Eigenvalue
The dynamic behavior of an electric network is governed by the Kirchhoff current law (KCL), the Kirchhoff voltage law (KVL), and the dynamical equations of any energy-storage elements (e.g., capacitors, inductors). The former two are algebraic equations, while the latter one is a first-order differential equation. When an electric network consists of components of lumped and time-invariant parameters, a single-input single-output (SISO) state-space model can be established [15]. Based on the above-mentioned three equations, where {A, T} together make up the system matrix of the state space; x is the vector of state variables, including the state of energy-storage components (for example, capacitor, voltage, and inductor currents); y i is the voltage of the i-th node; u j is the injected current at the j-th node; and both column vector b and row vector c are one-dimensional constant vectors. Applying the Laplace Transform to Equation (1), one can obtain sTx(s) = Ax(s) + bu j (s) y i (s) = cx(s).
Then the transfer impedance in the s-domain is where det(sT − A) and (sT − A) * represent the determinant and adjoint matrix of sT − A, respectively. According to the theory of linear algebra, det(sT − A) is the system's characteristic polynomial, and its zero-points are the system's eigenvalues. The system is stable if all the eigenvalues are located on the left half-plane of the complex plane. Note that as the network components are described by lumped and time-invariant parameters, the state-space model has a limited number of resonance modes.

Eigenvalue Acquisition Based on the s-Domain Nodal Admittance Matrix
Beside the state-space model, the network eigenvalues can also be acquired by the s-domain nodal admittance matrix. In the s-domain, the nodal voltage equation of the electric network is where Y(s) is the s-domain nodal admittance matrix, V node (s) and I node (s) are the vectors of node voltages and node injection currents, respectively. Considering the same SISO system as in the previous subsection, suppose that the j-th node injection current is i j (s) and the i-th node voltage is v i (s). Hence, Equation (4) can be rewritten as where the column vector b j has zero elements except for its j-th element which equals 1 (one); and similarly, row vector c i contains 1 only in the i-th element. The transfer impedance in the s-domain can be calculated as where det[Y(s)] and Y(s) * represent the determinant and adjoint matrix of Y(s), respectively. Note that not all the elements of Y(s) are polynomials of s, for example, the operational admittance of an inductor has s in the denominator (1/sL). However, there will exist a proper integer m to make s m det[Y(s)] a polynomial of s. Further compared with Equation (3), it can be seen that s m det[Y(s)] is also the system's eigenpolynomial. Therefore, det (sT − A) = 0 and det [Y(s)] = 0 will have identical roots (except for s = 0). The nonzero roots of det [Y(s)] = 0 are actually the system's eigenvalues, which means that the state-space model and the s-domain nodal admittance matrix are consistent and both solve for the system's eigenvalues.
It should also be noted that conditions for network components with lumped and time-invariant parameters are not required when establishing Y(s). Hence, compared to the state-space methods, the Y(s) method is more applicable in analyzing ENRS, especially for the general electric network with components of distributed or frequency-dependent parameters. It can be proven that an electric network with components of distributed parameters will have infinite resonance modes.

Case Verification
In this subsection, the consistency of using the state-space model and Y(s) in solving for a system's eigenvalues will be verified through a simple case.
For the RLC series circuit shown in Figure 1 (R < 2 √ L/C), the state-space model can be expressed as d dt The system's eigenpolynomial is hence In the s-domain, the nodal admittance matrix of the circuit is which has the corresponding determinant It is clear that Equations (8) and (10) have identical zero-points if s = 0, which suggests that an identical result will be acquired when analyzing ENRS by both methods.

Case Verification
In this subsection, the consistency of using the state-space model and ( ) in solving for a system's eigenvalues will be verified through a simple case.
For the RLC series circuit shown in Figure 1 ( < 2 / ), the state-space model can be expressed as The system's eigenpolynomial is hence In the s-domain, the nodal admittance matrix of the circuit is which has the corresponding determinant It is clear that Equations (8) and (10) have identical zero-points if ≠ 0, which suggests that an identical result will be acquired when analyzing ENRS by both methods.

Eigenvalue Solution Based on Y(s)
Based on the analysis above, the solution for the eigenvalues of the electric network can be transformed into the solution of the following equation: The solutions of Equation (11) can be found by the Newton-Raphson method, which is mathematically mature and its rapid convergence (in five or six iterations) has been widely recognized. The Newton-Raphson iteration expression in this case is given by [12]: The meaning of vector and are the same as Equation (5), and the indexes i and j can be chosen to be equal for efficiency. Note that the CPU time for building ( ) and ( )/ are practically the same, as their building rules are similar.

Eigenvalue Solution Based on Y(s)
Based on the analysis above, the solution for the eigenvalues of the electric network can be transformed into the solution of the following equation: The solutions of Equation (11) can be found by the Newton-Raphson method, which is mathematically mature and its rapid convergence (in five or six iterations) has been widely recognized. The Newton-Raphson iteration expression in this case is given by [12]: The meaning of vector c i and b j are the same as Equation (5), and the indexes i and j can be chosen to be equal for efficiency. Note that the CPU time for building Y(s) and dY(s)/ds are practically the same, as their building rules are similar.

Systematic Approach to Analyze ENRS
In this section, a systematic approach to analyze ENRS will be proposed, which includes the identification of the dominant resonance region and the determination of the key components related to the concerned resonance modes.

Identification of the Dominant Resonance Region
For a particular resonance mode, the determination of the dominant resonance region is based on two indexes, the nodal voltage mode shape and the participation factor matrix.
Suppose the eigenvalue of the k-th resonance mode of the electric network is s k , then det [Y(s k )] = 0. According to the matrix theory, the determinant of a matrix equals the product of all its eigenvalues, hence Y(s k ) must have one eigenvalue that equals 0 (zero). Further suppose that Y(s k ) is an n-order diagonalizable square matrix, with eigenvalues λ 1 , λ 2 , . . . , λ n (for simplicity, let λ 1 = 0). The corresponding right eigenvectors are R 1 , R 2 , . . . , R n , which satisfy the uniqueness condition Based on the assumptions above, using the technique introduced in where Λ = diag(λ 1 , λ 2 , . . . , λ n ) is the diagonal matrix of the system's eigenvalues; Substituting Equation (13) into (4) and letting s = s k , one obtains Defining the mode voltage and mode injection current as respectively, then Equation (14) can be rewritten as Namely  Given . Hence Equation (15) can be approximately expressed as Hence, R 1 can be defined as the nodal voltage mode shape of the k-th resonance mode, which represents the relative magnitude and phase of each node voltage under the k-th resonance mode.
The nodal voltage mode shape can be used to locate the dominant resonance region and resonance type. In R 1 , some of the elements may have significantly larger magnitude than others, which means that the corresponding nodes are much more affected and can hence be considered as the dominant resonance region of the k-th resonance mode. Additionally, elements of R 1 may have close phases (see Section 4.1 as an example) or inverse phases (see Section 4.2 as an example). The presence of close phases usually corresponds to parallel resonance, for example, resonance caused by a reactive load and a nearby distributed capacitance. The opposite phases, on the contrary, usually mean series resonance, for example, resonance caused by a series capacitor and nearby transmission lines.
Equation (19) can be further written as where, P k = R 1 R T 1 is defined as the participation factor matrix of the k-th resonance mode; its element p k(ij) represents the contribution of j-th node injection current to the i-th node voltage when the k-th resonance mode is excited.

Determination of the Key Components
Among all the resonance modes in an electric network, the divergent modes are of the most concern, hence the key components that greatly affect these divergent modes should be determined. The determination of the key components is based on mode sensitivity analysis with respect to the parameters of these components. Theoretically speaking, mode sensitivity analysis can be directly executed if the target resonance mode is identified. However, in a large-scale electric network, it is quite difficult to analyze the mode sensitivity with respect to each of the system parameters. Hence, the identification of the dominant resonance region in the previous subsection, which contributes to reducing the range of mode analysis, is necessary. There is reason to believe that components in the dominant resonance region have a relatively large influence on the target resonance mode.
Suppose that a component parameter p is assumed to vary, the s-domain nodal admittance matrix can be further considered as a function of p, namely Y(s, p). In reference [13], the sensitivity of the k-th mode with respect to p is given as where s k is the eigenvalue of the target resonance mode and p 0 is the initial parameter of the concerned component. The vectors c i and b j are defined in the same manner as in Equation (6). Theoretically, c i and b j can be chosen arbitrarily and they will not influence the calculation result. However, in an actual calculation, if the i and j do not correspond to the nodes in the dominant resonance region, the accuracy of the final result may be deteriorated. Hence, for better accuracy of numerical calculations, in the proposed method the i and j will be chosen to be the same as the most dominant node: the node having the largest magnitude in the nodal voltage mode shape (suppose it is the m-th node). Another advantage is that c m Y(s, p) −1 and Y(s, p) −1 b m are transposes of each other, which helps accelerate the calculation in reality. With this in mind, Equation (21) can be rewritten as The above sensitivity index is in actual units. Under certain situations, it is more desirable to have a normalized sensitivity index with relative units instead of actual ones, so that indexes calculated from different components are comparable. Normalized sensitivities can be expressed as The acquired normalized mode sensitivities can be ranked to determine the key components related to the target resonance mode. For a certain mode sensitivity, its real and imaginary parts can also be separately considered. A relatively larger real part means the component mainly affects the damping factor of the resonance mode, while a relatively larger imaginary part means that the resonance frequency is more affected.
After determining the key components related to concerned resonance modes, if needed, the relevant parameters can be adjusted, so that the system will have better resonance structures. This parameter adjustment method will be demonstrated in future work.

Case Studies
In this section, three cases will be analyzed to demonstrate the effectiveness of the proposed method in analyzing ENRS. The three cases are the original IEEE 39-bus system, the modified IEEE 39-bus system with an integrated wind farm, and a wind farm integration system with an MMC.

Original IEEE 39-Bus System
The diagram of the IEEE 39-bus system is shown in Figure 2. Other system data can be found in reference [17]. The results of ENRS analysis on this system with a frequency range of 0~1500 Hz is listed in Table 1, marked as original system. Note that all the transmission lines are treated with a distributed parameter model during processing. Table 1 shows that the system has nineteen resonance modes within the target frequency range, and the damping factor of each resonance mode is sufficiently positive, hence the system is resonantly stable. damping factor of the resonance mode, while a relatively larger imaginary part means that the resonance frequency is more affected. After determining the key components related to concerned resonance modes, if needed, the relevant parameters can be adjusted, so that the system will have better resonance structures. This parameter adjustment method will be demonstrated in future work.

Case Studies
In this section, three cases will be analyzed to demonstrate the effectiveness of the proposed method in analyzing ENRS. The three cases are the original IEEE 39-bus system, the modified IEEE 39-bus system with an integrated wind farm, and a wind farm integration system with an MMC.

Original IEEE 39-Bus System
The diagram of the IEEE 39-bus system is shown in Figure 2. Other system data can be found in reference [17]. The results of ENRS analysis on this system with a frequency range of 0~1500 Hz is listed in Table 1, marked as original system. Note that all the transmission lines are treated with a distributed parameter model during processing. Table 1 shows that the system has nineteen resonance modes within the target frequency range, and the damping factor of each resonance mode is sufficiently positive, hence the system is resonantly stable.   Take the 131.6 Hz resonance mode as an example. Figure 3 shows the nodal voltage mode shape, and considering the layout restriction, the participation factor matrix is not given. The magnitude analysis of the nodal voltage mode shape helps to locate the dominant resonance region of the 131.6 Hz mode, which has been marked by a dashed grey rectangle in Figure 2. The phase analysis reveals that the voltage phases of all the buses are nearly the same, indicating that it is a parallel resonance. Further data inspection demonstrates that the resonance mode is caused by reactive loads and the nearby distributed capacitance in the dominant resonance region.  Take the 131.6 Hz resonance mode as an example. Figure 3 shows the nodal voltage mode shape, and considering the layout restriction, the participation factor matrix is not given. The magnitude analysis of the nodal voltage mode shape helps to locate the dominant resonance region of the 131.6 Hz mode, which has been marked by a dashed grey rectangle in Figure 2. The phase analysis reveals that the voltage phases of all the buses are nearly the same, indicating that it is a parallel resonance. Further data inspection demonstrates that the resonance mode is caused by reactive loads and the nearby distributed capacitance in the dominant resonance region.

Modified IEEE 39-Bus System
In this subsection, a DFIG-based wind farm is connected to the original IEEE 39-bus system through a long-distance transmission line with 50% series compensation (Figure 4). This modified system aims at conceptually verifying the first case mentioned in Section 1, i.e., that subsynchronous resonance (SSR) can be triggered by a DFIG-based wind farm with a series-compensation line.

Modified IEEE 39-Bus System
In this subsection, a DFIG-based wind farm is connected to the original IEEE 39-bus system through a long-distance transmission line with 50% series compensation (Figure 4). This modified system aims at conceptually verifying the first case mentioned in Section 1, i.e., that subsynchronous resonance (SSR) can be triggered by a DFIG-based wind farm with a series-compensation line.  The wind farm is aggregated by DFIG-based wind generators with a total rated-capacity of 1500 MW. The detailed impedance model of a DFIG wind generator is given in reference [18]. The impedance-frequency characteristic of the aggregated wind farm for the frequency range 0-100 Hz is shown in Figure 5. It indicates that a negative resistance effect (∠ > 90°) appears in the frequency range of 25-45 Hz, which may lead to SSR if the series resonance frequency also lies in this frequency range. Strictly speaking, the impedance model of a DFIG wind generator cannot be represented in the form of Y(s), because the outer loop controls (OLCs) of the d and q axes in the voltage source converter are not symmetric. However, the response time of the OLCs is usually larger than 200 ms, which means that the OLCs can be ignored for the considered frequency range studied. The impedance model in reference [18] is deduced based on the simplification above. The results of ENRS analysis on the modified system for the frequency range 0-1500 Hz is also listed in Table 1, marked as modified system. Compared to the results from the original system, a new resonance mode with negative damping factor (divergent) appears at a frequency of 29.8 Hz. In Figure 6, the nodal voltage mode shape shows that bus-A1 and bus-A2 are the dominant participation buses related to the unstable resonance mode. As bus-A1 and bus-A2 are the connecting point of the wind farm and the series capacitor, respectively, the previous expectation that a DFIG-based wind farm with a series-compensation line may cause SSR can be verified.  The wind farm is aggregated by DFIG-based wind generators with a total rated-capacity of 1500 MW. The detailed impedance model of a DFIG wind generator is given in reference [18]. The impedance-frequency characteristic of the aggregated wind farm for the frequency range 0-100 Hz is shown in Figure 5. It indicates that a negative resistance effect (∠Z DFIG > 90 • ) appears in the frequency range of 25-45 Hz, which may lead to SSR if the series resonance frequency also lies in this frequency range. Strictly speaking, the impedance model of a DFIG wind generator cannot be represented in the form of Y(s), because the outer loop controls (OLCs) of the d and q axes in the voltage source converter are not symmetric. However, the response time of the OLCs is usually larger than 200 ms, which means that the OLCs can be ignored for the considered frequency range studied. The impedance model in reference [18] is deduced based on the simplification above. The wind farm is aggregated by DFIG-based wind generators with a total rated-capacity of 1500 MW. The detailed impedance model of a DFIG wind generator is given in reference [18]. The impedance-frequency characteristic of the aggregated wind farm for the frequency range 0-100 Hz is shown in Figure 5. It indicates that a negative resistance effect (∠ > 90°) appears in the frequency range of 25-45 Hz, which may lead to SSR if the series resonance frequency also lies in this frequency range. Strictly speaking, the impedance model of a DFIG wind generator cannot be represented in the form of Y(s), because the outer loop controls (OLCs) of the d and q axes in the voltage source converter are not symmetric. However, the response time of the OLCs is usually larger than 200 ms, which means that the OLCs can be ignored for the considered frequency range studied. The impedance model in reference [18] is deduced based on the simplification above. The results of ENRS analysis on the modified system for the frequency range 0-1500 Hz is also listed in Table 1, marked as modified system. Compared to the results from the original system, a new resonance mode with negative damping factor (divergent) appears at a frequency of 29.8 Hz. In Figure 6, the nodal voltage mode shape shows that bus-A1 and bus-A2 are the dominant participation buses related to the unstable resonance mode. As bus-A1 and bus-A2 are the connecting point of the wind farm and the series capacitor, respectively, the previous expectation that a DFIG-based wind farm with a series-compensation line may cause SSR can be verified.  The results of ENRS analysis on the modified system for the frequency range 0-1500 Hz is also listed in Table 1, marked as modified system. Compared to the results from the original system, a new resonance mode with negative damping factor (divergent) appears at a frequency of 29.8 Hz. In Figure 6, the nodal voltage mode shape shows that bus-A1 and bus-A2 are the dominant participation buses related to the unstable resonance mode. As bus-A1 and bus-A2 are the connecting point of the wind farm and the series capacitor, respectively, the previous expectation that a DFIG-based wind farm with a series-compensation line may cause SSR can be verified. The wind farm is aggregated by DFIG-based wind generators with a total rated-capacity of 1500 MW. The detailed impedance model of a DFIG wind generator is given in reference [18]. The impedance-frequency characteristic of the aggregated wind farm for the frequency range 0-100 Hz is shown in Figure 5. It indicates that a negative resistance effect (∠ > 90°) appears in the frequency range of 25-45 Hz, which may lead to SSR if the series resonance frequency also lies in this frequency range. Strictly speaking, the impedance model of a DFIG wind generator cannot be represented in the form of Y(s), because the outer loop controls (OLCs) of the d and q axes in the voltage source converter are not symmetric. However, the response time of the OLCs is usually larger than 200 ms, which means that the OLCs can be ignored for the considered frequency range studied. The impedance model in reference [18] is deduced based on the simplification above. The results of ENRS analysis on the modified system for the frequency range 0-1500 Hz is also listed in Table 1, marked as modified system. Compared to the results from the original system, a new resonance mode with negative damping factor (divergent) appears at a frequency of 29.8 Hz. In Figure 6, the nodal voltage mode shape shows that bus-A1 and bus-A2 are the dominant participation buses related to the unstable resonance mode. As bus-A1 and bus-A2 are the connecting point of the wind farm and the series capacitor, respectively, the previous expectation that a DFIG-based wind farm with a series-compensation line may cause SSR can be verified.  Further sensitivity analysis of the newly-appeared unstable resonance mode is listed in Table 2. It shows that the stability (damping factor) is most influenced by the line resistance. Considering that adjusting R L is not applicable in real projects, applying a virtual impedance control at a certain frequency may be an optional choice to increase the system's stability. Table 2. Mode sensitivity analysis of the modified system.

Parameters
Sensitivity Apart from the new resonance mode, the existing resonance modes may be changed after the connection of the DFIG-based wind farm. However, the changes are negligible, hence they do not influence the system's stability.

Wind Farm Integration System with an MMC
In this subsection, the scenario of wind power integration through a modular multilevel converter (MMC) will be studied. The relevant case is shown in Figure 7, which is part of a newly planned renewable power exploitation project in China. In this case, four clustered wind power plants, with rated capacities of 400 MW, 100 MW, 300 MW and 400 MW, respectively, form an islanded sending system and are connected to the MMC. The wind power plants are assumed to have a DFIG structure (impedance characteristic is shown in Figure 5). Line parameters and MMC parameters are listed in Tables 3 and 4, respectively. Further sensitivity analysis of the newly-appeared unstable resonance mode is listed in Table 2. It shows that the stability (damping factor) is most influenced by the line resistance. Considering that adjusting RL is not applicable in real projects, applying a virtual impedance control at a certain frequency may be an optional choice to increase the system's stability. Apart from the new resonance mode, the existing resonance modes may be changed after the connection of the DFIG-based wind farm. However, the changes are negligible, hence they do not influence the system's stability.

Wind Farm Integration System with an MMC
In this subsection, the scenario of wind power integration through a modular multilevel converter (MMC) will be studied. The relevant case is shown in Figure 7, which is part of a newly planned renewable power exploitation project in China. In this case, four clustered wind power plants, with rated capacities of 400 MW, 100 MW, 300 MW and 400 MW, respectively, form an islanded sending system and are connected to the MMC. The wind power plants are assumed to have a DFIG structure (impedance characteristic is shown in Figure 5). Line parameters and MMC parameters are listed in Tables 3 and 4, respectively.     The impedance-frequency characteristic is acquired by using the test signal method [19]. Considering that the result of the test signal method is discrete, numerical fitting is executed so that the result can be used in the proposed method. Figure 8 shows the results of the test signal method and the fitting curves.  The impedance-frequency characteristic is acquired by using the test signal method [19]. Considering that the result of the test signal method is discrete, numerical fitting is executed so that the result can be used in the proposed method. Figure 8 shows the results of the test signal method and the fitting curves. The results of ENRS analysis on this system are listed in Table 5. As the system topology is relatively simple, only two resonance modes are found. The 126.6 Hz mode is sufficiently stable, while the 42.7 Hz mode is slightly divergent. Figure 9 shows the nodal voltage mode shape of the 42.7 Hz mode. Table 6 further shows its per-unit participation matrix (only the upper triangular components are listed considering the symmetry). Both Figure 9 and Table 6 indicate that this resonance mode appears significantly among all the six buses with similar phases. Namely, each connection point of the DFIG or MMC participates in, and is influenced by, this resonance mode. The mechanism of this kind of resonance is quite complex and the relevant suppression method will be demonstrated in future work that combine sensitivity analysis.  The results of ENRS analysis on this system are listed in Table 5. As the system topology is relatively simple, only two resonance modes are found. The 126.6 Hz mode is sufficiently stable, while the 42.7 Hz mode is slightly divergent. Figure 9 shows the nodal voltage mode shape of the 42.7 Hz mode. Table 6 further shows its per-unit participation matrix (only the upper triangular components are listed considering the symmetry). Both Figure 9 and Table 6 indicate that this resonance mode appears significantly among all the six buses with similar phases. Namely, each connection point of the DFIG or MMC participates in, and is influenced by, this resonance mode. The mechanism of this kind of resonance is quite complex and the relevant suppression method will be demonstrated in future work that combine sensitivity analysis.