Journal of Modern Power Systems and Clean Energy

ISSN 2196-5625 CN 32-1884/TK

网刊加载中。。。

使用Chrome浏览器效果最佳,继续浏览,你可能不会看到最佳的展示效果,

确定继续浏览么?

复制成功,请在其他浏览器进行阅读

Smallest Eigenvalues Based Logarithmic Derivative Method for Computing Dominant Oscillation Modes in Large-scale Power Systems  PDF

  • Linguang Wang 1 (Graduated Student Member, IEEE)
  • Xiaorong Xie 1 (Fellow, IEEE)
  • Wenkai Dong 1
  • Yong Mei 2
  • Aoyu Lei 2
1. State Key Laboratory of Power System, Department of Electrical Engineering, Tsinghua University, Beijing 100084, China; 2. Power Dispatching and Control Center of China Southern Power Grid, Guangzhou 510663, China

Updated:2025-05-21

DOI:10.35833/MPCE.2024.000630

  • Full Text
  • Figs & Tabs
  • References
  • Authors
  • About
CITE
OUTLINE

Abstract

With the rapid integration of renewable energy, wide-band oscillations caused by interactions between power electronic equipment and grids have emerged as one of the most critical stability issues. Existing methods are usually studied for local power systems with around one hundred nodes. However, for a large-scale power system with tens of thousands of nodes, the dimension of transfer function matrix or the order of characteristic equation is much higher. In this case, the existing methods such as eigenvalue analysis method and impedance-based method have difficulty in computation and are thus hard to utilize in practice. To fill this gap, this paper proposes a novel method named the smallest eigenvalues based logarithmic derivative (SELD) method. It obtains the dominant oscillation modes by the logarithmic derivative of the k-smallest eigenvalue curves of the sparse extended nodal admittance matrix (NAM). An oscillatory stability analysis tool is further developed based on this method. The effectiveness of the method and the tool is validated through a local power system as well as a large-scale power system.

I. Introduction

OSILLATION issues are of great concern in modern power systems with massive integration of renewable energy [

1]-[4]. The oscillations caused by interactions between power electronic converters and grids have posed severe threatens to the safe and stable operation of the system. Therefore, it is necessary to evaluate the risk of oscillation in a large-scale power system. For this purpose, the existing methods have been developed in either time domain [5]-[11] or frequency domain [12]-[18]. However, these methods have difficulties in dealing with large-scale power systems with tens of thousands of nodes. Their shortages are specifically introduced below.

The time-domain methods mainly include the eigenvalue analysis and electromagnetic transient (EMT) simulation [

5]-[11]. The former has the problem named “dimensional disaster”. The order of characteristic equation will become very high as the system expands, making it hard to calculate the eigenvalues accurately [19]. The latter also has an extremely high computational burden on numerical integration when there is a large amount of power electronics in the system.

The widely recognized frequency-domain methods are based on the impedance concept, for instance, the source-load method [

9]-[14], the frequency-domain modal analysis [15], [16], and the impedance aggregation method [17], [18]. The source-load method divides the system into converter side and grid side, and then uses the generalized Nyquist criterion (GNC) to judge the stability. However, this method has too many equivalences on the grid side and some oscillation modes may not be observable. It is normally used in the oscillation analysis of the local power system [14]. The frequency-domain modal analysis involves the construction of nodal admittance matrix (NAM) and the solution of oscillation modes. However, for a large-scale power system, the dimension of NAM would be too high to fulfil singular value decomposition (SVD) or matrix inversion operation [19]. The impedance aggregation method aggregates the impedances into a node. However, this method needs to calculate the inverse of the NAM in the aggregation of the impedances [17], which is very time-consuming for a high-dimensional NAM. Therefore, it is still an urgent necessity to develop an efficient method to evaluate oscillatory stability for extremely large-scale power systems. Furthermore, there is a need to develop analysis tools for actual systems.

To solve the problem of oscillatory stability analysis in practical large-scale power systems, this paper proposes a novel method named the smallest eigenvalues based logarithmic derivative (SELD) method. It has the following features:

1) For the modeling of large-scale power systems, the frequency-coupled admittance network model (FCANM) represented with sparse extended NAM is used to improve the performance in storage.

2) The implicit restarted Arnoldi method (IRAM) is applied to calculate the k-smallest eigenvalues of the sparse matrix with the great improvement in time complexity of computation, thus making the proposed method applicable to extremely large-scale power systems.

3) A segmentation function based logarithmic derivative criterion is proposed to identify the dominant oscillation modes based on the curves of the k-smallest eigenvalues.

An oscillatory stability analysis tool is further developed to implement the proposed method for an actual large-scale power system. The effectiveness of the proposed method and tool is verified in a local power system and a large-scale power system, respectively.

The rest of the paper is organized as follows: Section II introduces the problem of oscillation analysis in a practical large-scale power system. Section III proposes the SELD method for analyzing the oscillatory stability of large-scale power systems and develops an oscillation stability analysis tool. Section IV validates the effectiveness of the method and tool in two practical systems with different sizes. Section V gives a brief conclusion.

II. Problem Statement

A. Description of Large-scale Power System

The oscillatory stability of a practical large-scale power system in China Southern Power Grid is analyzed. The diagram of the large-scale power system is shown in Fig. 1. This system includes 4 regional power systems (Yunnan grid, Guizhou grid, Guangxi grid, Guangdong grid). The numbers of buses in Yunnan grid, Guizhou grid, Guangxi grid, and Guangdong grid are 5109, 1638, 2028, and 5203, respectively. Yunnan Grid is asynchronously connected to other grids by high-voltage direct current (HVDC) transmissions. Other grids are interconnected by HVDC and high-voltage alternative current (HVAC) transmissions. The renewable energy sources are mainly integrated into Yunnan grid and Guangdong grid. The total capacities of renewable energy sources in Yunnan grid and Guangdong grid are approximately 8000 MVA and 10000 MVA, respectively.

Fig. 1  Diagram of a large-scale power system.

The order of the characteristic equation becomes extremely high due to the large number of electrical equipment. The total numbers of nodes and branches are 14647 and 20309, respectively. There are approximately 6872 transmission lines, 11804 transformers, and 1459 thermal units or synchronous condensers in this system. Besides, there are around 86 wind farms or photovoltaic power plants in the system. The total capacity of renewable energy sources is about 18000 MVA. The transmission lines are considered as a π-equivalence model, which is a three-order circuit. The transformers can be equated to the inductors. The thermal units or synchronous condensers can be modeled as an 18-order state space model. For simplification, each wind farm is aggregated as the wind turbine generators (WTGs) in parallel connection. The order of the WTG model is 8 at least [

20]. Therefore, the total order of the characteristic equation in this system can be estimated as 6872×3×2+11804×1×2+1459×18+86×18=92650, which is extremely large. Note that since the frequency-coupled admittance model (FCAM) is applied, the “×2” is needed for transmission lines and transformers. The high order of the system obstacles to oscillatory stability analysis.

B. Computational Problem in Existing Analysis Methods

For the eigenvalue analysis method, according to the estimated order of the characteristic equation, the dimension of matrix A of the state space equation is approximately 105. The calculation of the eigenvalue of matrix A is impossible to achieve, as the eigenvalue computation requires O(n3) floating point operations (FLOPs), where O denotes the upper bound of a function within a constant factor; and n is the dimension of matrix A [

19]. Thus, the computational effort of eigenvalue analysis method is unbearable for analyzing the large-scale power system in Fig. 1.

Conventional impedance-based methods analyze the stability by the closed-loop transfer function matrix. There are mainly two schemes. The first scheme uses the stability theory of the multiple-input multiple-output (MIMO) system [

18]. The closed-loop transfer function matrix of the system is denoted as Y^(s), and the oscillation modes are the zero elements of det(Y^(s)), where det() denotes the determinant. The second scheme uses the MIMO system of the voltages and currents of an aggregated node [17]. The aggregated impedance ZΣ is defined and obtained by (1) and (2). The oscillation modes are the zero elements of the det(ZΣ(s)). For simplicity, (s) is omitted.

Ui=ZΣIi (1)
ZΣ=(Y^-1)ii (2)

where i is the node number chosen to aggregate; Ui and Ii are the voltage and current with coupled frequencies of node i, respectively; and subscript ii means the ith row and the ith column of the matrix by 2×2 partition.

However, both schemes have difficulties in computation. For the first one, the algorithm of determinant calculation is normally based on the LU decomposition. The determinant is the production of diagonal elements of L and U, where L and U are the decomposed lower and upper triangular matrices, respectively. The LU decomposition is based on Gaussian elimination. The transformation of the original matrix into an upper triangular matrix via the Gaussian elimination needs O(n3) FLOPs [

19]. For the large-scale power system, n is normally very large, O(n3) FLOPs are hard to achieve for computers. For the second one, the calculation of inverse of matrix also needs O(n3) FLOPs, which is also very time-consuming.

Therefore, it is necessary to develop a method with high computational efficiency to analyze oscillatory stability for extremely large-scale power systems.

III. Oscillatory Stability Analysis by SELD Method

A. FCANM Represented with Sparse Extended NAM

According to the impedance-based method, each component in a power system can be modeled as an FCAM represented by a 2×2 matrix [

21], which is described as:

ΔIp(ω)ΔIpc*(ω)=Y11(ω)Y12(ω)Y21(ω)Y22(ω)ΔUp(ω)ΔUpc*(ω) (3)

where ΔUp, ΔUpc, ΔIp, and ΔIpc are the small-signal voltage and current phasors of two coupled frequencies fp and 2f1-fp, respectively, and f1 and fp are the fundamental and non-fundamental frequencies, respectively; superscript * means the conjugation of phasor; ω is the angular frequency of fp; and Yij is the element of FCAM.

The FCANM can be described by the extended NAM Y^(ω). Then, the network is expressed as:

I^n(ω)=Y^(ω)U^n(ω) (4)

where U^n(ω) is the nodal voltage vector; and I^n(ω) is the nodal current vector.

The FCANM can be represented by:

Y^(ω)=D^Y^D(ω)D^T (5)
D^=DI2 (6)

where Y^D(ω) is the diagonal matrix with each diagonal element representing the FCAM; D^ is the extended node-branch incidence matrix; I2 is a two-dimensional identity matrix; denotes the Kronecker product; and D is the node-branch incidence matrix generated by the system topology.

The huge number of AC nodes in a large-scale power system would result in problems in storage and calculation of FCANM. The cost of the matrix multiplication in (5) is O(nf(b2nv+bnv2)) FLOPs, where nv is the number of nodes; b is the number of branches; and nf is the number of frequency points. It is extremely time-consuming for a large-scale power system.

The sparse matrix technique is very effective in addressing this issue. Similar to the power flow calculation, there are very few non-zero elements in D^ (e.g., 0.1% non-zeros elements for a practical system). Besides, Y^D(ω) is a diagonal matrix. Therefore, both the matrices are sparse.

In (5), the operation of the zero elements in the matrix is not necessary. Meanwhile, the storage of these zero elements is also redundant. The sparse matrix uses three vectors u, v, w to store and operate the non-zero elements in the matrix only, where u, v, and w are the vectors of horizontal coordinates, vertical coordinates, and values of non-zero elements, respectively. As a result, the massive storage and calculation of zero elements are eliminated. Using sparse matrix can accelerate the calculation by tens or even hundreds of times. The memory requirement can also be greatly reduced. The larger the scale of the power system, the more significant the benefits. Therefore, the computational efficiency in constructing FCANM is highly improved by the sparse extended NAM, which is the extended NAM Y^(ω) with sparse matrix representation.

B. Improvement of Computational Efficiency Using k-smallest Eigenvalues

To improve the computational efficiency, the k-smallest eigenvalues of the matrix Y^ rather than the determinant or the inverse of the matrix is used, where Y^ is Y^(ω) at any ω. The k-smallest eigenvalues are the k eigenvalues with the smallest modulus. By using them, the computational efficiency of oscillatory stability analysis is improved from O(n3) to O(n) FLOPs, which will be specifically introduced in Section III-C. The reason for selecting the k-smallest eigenvalues is introduced as follows.

The SVD of Y^ is given by:

Y^=TTΛT (7)
Λ=λ1000λ2000λn (8)

where Λ and T are the orthogonal matrices, and the rows of Λ and the columns of T are the eigenvectors; and λi is the eigenvalue.

For λi on the diagonal of Λ, these exists:

λi-1xiI^n=xiU^n (9)

where xi is the eigenvector of λi.

The modal impedance is defined as λi-1. Each peak value of the modal impedance is usually corresponding to an oscillatory mode [

15], [16]. Meanwhile, the peak values of modal impedance are usually accompanied with small eigenvalues λi. Therefore, the oscillation modes are more easily identified by the small eigenvalue curves λi(ω). To balance both the accuracy and the computational efficiency, the k-smallest eigenvalue curves among those eigenvalue curves are selected and the frequency-domain stability criterion is applied on them.

C. Calculation of k-smallest Eigenvalues

The IRAM is applied to calculate the k-smallest eigenvalues of the n-dimensional matrix A. The key idea of IRAM is to decrease the dimension. The specific process is introduced as follows.

1) Krylov Subspace

Define the Krylov subspace of the n-dimensional matrix A as:

Km(A,x)=span{x,Ax,...,Am-1x} (10)

where x is a n-dimensional vector; m is an integer; and span{} is the set of all possible vectors with linear combinations of the vectors.

2) Arnoldi Process

The Arnoldi process is an efficient way to reduce the dimension of matrix A, thereby improving the computational efficiency. It transfers the n-dimensional matrix A to a m-dimensional matrix Hm by Krylov subspace, where m<<n.

The first step of the Arnoldi process is to find the orthonormal basis q1,q2,,qm (11) of the Krylov subspace (10). The Gram-Schmidt process is used in this step, i.e.,

span{x,Ax,...,Am-1x}=span{q1,q2,...,qm} (11)
q1=xx (12)
qi+1=1hi+1,iAqi-j=1iqjhj,i (13)

where hj,i is the Gram-Schmidt coefficient.

hj,i=qjTAqi (14)

The Gram-Schmidt coefficients form the (m+1)×m upper Hessenberg matrix H¯m.

H¯m=h1,1h1,2h1,m-1h1,mh2,1h2,2h2,m-1h2,m0h3,2h3,m-1h3,m000hm+1,m (15)
AQm+1=Qm+1H¯m (16)

where Qm+1=[q1,q2,,qm+1].

The transformation of Arnoldi process can be written as:

AQm=QmHm+hm+1,mqm+1emT (17)

where Hm is H¯m without the last row; and emT=[0,0,,0,1].

AQmuiQmHmui=λiQmui (18)
vi=Qmui (19)

According to (18), the eigenvalues {λi} and eigenvectors {ui} of Hm are calculated to approximate the eigenvalues {λi} and eigenvectors {vi} of A. This approximation holds when the term hm+1,mqm+1emTui is small or:

hm+1,mqm+1emTuihm+1,mui,m<ε (20)

where ε is the tolerance; and ui,m is the mth component of ui.

3) Procedure of IRAM

In the Arnoldi process, the restriction (20) does not usually hold. Therefore, the restart process is used to iterate the initial value of the algorithm until (20) holds. The IRAM improves the Arnoldi process by adding the restart process and using the implicit QR decomposition.

The focus is on the k-smallest eigenvalues, while leaving p=m-k eigenvalues ignored. The Arnoldi process begins with a unit vector v and computes the m step Arnoldi decomposition (13). Then, the eigenvalues and eigenvectors of Hm are calculated by the reduction to the upper Hessenberg form. The eigenvalues are sorted from the largest to the smallest by the magnitude. The m-k largest eigenvalues are used as a shift to enhance the calculation of the required k-smallest eigenvalues. The restart is chosen by the filtering polynomial pr(Hm) in (21). The pr(Hm) shifts away from the undesirable m-k eigenvalues, as opposed to shifts closer to the desirable k eigenvalues.

pr(Hm)=(Hm-λk+1I)(Hm-λk+2I)(Hm-λmI) (21)

where I is the m-dimensional identity matrix.

Apply the following implicit QR decomposition to these m-k eigenvalues.

Hm-λiI=QiRi (22)

where Qi is the orthogonal matrix; and Ri is the upper triangle matrix.

The Hm is iterated to H˜m by:

H˜m=Q˜THmQ˜ (23)
Q˜=i=k+1mQi (24)

Then, the restart initial vector v˜ is derived as [

19]:

v˜=Qk+1H˜m,k,k+1+v¯Q˜m,k (25)

where Qk+1 is the (k+1)th column of Qm in the previous iteration; H˜m,k,k+1 is the element at the kth row and (k+1)th column of H˜m; v¯ is the previous initial vector; and Q˜m,k is the element at the mth row and kth column of Q˜.

We use v˜ as the new initial vector and apply the Arnoldi process to obtain the eigenvalues and eigenvectors. If (20) holds for the k-smallest eigenvalues, return the k-smallest eigenvalues. Otherwise, repeat the above restart process.

4) Computational Efficiency

Since the Arnoldi process uses the Gram-Schmidt process, the approximate number of FLOPs for the algorithm is 2m2n [

19]. In the procedure of IRAM, the filter polynomial (21) is evaluated by the implicit QR decomposition. The implicit QR decomposition only costs O(m2) FLOPs. Therefore, the IRAM requires O(m2n) FLOPs. Since m is a very small value in the case of calculating the k-smallest eigenvalues, O(m2n) is equivalent to O(n), which is a huge improvement comparing to O(n3) FLOPs in the LU decomposition. For the oscillatory stability analysis in an extremely large-scale power system, the time complexity of this method O(n) is sufficient.

D. Identification of Dominant Oscillation Modes by Segmentation Function Based Logarithmic Derivative Criterion

The dominant oscillation modes of a power system are identified by the logarithmic derivative of the k-smallest eigenvalue curves. Due to the utilization of the small-signal linearized model, the general form of eigenvalue curves λi(ω) takes the following expression:

λi(ω)=Ki=1Njω-zii=1Mjω-pi (26)

where zi=αzi+jωzi and pi=αpi+jωpi are the zeros and poles of the λi(s), respectively; N and M are the numbers of zeros and poles, respectively; and K is a constant coefficient.

The jth smallest eigenvalue λminj(ω) (j=1,2,,k) is the segmentation function of λi(ω). Except for the non-differentiable points, λminj(ω) follows (26). The logarithmic derivative criterion is used to identify the dominant zeros of λminj(ω).

The logarithmic derivative of a function f(ω) is defined as the derivative of the logarithm of the function (27). Its differential form is given in (28), which is more pragmatic.

DL(f)=dIn(f)d(jω)=f'(ω)jf(ω) (27)
DL(f)=1jΔωf(ω+Δω)-f(ω)f(ω) (28)

where f'(ω) is the derivative of f(ω).

The logarithmic derivative of λi(ω) in (26) is presented as:

DL(λi)=i=1N1jω-zi-i=1M1jω-pi (29)

From (29), it can be observed that DL(λi) is the sum of DL(Gzi) minus the sum of DL(Gpi). DL(Gzi) and DL(Gpi) have the same expression, as shown in (30). The real- and imaginary-part curves of DL(Gzi) are shown in Fig. 2.

DL(Gzi)=1jω-ziDL(Gpi)=1jω-pi (30)

Fig. 2  Real- and imaginary-part curves of DL(Gzi). (a) Real-part curve. (b) Imaginary-part curve.

According to (30) and Fig. 2, DL(Gzi) and DL(Gpi) have the following properties.

1) The real part of the logarithmic derivative curve is larger only near the oscillation frequency ωi and decays sharply at the inverse square rate away from ωi. ωi is the extreme point whose value of the real part of the logarithmic derivative is the negative inverse of the real part of zi.

Re(DL(Gzi))ωi=-1αi (31)

where αi is the divergence rate.

2) The imaginary part of the logarithmic derivative curve is approximately symmetric about the center of the point (ωi,0), and again has a large absolute value only near ωi, while it decays at an inverse rate away from ωi. As the frequency increases, the imaginary part of the logarithmic derivative curve rises, then falls, then rises again and has a negative slope -1/αi2 at ωi.

The properties of the real- and imaginary-part curves of the DL(λi) can be obtained by (29) and the properties of DL(Gzi) and DL(Gpi). Suppose there is a dominant zero zi=αi+jωi.

Firstly, since the real part of DL(Gzi) decays sharply at the inverse square rate away from ωi, the influence of other zero-poles on the real part of DL(λi) can be approximately neglected near the frequency ωi of the dominant zero zi.

Secondly, due to that the real-part curve has an extreme value on ωi, the stationary points of the real part of DL(λi) can be used to identify the dominated zeros.

Thirdly, based on the surge of the imaginary-part curve at ωi, the slope of the imaginary part of the logarithmic derivative curve can be used to distinguish the zeros from the poles. The slope of the imaginary part of the zero (pole) term at the frequency ωi is negative (positive).

Finally, the real part of the logarithmic derivative is solved for the real part of the zero-pole based on the extreme value of the real part of the logarithmic derivative [

22].

According to the properties of the real- and imaginary-part curves of DL(λi), the logarithmic derivative criterion is given as follows.

If Re(DL(λi)) exhibits as a stationary point (Re(DL(λi))'=0) at the oscillation frequency ω and the slope of Im(DL(λi)) at ω is negative, there exists a zero z=α+jω, and the real part α can be solved by:

α=-1Re(DL(λi))ω (32)

Examine all the identified zeros. If there is a zero z=α+jω with a positive real part, it corresponds to an unstable oscillation mode. Conversely, if the real parts of all the obtained zeros are less than 0, the system is considered stable.

For those non-differentiable points, both the real- and imaginary-part curves charge abruptly at this point, thus the slope cannot be calculated, and the above criterion can ignore those points.

E. Oscillatory Stability Analysis Tool Based on SELD Method

1) Process of Oscillatory Stability Analysis Tool

A platform of data collection and its interfaces are developed. This platform collects the original data of the power system. The original data include both the static and dynamic data. The static data include the parameters of transmission lines, WTGs, thermal units, transformers, etc. The dynamic data include wind speed, system load, operating condition of power system, etc. These data are transformed into standard data files, which are the inputs of the tool.

The process of the oscillatory stability analysis tool is given in Fig. 3. The original data are collected to obtain the parameters of electrical equipment, initial value of power flow, and the system topology. The power flow is firstly calculated by these data to determine the operating condition of the electrical equipment. Wide-band FCAMs are established based on the parameters of electrical equipment and the power flow. Then, the FCANM represented with sparse matrix is built by the FCAMs. Finally, the proposed method is applied to identify the dominant oscillation modes.

Fig. 3  Process of oscillatory stability analysis tool.

2) Algorithm

The algorithm of the proposed method in the oscillatory stability analysis tool is given in Algorithm 1.

Algorithm 1  : proposed SELD method

Input: data collected by the platform.

1. Initialize the array of “branch” and “node” by the data. The parameters of the equipment xj (j=1,2,,b), initial value of power flow, and node-branch incidence matrix D are acquired.

2. The nodal Ui, δi, Pi, Qi (i=1,2,,n) are calculated by the power flow analysis.

3. The FCAMs Yj of each branch are obtained by the parameters xj and the nodal Ui, δi, Pi, Qi.

4. The sparse-matrix-based FCANM Ŷ(ω) is constructed by (5).

5. The k-smallest eigenvalue curves λl(ω) are obtained by IRAM, where l=1,2,,k.

6. The dominant oscillation modes zi=αi+jωi are identified by DL(λl(ω)) (28), (32).

Output: the dominant oscillation modes zi of the system.

F. A Demonstrative Example for Proposed Method

A demonstrative example of two parallel RLC circuits, as shown in Fig. 4, is used to help understand the implementation of the proposed method. Obviously, the circuit has two parallel resonance modes 1 and 2 formed by L1, C1 and L2, C2, respectively.

Fig. 4  Demonstrative example of two parallel RLC circuits.

The NAM is derived as (33). The circuit parameters are: R1=1 Ω, L1=1 H, C1=10 F, R2=1 Ω, L2=0.1 H, C2=5 F, R0=100 Ω. The resonance modes are: -0.1010±j0.2245×2π and -0.0505±j0.0497×2π.

Y(s)=1R1+1sL1+sC1+1R0-1R0-1R01R2+1sL2+sC2+1R0 (33)

The eigenvalue of Y(ω) is calculated. The two modal impedance curves 1 and 2 are obtained by calculating the inverse of absolute value of eigenvalues. The extreme points of modal impedance curves indicate the oscillation modes. The curve of the largest modal impedance is selected to contain the extreme points. The modes 1 and 2 are clearly identified by the peak value of the curve of the largest modal impedance in Fig. 5. The curve of the largest modal impedance also corresponds to the curve of the smallest eigenvalue. Therefore, the smallest eigenvalue curve contains more information about the oscillations and should be selected. The IRAM method is used to calculate the k-smallest eigenvalues with high computational efficiency.

Fig. 5  Modal impedance curves.

The smallest eigenvalue curve has the expression of (26). The logarithmic derivative criterion is adopted to the smallest eigenvalue curve to obtain the oscillation modes. The logarithmic derivative curves obtained by the proposed method are shown in Fig. 6. It is shown that there are two stable modes at the frequency of 0.0497 Hz and 0.02244 Hz, respectively, which is consistent with results of circuit analysis.

Fig. 6  Logarithmic derivative curves obtained by proposed method.

IV. Applications in Two Practical Systems

Two practical systems are used to validate the effectiveness of the proposed method and the developed tool. A local power system is used to verify the correctness. The large-scale power system in Fig. 1 with 14000 nodes is used to validate the computational efficiency.

A. Local Power System

1) Local Power System Modeling

The simplified diagram of a real offshore wind-integrated local power system is shown in Fig. 7.

Fig. 7  Simplified diagram of a real offshore wind-integrated local power system.

The WTGs in these wind farms are all the same type of direct-drive permanent magnet synchronous generator (PMSG), with a total capacity of 2734.2 MVA. The local load in this area is very small. The power from QZ, YY, and PT wind farms is collected to the 500 kV PL station via 220 kV transmission lines, and then to the EH2 station via a 500 kV transmission line (EH2-PL). The power from QZ6 wind farm is collected to the 500 kV EH1 station via several transmission lines. The EH1 station and EH2 station are connected to the external network equivalence via 500 kV transmission lines. This local power system is equipped with five thermal units, all of which have a rated capacity of 1112 MVA. The BH1-3 thermal power plant is equipped with three thermal units. The BH4 thermal power plant and the YX thermal power plant are each equipped with one thermal unit. The numbers of nodes and branches are 64 and 93, respectively, in this local power system.

The FCAMs of components in this power system are established based on the collected parameters. The FCAMs of the PMSGs are derived in [

20]. The FCAM of each wind farm is created by aggregating the FCAMs of the PMSGs. The FCAMs of thermal units are derived by the state-space equations. Additionally, other components such as transmission lines and transformers are also modeled as FCAMs. The external system is represented as voltage-behind-reactance model based on the short-circuit capacities at the WL and CJ stations, respectively. The sparse-matrix-based impedance network model is constructed in Section III-A.

2) Oscillatory Stability Analysis in Local Power System

To verify the effectiveness of the proposed method, eigenvalue analysis, zero-crossing stability criterion [

17], and for comparison purposes, the proposed method are applied.

With the derived FCAMs of components, the s-domain FCANM is constructed by (5) and next transferred into state-space equations. Then, the eigenvalue analysis can be carried out. It is observed that there is an eigenvalue with a positive real part, 2.7497+j46.23×2π. This indicates an unstable oscillation mode at the frequency of 46.23 Hz.

Apply the zero-crossing stability criterion to the extended NAM det(Y^(ω)). Its real and imaginary parts with the frequency around 46 Hz are plotted in Fig. 8. The real-part curve has a zero-crossing point at the frequency of 46.24 Hz. The product of the slope at the zero-crossing point and the value of the imaginary part at this frequency is larger than 0. This analysis reveals the presence of an unstable oscillation mode characterized by an oscillation frequency of 46.24 Hz.

Fig. 8  Real and imaginary parts of extended NAM.

Next, the proposed method is employed and the smallest eigenvalue λmin(ω) is calculated by the IRAM. The real and imaginary parts of the logarithmic derivative within the frequency range of 45.5-47 Hz are plotted in Fig. 9.

Fig. 9  Real and imaginary parts of logarithmic derivative within frequency range of 45.5-47 Hz.

There are several stationary points within the colored intervals of the real-part curve. Since the slope of the imaginary-part curve is positive in the blue interval, the stationary point corresponds to a pole of the transfer function. The stationary points in the green intervals are the non-differentiable points of the smallest eigenvalue and cannot be identified by the criterion in Section III-D. Since the slope of the imaginary-part curve is negative in the red interval, the stationary point corresponds to a dominant oscillation mode. By calculating the real part of this mode (32), an unstable mode is identified, and the oscillation frequency is 46.31 Hz.

Comparing the results obtained from the proposed method with the zero-crossing stability criterion, the relative error of the latter in the frequency is 0.15%. Such small relative errors demonstrate the accuracy of the proposed method, as the results closely align with those obtained from the zero-crossing stability criterion.

The proposed and existing methods were also validated in other scenarios, and the results verify the effectiveness of the proposed method as well.

3) Validation by Time-domain Simulation

The accuracy of the proposed method is further validated by comparing its results with commercial EMT simulation. When the output reactive power of the WTGs is increased to 0.3 p.u., sub-/super-synchronous oscillation is excited. The output active power of the WTGs is plotted in Fig. 10.

Fig. 10  Output active power of WTGs.

The frequency spectrum of the output current in the sub-/super-synchronous band are presented in Fig. 11, clearly indicating the presence of distinct sub-/super-synchronous oscillations with frequencies of 46.27 Hz and 53.73 Hz. Comparing these results, it is observed that the oscillation frequencies obtained through the EMT simulation are consistent with the results of the proposed method, which further confirms the accuracy and reliability of the proposed method in capturing and analyzing the oscillatory behavior of the system under investigation.

Fig. 11  Frequency spectrum of output current in sub-/super-synchronous band.

B. Large-scale Power System

1) Large-scale Power System Modeling

The large-scale power system in Fig. 1 is modeled for the oscillation analysis.

Firstly, the original data are collected using the dynamic simulation program (DSP) developed by Electric Power Research Institute of China Southern Power Grid. The data are converted into the standard format by our developed tool.

Secondly, the FCAMs are constructed by the data. Due to the large number and type of new energy units in the actual power system, it is difficult to obtain all the new energy unit models from the manufacturers. For simplicity, we use typical models of new energy units in the analysis. The direct-driven PMSGs and photovoltaic units are equivalent by the admittance model derived in [

20]. The utilized FCAMs of both the doubly-fed induction generators (DFIGs) and the thermal units are derived by the state-space equations. The transformers are modeled as resistance and reactance in series connection by the leakage reactance and losses. The transmission lines are treated as the π-equivalent circuit.

Then, the extended node-branch incidence matrix D^ is obtained by the system topology and stored by the sparse matrix. The diagonal matrix Y^D(ω) is obtained by the FCAMs of the branches and stored by the sparse matrix.

Finally, the FCANM is constructed by (5), and the k-smallest eigenvalue curves are obtained by IRAM (k=10).

2) Oscillatory Stability Analysis in Large-scale Power System

The k-smallest eigenvalue curves are used to identify the dominant zeros by the improved logarithmic derivative criterion. Figure 12 displays the logarithmic derivative curves of the 4th smallest eigenvalue curve in the frequency range of 45-49 Hz. By applying the evaluation criterion of the proposed method, the oscillation mode is identified as unstable and the oscillation frequency is identified as 46.94 Hz.

Fig. 12  Logarithmic derivative curves of the 4th smallest eigenvalue curve.

3) Computational Efficiency

The above results indicate that the oscillation analysis in the extremely large-scale power system is achievable. The total time consumed for the proposed method is 410.57 s.

As a comparison, the existing methods are applied to this large-scale power system. According to the analysis in Section II-B, the computational efficiency of both the eigenvalue analysis method and conventional impedance-based method is very low. The number of FLOPs of the eigenvalue analysis method is (105)3=1015 at least, where 105 is the estimated order of the characteristic equation. The CPU used has a base frequency of 2.90 GHz, so the computational time is 1015/2.90/109/3600=95.8 hours, which is too long, not to mention the huge amount of time consumed in memory access. As for the conventional impedance-based method, assuming that the number of selected frequency points is 40, the estimated computational time is 40×(14000×2)3/2.90/109/ 3600=84.1 hours, which is also too long. In practical applications, the eigenvalue analysis method and the conventional impedance-based method are unable to produce results in 12 hours. This indicates that the proposed method has an advantage over the above two methods in terms of computational efficiency.

The aforementioned computation was performed on an Intel Core i7-10700 CPU with a base frequency of 2.90 GHz, 16.0 GB RAM, and Windows 10 operating system.

V. Conclusion

This paper proposes an SELD method to analyze the oscillatory stability of the large-scale power systems. Compared with the existing methods, the proposed method has following advantages.

1) The sparse extended NAM based modeling method has better performance in storage than the conventional method in the modeling of large-scale power systems.

2) The k-smallest eigenvalues of the sparse extended NAM are used to represent the dynamics of the power system. The IRAM significantly improves the time complexity of computation from O(n3) to O(n).

3) The segmentation function based logarithmic derivative criterion is proposed to quickly and accurately identify the dominant oscillation modes based on the k-smallest eigenvalues.

The accuracy and computational efficiency of the proposed method and the developed oscillatory stability analysis tool are demonstrated through case studies on two practical systems. In the case of a local power system, the results of the proposed method are equivalent to those of the impedance-based method and EMT simulations. In the case of a large-scale power system with 14000 nodes, the total time consumed for the proposed method is 410.57 s. The proposed method has better computational efficiency over the eigenvalue analysis method and the conventional impedance-based method.

References

1

J. Wachter, L. Gröll, and V. Hagenmeyer, “Survey of real-world grid incidents – opportunities, arising challenges and lessons learned for the future converter dominated power system,” IEEE Open Journal of Power Electronics, vol. 5, pp. 50-69, Dec. 2024. [Baidu Scholar] 

2

N. Modi, E. M. Farahani, A. Jalali et al., “Replication of real-world sub-synchronous oscillations in inverter-based resources dominated grid,” IEEE Transactions on Power Delivery, vol. 39, no. 3, pp. 1399-1406, Jun. 2024. [Baidu Scholar] 

3

L. Fan, Z. Miao, and Z. Wang, “A new type of weak grid IBR oscillations,” IEEE Transactions on Power Systems, vol. 38, no. 1, pp. 988-991, Jan. 2023. [Baidu Scholar] 

4

Q. Chen, S. Bu, and C. Y. Chung, “Small-signal stability criteria in power electronics-dominated power systems: a comparative review,” Journal of Modern Power Systems and Clean Energy, vol. 12, no. 4, pp. 1003-1018, Jul. 2024. [Baidu Scholar] 

5

L. P. Kunjumuhammed, B. C. Pal, R. Gupta et al., “Stability analysis of a PMSG-based large offshore wind farm connected to a VSC-HVDC,” IEEE Transactions on Energy Conversion, vol. 32, no. 3, pp. 1166-1176, Sept. 2017. [Baidu Scholar] 

6

J. Liu, W. Yao, J. Wen et al., “Impact of power grid strength and PLL parameters on stability of grid-connected DFIG wind farm,” IEEE Transactions on Sustainable Energy, vol. 11, no. 1, pp. 545-557, Jan. 2020. [Baidu Scholar] 

7

S. Bu, W. Du, and H. Wang, “Model validation of DFIGs for power system oscillation stability analysis,” IET Renewable Power Generation, vol. 11, no. 6, pp. 858-866, May 2017. [Baidu Scholar] 

8

Z. Zhang, X. Zhao, L. Fu et al., “Stability and dynamic analysis of the PMSG-based WECS with torsional oscillation and power oscillation damping capabilities,” IEEE Transactions on Sustainable Energy, vol. 13, no. 4, pp. 2196-2210, Oct. 2022. [Baidu Scholar] 

9

W. Zhou and J. Beerten, “Black box-based incremental reduced-order modeling framework of inverter-based power systems,” IEEE Open Journal of the Industrial Electronics Society, vol. 4, pp. 506-518, Nov. 2023. [Baidu Scholar] 

10

C. Collados-Rodriguez, M. Cheah-Mane, E. Prieto-Araujo et al., “Stability analysis of systems with high VSC penetration: where is the limit?” IEEE Transactions on Power Delivery, vol. 35, no. 4, pp. 2021-2031, Aug. 2020. [Baidu Scholar] 

11

Y. Han, H. Sun, B. Huang et al., “Discrete-time domain modal analysis of oscillatory stability of renewables integrated power systems,” IEEE Transactions on Power Delivery, vol. 37, no. 5, pp. 4248-4260, Oct. 2022. [Baidu Scholar] 

12

J. Sun, “Impedance-based stability criterion for grid-connected inverters,” IEEE Transactions on Power Electronics, vol. 26, no. 11, pp. 3075-3078, Nov. 2011. [Baidu Scholar] 

13

B. Wen, D. Boroyevich, R. Burgos et al., “Analysis of D-Q small-signal impedance of grid-tied inverters,” IEEE Transactions on Power Electronics, vol. 31, no. 1, pp. 675-687, Jan. 2016. [Baidu Scholar] 

14

W. Cao, Y. Ma, and F. Wang, “Sequence-impedance-based harmonic stability analysis and controller parameter design of three-phase inverter-based multibus AC power systems,” IEEE Transactions on Power Electronics, vol. 32, no. 10, pp. 7674-7693, Oct. 2017. [Baidu Scholar] 

15

W. Xu, Z. Huang, Y. Cui et al., “Harmonic resonance mode analysis,” IEEE Transactions on Power Delivery, vol. 20, no. 2, pp. 1182-1190, Apr. 2005. [Baidu Scholar] 

16

S. Chou, X. Wang, and F. Blaabjerg, “Frequency-domain modal analysis for power-electronic-based power systems,” IEEE Transactions on Power Electronics, vol. 36, no. 5, pp. 4910-4914, May 2021. [Baidu Scholar] 

17

H. Liu, X. Xie, and W. Liu, “An oscillatory stability criterion based on the unified dq-frame impedance network model for power systems with high-penetration renewables,” IEEE Transactions on Power Systems, vol. 33, no. 3, pp. 3472-3485, May 2018. [Baidu Scholar] 

18

L. Wang, X. Xie, J. Shair et al., “Frequency-domain admittance network model (FANM) based oscillatory stability analysis of hybrid AC-DC systems with MMCs,” IEEE Transactions on Power Systems, vol. 39, no. 2, pp. 3444-3458, Mar. 2024. [Baidu Scholar] 

19

W. H. Ford, Numerical Linear Algebra with Applications: Using MATLAB. San Diego: Academic Press, pp. 541-543, 2015. [Baidu Scholar] 

20

W. Liu, X. Xie, J. Shair et al., “A nearly decoupled admittance model for grid-tied VSCs under variable operating conditions,” IEEE Transactions on Power Electronics, vol. 35, no. 9, pp. 9380-9389, Sept. 2020. [Baidu Scholar] 

21

W. Liu, X. Xie, X. Zhang et al., “Frequency-coupling admittance modeling of converter-based wind turbine generators and the control-hardware-in-the-loop validation,” IEEE Transactions on Energy Conversion, vol. 35, no. 1, pp. 425-433, Mar. 2020. [Baidu Scholar] 

22

L. Wang, X. Xie, T. Xu et al., “Oscillatory stability analysis method based on the logarithmic derivative calculation of aggregated impedance determinant,” Proceedings of the CSEE, vol. 43, no. 11, pp. 4254-4261, Jun. 2023. [Baidu Scholar]