A Closed-form Formulation of Eigenvalue Sensitivity Based on Matrix Calculus for Small-signal Stability Analysis in Power System

Peijie Li; Yucheng Wei; Junjian Qi; Xiaoqing Bai; Hua Wei

网刊加载中。。。

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

确定继续浏览么?

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

A Closed-form Formulation of Eigenvalue Sensitivity Based on Matrix Calculus for Small-signal Stability Analysis in Power System PDF

- ORCID：
Peijie Li
✉
- ORCID：
Yucheng Wei
✉
- ORCID：
Junjian Qi
✉
- ORCID：
Xiaoqing Bai
✉
- ORCID：
Hua Wei
✉

School of Electrical Engineering, Guangxi University, Nanning 530004, China； Department of Electrical and Computer Engineering, Stevens Institute of Technology, Hoboken, NJ 07030, USA

Updated：2021-11-23

DOI：10.35833/MPCE.2019.000225

OUTLINE

Abstract

With the rapid development of power-electronics-enabled power systems, the new converter-based generators are deteriorating the small-signal stability of the power system. Although the numerical differentiation method has been widely used for approximately calculating the eigenvalue sensitivities, its accuracy has not been carefully investigated. Besides, the element-based formulation for computing closed-form eigenvalue sensitivities has not been used in any commercial software due to the average efficiency, complicated formulation, and error-prone characteristics. Based on the matrix calculus, this paper proposes an easily manipulated formulation of the closed-form eigenvalue sensitivities with respect to the power generation. The distinguishing feature of the formulation is that all the formulas consist of vector and matrix operations, which can be performed by developed numerical algorithms to take full advantages of architectural features of the modern computer. The tests on WSCC 3-machine 9-bus system, New England 10-machine 39-bus system, and IEEE 54-machine 118-bus system show that the accuracy of the proposed formulation is superior to the numerical differentiation method and the efficiency is also greatly improved compared to the element-based closed-form formulation. The proposed formulation will be helpful to perform a more accurate and faster stability analysis of a power grid with converter-based devices.

Keywords

Closed-form formulation; converter-based devices; eigenvalue sensitivity; matrix calculus; small-signal stability

I. Introduction

AS more and more power electronic converter-interfaced devices are integrated into the network on both the supply and demand sides, the power system is changing into the power-electronics-enabled system. Displacing synchronous machines with large inertia by converter-based generators with low inertia will result in a significant decrease in the system total inertia [

1]. Incidents of oscillation induced by the integration of converter-interfaced generators have been reported [2], [3]. These instability phenomena indicate that the small-signal stability is deteriorating in the power-electronics-enabled system.

Technically, eigenvalue analysis is a dominant method to small-signal stability problems of the power system. The participation factor formed by eigenvectors can help identify the dynamic variables that significantly affect a given mode or eigenvalue and has been widely used in commercial software for power system analysis [

4]. However, it cannot determine the increment or decrement of system parameters [5]. By contrast, the eigenvalue sensitivities with respect to system parameters can provide this critical information for both offline study and real-time control, which can help coordinate controller tuning or take remedial actions to suppress oscillations [6]-[12]. Besides, eigenvalue sensitivities can be treated as derivatives in each iteration when used in a mathematical optimization method [9], [10], [13]. The numerical differentiation method has been widely used for calculating eigenvalue sensitivity. It is easy to implement and can obtain an approximate solution [6]-[8], [10], [11], [14]. However, it is time-consuming due to the repetitive procedures. More importantly, it cannot obtain the eigenvalue sensitivities with respect to some operation parameters such as the power generation of the slack bus because the perturbation of them is invalid.

Alternatively, the closed-form eigenvalue sensitivity (CFES) only needs to be derived once and thus has higher efficiency. However, solving eigenvalue sensitivity needs to differentiate a complex implicit function in some cases, which makes the formulation of the closed-form very challenging. In [

15], an element-based closed-form formulation is proposed for eigenvalue sensitivities with respect to any arbitrary system parameter. The CFES has been applied to generation re-dispatch [9],[12], parameter optimization of power system stabilizer (PSS) [6], and stochastic analysis of grid-connected photovoltaic (PV) systems [16]. To the best of our knowledge, however, this formulation has not been used in any commercial software for eigenvalue analysis due to the average efficiency, complicated formulation, and error-prone characteristics. In addition, neither the special consideration for PV buses or the slack bus nor the evaluation of accuracy and efficiency between the numerical eigenvalue sensitivity (NES) and CFES has been discussed in the previous studies.

In this paper, an easily manipulated formulation of the CFES with respect to any power system parameter is proposed. The major contributions of the proposed formulation include: ① reformulating all the eigenvalue sensitivity formulas as combinations of only vector and matrix operations based on the matrix calculus, which can be performed by well-developed numerical algorithms to get high computation efficiency owing to reuse of data in the computer cache; ② transforming a three-dimensional matrix to a two-dimensional matrix by multiplying a constant vector preliminarily when formulating derivatives of state matrix with respect to variables, which can greatly improve the computation efficiency and make the formulation more intuitive. The proposed formulation also makes the derivation less error-prone owing to avoiding complex index in all formulas. Furthermore, a comprehensive comparison is performed between the CFES and the NES to show some interesting findings, such as the accuracy of NES declining as the size of the system increases.

All codes and detailed results of three benchmark systems and the test data have been made publicly available in GitHub [

17], which can be used to verify the correctness of the solutions.

The remainder of this paper is organized as follows. Section II introduces the small-signal stability model. Section III discusses the NES and the CFES and performs a theoretical comparison between them. Section IV describes the general mathematical formulation of the CFES based on matrix calculus. A detailed formulation through an example for calculating eigenvalue sensitivity with respect to the active power generation is presented. In Section V, the comparative analysis of CFES with the NES is performed on three benchmark systems. Finally, the conclusion is drawn in Section VI.

II. Small-signal Stability Model

A. State-space Analysis Model for Converter-fed Power Systems

The state-space approach is widely used for small-signal stability analysis for which the system is described by a set of differential algebraic equations (DAEs). The dynamic devices, which include synchronous machines and their regulators, converter-based devices such as flexible AC transmission system (FACTS) devices, high-voltage direct current (HVDC) devices and wind turbines, can be expressed as [

18]:

\dot{x} = F_{d} (x, y, u)

(1)

0 = G_{d} (x, y)

(2)

where $F_{d}$ is the vector of the differential equations; $G_{d}$ is the vector of the algebraic equations; $x$ is the vector of the device state variables; $y$ is the vector of device algebraic variables; and $u$ is the vector of the device input. For example, the internal current $I_{d}$ at d-axis, $I_{q}$ at q-axis and internal voltage $V_{d}$ at d-axis, $V_{q}$ at q-axis are often selected as algebraic variables for rotating machines such as synchronous machines and doubly-fed asynchronous wind turbines. Then they can be modeled as [

19]:

\dot{x} = F_{d} (x, I_{d}, I_{q}, V_{d}, V_{q}, u)

(3)

0 = G_{d} (x, I_{d}, I_{q}, V_{d}, V_{q})

(4)

As the model of individual converters or rotating machines is referring to its own rotating frame, a common reference frame is needed for all converters and rotating machines. An interface block is usually used to reflect the machine (converter)-network transformation. The transformation is defined as [

18]:

[\begin{matrix} F_{D} \\ F_{Q} \end{matrix}] = [\begin{matrix} s i n δ & c o s δ \\ - c o s δ & s i n δ \end{matrix}] [\begin{matrix} F_{d} \\ F_{q} \end{matrix}]

(5)

where $F$ may be either $I$ or $V$ . Thus

0 = V_{I} * s i n (δ - θ_{I}) - V_{d}

(6)

0 = V_{I} * c o s (δ - θ_{I}) - V_{q}

(7)

where $*$ denotes Hadamard product; $δ$ is the vector of the rotor angle; and $V_{I}$ and $θ_{I}$ are vectors of the voltage magnitude and phase angle for interface buses in common reference frame, respectively. Moreover, the power injection for the interface can be obtained through the following equations [

19]:

G_{p} = P_{I} - V_{d} * I_{d} - V_{q} * I_{q}

(8)

G_{q} = Q_{I} - V_{q} * I_{d} + V_{d} * I_{q}

(9)

where $P_{I}$ and $Q_{I}$ are the active and reactive power injections from the interface, respectively; $G_{p}$ is the vector of the active power injector equations; and $G_{q}$ is the vector of the reactive power injector equations.

Therefore, on a common reference frame, the network equations are given as:

P_{I} - P_{L} (V) - d i a g (V) (Y * c o s φ) V = 0

(10)

Q_{I} - Q_{L} (V) - d i a g (V) (Y * s i n φ) V = 0

(11)

where $P_{L} (V)$ and $Q_{L} (V)$ are the vectors of active and reactive load of all buses, respectively, which may relate to bus voltage amplitude $V$ ; $Y$ is the admittance matrix of the system; $d i a g (V)$ converts vector $V$ into a diagonal matrix; and $φ$ is the matrix of the composite phase angle. The element of $φ$ is $φ_{i j} = θ_{i} - θ_{j} - α_{i j}$ , where $α_{i j}$ is the admittance angle of branch $j$ connecting $i$ bus, and $θ_{i}$ and $θ_{j}$ are the bus voltage angles of bus $i$ and bus $j$ , respectively.

For small-signal stability analysis of conventional power systems, the axis transformation equations (6) and (7) for rotating machines can be substituted into (3) and (4) to obtain (12) and (13).

\dot{x} = F_{g} (x, I_{d}, I_{q}, V_{I}, θ_{I}, u)

(12)

0 = G_{g} (x, I_{d}, I_{q}, V_{I}, θ_{I})

(13)

where $F_{g}$ is the vector of the differential equations of the generator; and $G_{g}$ is the vector of the algebraic equations of the generator.

Substituting the axis transformation equations (6) and (7) into (8) and (9), (14) and (15) can be obtained.

P_{I} - V_{I} * s i n (δ - θ_{I}) * I_{d} - V_{I} * c o s (δ - θ_{I}) * I_{q} = 0

(14)

Q_{I} - V_{I} * c o s (δ - θ_{I}) * I_{d} + V_{I} * s i n (δ - θ_{I}) * I_{q} = 0

(15)

At last, combining the dynamic device equations (1) and (2), the power interface equations (14) and (15) for rotating machines and network equations (10) and (11), the state space model for conventional power system analysis can be described by the following DAEs:

\dot{x} = F (x, y, u)

(16)

0 = G (x, y)

(17)

where $F$ is the vector of the differential equations for all the devices; and $G$ is the vector of the algebraic equations for the devices, interface and network.

B. Linearization Technique and Initial Conditions

The key step in small-signal stability analysis is the linearization of DAEs. Linearizing (16) and (17), (18) can be obtained.

[\begin{matrix} Δ \dot{x} \\ 0 \end{matrix}] = [\begin{matrix} \tilde{A} & \tilde{B} \\ \tilde{C} & \tilde{D} \end{matrix}] [\begin{matrix} Δ x \\ Δ y \end{matrix}] + [\begin{matrix} E \\ 0 \end{matrix}] Δ u

(18)

where $\tilde{A}$ , $\tilde{B}$ , $\tilde{C}$ , and $\tilde{D}$ are matrices of all numbers with the values calculated under the initial conditions; $Δ x$ is the vector of the state variables after linearization; $Δ y$ is the vector of the algebraic variables after linearization; $Δ u$ is the vector of the input variables after linearization; and $E$ is the identity matrix.

Eliminating $Δ y$ , (19) can be obtained.

Δ \dot{x} = A Δ x + E Δ u

(19)

where A is commonly known as the state matrix and is given as:

A = \tilde{A} - \tilde{B} {\tilde{D}}^{- 1} \tilde{C}

(20)

Similar to (17), setting the $\dot{x}$ of (16) to be zero, (21) can be obtained.

0 = F (x, y, u)

(21)

Then the initial conditions of the variables for the small-signal stability model are calculated at an equilibrium point.

If all the real parts of the eigenvalues of $A$ are negative, the system is stable in small-signal stability sense according to Lyapunov theory. Usually, an index $η$ called spectral abscissa is introduced to describe the security margin:

η (A) = m a x \{R e (λ) : λ \in λ (A)\} = R e (λ_{η})

(22)

where $λ (A)$ represents all of the eigenvalues of $A$ ; $R e (λ)$ is the real part of an eigenvalue $λ$ ; and $λ_{η}$ is the eigenvalue with the largest real part.

III. NES and CFES

A. NES

According to the numerical differentiation method, eigenvalue analysis is performed to obtain the eigenvalue $λ (A)$ of the state matrix at the equilibrium point and then a system parameter P is varied by a small quantity $ε$ to get the perturbed state matrix $A_{ε}$ and its eigenvalue $λ (A_{ε})$ . The NES with respect to parameter P can be approximated by:

\frac{\partial λ}{\partial P} \approx \frac{λ (A_{ε}) - λ (A)}{ε}

(23)

Generally, the NES only cares about spectral abscissa sensitivity. It is the real part of the eigenvalue sensitivity, i.e.,

\frac{\partial η}{\partial P} = R e (\frac{\partial λ}{\partial P})

(24)

B. CFES

In fact, the CFES is a mathematical eigenvalue derived at an equilibrium point. Its formulation should base on the following formula in the mathematical theorem [

14], [20], [21]:

\frac{\partial λ_{i}}{\partial P} = \frac{ψ_{i}^{T} \frac{\partial A}{\partial P} ϕ_{i}}{ψ_{i}^{T} ϕ_{i}}

(25)

where $λ_{i}$ is the $i$ ^th eigenvalue; $ψ_{i}$ and $ϕ_{i}$ are the left and right eigenvectors of $λ_{i}$ at the equilibrium point, respectively, which can be calculated with the $λ_{i}$ .

Furthermore, if the state matrix $A$ is the explicit function of the parameter $P$ , the CFES with respect to the parameter can be calculated directly based on (25). The parameters of converters fall into this category, including the integral gain of the current control, time constant of the voltage control loop in wind turbines and photovoltaic cell regulator [

22]. Conversely, if the state matrix

A

is the implicit function of the parameter

P

, the formulation based on the chain rule will be a complicated process. Some operation parameters such as the active power output of the generators should be formulated according to this rule.

C. Comparison Between Numerical Spectral Abscissa Sensitivity and Closed-form Spectral Abscissa Sensitivity

Here the CFES with respect to a system parameter is compared with the numerical sensitivity. The spectral abscissa can be expressed with a function of the parameter vector $P$ as $η (P)$ , where $P$ consists of independent variables. According to the Taylor series expansion, the following relationship is obtained:

Δ η = \sum_{i \in S} \frac{\partial η}{\partial P_{i}} Δ P_{i} + \frac{\sum_{i \in S} \frac{\partial^{2} η}{\partial P_{i}^{2}}}{2!} Δ P_{i}^{2} + \dots + \frac{\sum_{i \in S} \frac{\partial^{n} η}{\partial P_{i}^{n}}}{n!} Δ P_{i}^{n} + \dots

(26)

where $S$ is the index set of the parameter vector. If the $i$ ^th parameter $P_{i}$ is varied by a small quantity while the other parameters remain unchanged, then (27) can be obtained.

Δ η = \frac{\partial η}{\partial P_{i}} Δ P_{i} + \frac{\frac{\partial^{2} η}{\partial P_{i}^{2}}}{2!} Δ P_{i}^{2} + \dots + \frac{\frac{\partial^{n} η}{\partial P_{i}^{n}}}{n!} Δ P_{i}^{n} + \dots

(27)

Dividing $Δ P_{i}$ on both sides of (27), (28) can be obtained.

\frac{Δ η}{Δ P_{i}} = \frac{\partial η}{\partial P_{i}} + \frac{\frac{\partial^{2} η}{\partial P_{i}^{2}}}{2!} Δ P_{i}^{2} + \dots + \frac{\frac{\partial^{n} η}{\partial P_{i}^{n}}}{n!} Δ P_{i}^{n - 1} + \dots

(28)

In fact, $Δ η / Δ P_{i}$ is the NES, while the mathematical eigenvalue derivative $\partial η / \partial P_{i}$ is the CFES in general. From (28), the NES which only considers the first-order Taylor series expansion in this case is just an approximation of the CFES.

Another case is that when one parameter is varied, it will cause changes in another two or more parameters. For instance, if the active power output of the $i$ ^th generator is varied with a small quantity, the active power output of the slack bus will also change in order to guarantee the power balance. Meanwhile, the reactive power output of the slack bus and all PV buses will also have slight variations. Thus, (29) can be obtained.

\frac{Δ η}{Δ P_{G i}} = \frac{\partial η}{\partial P_{G i}} + \frac{\partial η}{\partial P_{G s}} \frac{Δ P_{G s}}{Δ P_{G i}} + \frac{\partial η}{\partial Q_{G s}} \frac{Δ Q_{G s}}{Δ P_{G i}} + \sum_{j \in S_{P V}} \frac{\partial η}{\partial Q_{G j}} \frac{Δ Q_{G j}}{Δ P_{G i}} + \dots

(29)

where $Δ P_{G s}$ and $Δ Q_{G s}$ are the active and reactive power changes of the slack bus, respectively; and $S_{P V}$ is the set of PV buses. In this situation, the NES will have a bigger difference with the CFES.

Overall, the NES is calculated based on the numerical method, while the closed-form sensitivity is derived based on rigorous mathematical formulation. The NES is just an approximation of the closed-form sensitivity. However, the approximation may cause the optimality and convergence problem in optimization due to the inaccurate descent direction, which is shown in Section V. The accuracy is very important for controller parameter coordination and power generation redispatch.

IV. Mathematical Formulation of CFES Based on Matrix Calculus

With the advent of numerical programming, a set of well-developed algorithms for performing common linear algebra operations such as vector addition, scalar multiplication, Hadamard products, linear combinations, and matrix multiplication are becoming the de facto standard routines for linear algebra [

23], [24]. These algorithms leverage the idea of blocking to limit the amount of bus traffic in favor of high reuse of the data that is presented at a higher level and in faster memories. Therefore, they have extremely high computation efficiency. To take advantage of the architectural features, they force our engineers to reformulate our models to the new ones with only vector or matrix operations.

A. General Formulation of CFES

The proposed formulation of CFES with respect to parameter vector based on matrix calculus is first described generally. Depending on whether or not the state matrix $A$ is the explicit function of the parameter vector $P$ , two ways to calculate the CFES are presented.

1)　 $A$ is Explicit Function of $P$

According to (25), the eigenvalue sensitivity with respect to parametric vector $P$ is expressed as:

\frac{\partial λ_{i}}{\partial P} = \frac{ψ_{i}^{T} \frac{\partial A}{\partial P} ϕ_{i}}{ψ_{i}^{T} ϕ_{i}}

(30)

To solve $\partial λ_{i} / \partial P$ , the most important procedure is to formulate $\partial A / \partial P$ in (30). By substituting (20) into (30), (31) can be obtained.

\frac{\partial A}{\partial P} = \frac{\partial \tilde{A}}{\partial P} - \frac{\partial \tilde{B}}{\partial P} {\tilde{D}}^{- 1} \tilde{C} + \tilde{B} {\tilde{D}}^{- 1} \frac{\partial \tilde{D}}{\partial P} {\tilde{D}}^{- 1} \tilde{C} - \tilde{B} {\tilde{D}}^{- 1} \frac{\partial \tilde{C}}{\partial P}

(31)

Equation (31) converts the derivative of state matrix $A$ , whose elements are hard to represent, to that of some available block matrices, whose elements are easy to be derived from (18). In element-based formulation in [

15], the derivative for

A

with respect to the

j

^th parameter

P_{j}

is calculated by (32).

\begin{array}{l} \frac{\partial A}{\partial P_{j}} = \frac{\partial \tilde{A}}{\partial P_{j}} - \frac{\partial \tilde{B}}{\partial P_{j}} {\tilde{D}}^{- 1} \tilde{C} + \tilde{B} {\tilde{D}}^{- 1} \frac{\partial \tilde{D}}{\partial P_{j}} {\tilde{D}}^{- 1} \tilde{C} - \tilde{B} {\tilde{D}}^{- 1} \frac{\partial \tilde{C}}{\partial P_{j}} \\ j = 1,2, \dots, k \end{array}

(32)

where $k$ is the dimension of the parameter vector $P$ ; $\partial \tilde{A} / \partial P_{j}$ , $\partial \tilde{B} / \partial P_{j}$ , $\partial \tilde{C} / \partial P_{j}$ , and $\partial \tilde{D} / \partial P_{j}$ are all two-dimension matrices and the calculation of $\partial A / \partial P_{j}$ will need matrix operations of these two-dimensional matrices. Then the result of $\partial A / \partial P_{j}$ will be substituted into (25) to get $\partial λ_{i} / \partial P_{j}$ . To get $\partial λ_{i} / \partial P$ , the calculation of $\partial A / \partial P_{j}$ in (32) and $\partial λ_{i} / \partial P_{j}$ in (25) have to be executed for $k$ times. The element-based formulation can be implemented directly when it is programmed. But its performance is barely satisfactory as the formulation involves loops of matrix operations.

If $\partial λ_{i} / \partial P$ could be directly calculated instead of the loop calculation of $\partial λ_{i} / \partial P_{j}$ , the calculation efficiency will be increased sharply. However, $\partial A / \partial P$ is a three-dimension matrix, whose operations cannot be supported by the low-level routines. Since a row of left eigenvectors which are further obtained after the calculation of the $λ_{i}$ could be treated as a constant [

15] during the derivation, they are multiplied by the state matrix

A

to make

ψ_{i}^{T} A

become a row vector. Hence, as a two-dimensional matrix,

\partial (ψ_{i}^{T} A) / \partial P

is easy to express. Then (30) can be rewritten as:

\frac{\partial λ_{i}}{\partial P} = \frac{\frac{\partial (ψ_{i}^{T} A)}{\partial P} ϕ_{i}}{ψ_{i}^{T} ϕ_{i}}

(33)

Letting $τ^{T} = ψ_{i}^{T} \tilde{B} {\tilde{D}}^{- 1}$ which is also a constant vector and substituting it into (31), the formulation is obtained as:

\frac{\partial (ψ_{i}^{T} A)}{\partial P} = \frac{\partial (ψ_{i}^{T} \tilde{A})}{\partial P} - \frac{\partial (ψ_{i}^{T} \tilde{B})}{\partial P} {\tilde{D}}^{- 1} \tilde{C} + \frac{\partial (τ^{T} \tilde{D})}{\partial P} {\tilde{D}}^{- 1} \tilde{C} - \frac{\partial (τ^{T} \tilde{C})}{\partial P}

(34)

For practical computation, explicit evaluation of the inverse is avoided and the methods for solving sparse linear equations are used. For example, the calculation of $ψ_{i}^{T} \tilde{B} {\tilde{D}}^{- 1}$ just needs to solve (35) to get $τ$ .

{\tilde{D}}^{T} τ = {\tilde{B}}^{T} ψ_{i}

(35)

After the above formulation, the next step is to calculate $\partial (ψ_{i}^{T} \tilde{A}) / \partial P$ , $\partial (ψ_{i}^{T} \tilde{B}) / \partial P$ , $\partial (τ^{T} \tilde{C}) / \partial P$ , and $\partial (τ^{T} \tilde{D}) / \partial P$ , respectively. According to the parameter location in the structure of the matrices $\tilde{A}$ , $\tilde{B}$ , $\tilde{C}$ , and $\tilde{D}$ . If the structure of a matrix $\tilde{D}$ is :

\tilde{D} = [\begin{matrix} D_{11} & D_{12} & \dots & D_{1 n} \\ D_{21} & D_{22} & \dots & D_{2 n} \\ ⋮ & ⋮ & ⋮ \\ D_{m 1} & D_{m 2} & \dots & D_{m n} \end{matrix}]

(36)

And let $τ^{T} = [τ_{1}^{T}, τ_{2}^{T}, \dots, τ_{m}^{T}]$ , where the number of rows of $τ_{i}$ is equal to the number of rows of $D_{i 1}$ for $i = 1,2, \dots, m$ , then $\partial (τ^{T} \tilde{D}) / \partial P$ can further be expressed as:

\frac{\partial (τ^{T} \tilde{D})}{\partial P} = {[\begin{matrix} {(\frac{\partial (τ_{1}^{T} D_{11})}{\partial P} + \frac{\partial (τ_{2}^{T} D_{21})}{\partial P} + \dots + \frac{\partial (τ_{m}^{T} D_{m 1})}{\partial P})}^{T} \\ {(\frac{\partial (τ_{1}^{T} D_{12})}{\partial P} + \frac{\partial (τ_{2}^{T} D_{22})}{\partial P} + \dots + \frac{\partial (τ_{m}^{T} D_{m 2})}{\partial P})}^{T} \\ ⋮ \\ {(\frac{\partial (τ_{1}^{T} D_{1 n})}{\partial P} + \frac{\partial (τ_{2}^{T} D_{2 n})}{\partial P} + \dots + \frac{\partial (τ_{m}^{T} D_{m n})}{\partial P})}^{T} \end{matrix}]}^{T}

(37)

Through the above transformation, the computation efficiency will be greatly improved because all formulas are operations of two-dimensional matrices without any loop computation. Actually, this transformation has preliminarily performed some matrix multiplications during the derivation.

2)　 $A$ is Implicit Function of $P$

Select an appropriate variable vector $X$ of which the state matrix is the explicit function, and then write the $i$ ^th eigenvalue $λ_{i}$ as a function of parameter $X$ and $P$ as $λ_{i} (X (P), P)$ according to the derivation in Section II. The initial conditions of the variables for the small-signal stability model of power system in (17) and (21) are rewritten as:

0 = H (X (P), P)

(38)

where $H$ is the function vectors in the initial condition equations. Thus, $λ_{i}$ is an implicit function of $P$ .

To differentiate the function $λ_{i}$ , it is generally impossible to solve it explicitly with $X$ eliminated. Instead, the defined function $λ_{i}$ can be differentiated implicitly by the following implicit differentiation and the chain rule:

\frac{\partial λ_{i}}{\partial P} = {(\frac{\partial X}{\partial P})}^{T} \frac{\partial λ_{i}}{\partial X}

(39)

Here $\partial λ_{i} / \partial X$ can be calculated by the formulation proposed in Section VI. $\partial X / \partial P$ in (39) can be calculated by solving the following equations that are obtained by differentiating both sides of (38) with respect to $P$ :

0 = {(\frac{\partial X}{\partial P})}^{T} \frac{\partial H}{\partial X} + \frac{\partial H}{\partial P}

(40)

Therefore, the formulation for the implicit function of the parameter includes the one for the explicit function of the parameter. The detailed formulation will be discussed in Section IV through an example.

B. Formulation of $\partial λ_{i} / \partial P_{G}$

Since the eigenvalue sensitivity with respect to active power generation $P_{G}$ can provide useful information for remedial actions before or during the oscillation incident, the calculation of the eigenvalue sensitivity with respect to active power generation is important in the research of small-signal stability. The state matrix $A$ is the implicit function of active power generation $P_{G}$ . By the formulation of $\partial λ_{i} / \partial P_{G}$ for WSCC 3-machine 9-bus system which has a particular description in [

18], this subsection will present the details of the proposed formulation.

1)　Formulation of $\partial X / \partial P_{G}$

The appropriate variable vector $x$ can be obtained as follows:

x = {[\begin{matrix} δ^{T} & ω^{T} & {(E_{d}^{'})}^{T} & {(E_{q}^{'})}^{T} & E_{f d}^{T} & V_{R}^{T} & R_{F}^{T} \end{matrix}]}^{T}

(41)

where $ω$ is the rotor speed vector; $E_{d}^{'}$ and $E_{q}^{'}$ are the internal voltage vectors at d-axis and q-axis, respectively; $E_{f d}$ is the DC generator output voltage vector; $V_{R}$ is the voltage regulator output vector; and $R_{F}$ is the exciter rate feedback vector.

Let $X = {[x^{T}, I_{d}^{T}, I_{q}^{T}, V^{T}, θ^{T}]}^{T}$ , the differentiations of (10) with respect to $V$ and $θ$ are expressed as:

\frac{\partial P_{I}}{\partial V} - d i a g ((Y * c o s φ) V) - d i a g (V) (Y * c o s φ) = 0

(42)

\begin{array}{l} \frac{\partial P_{I}}{\partial θ} + d i a g (V) (d i a g ((Y * c o s φ) V) - \\ (Y * s i n φ) d i a g (V)) = 0 \end{array}

(43)

The differentiation of (11) can be derived similarly. Then $\partial V / \partial P_{G}$ and $\partial θ / \partial P_{G}$ are obtained by (42) and (43).

The differentiation of (14) with respect to $P_{G}$ is expressed as:

E + {(\frac{\partial V_{G}}{\partial P_{G}})}^{T} \frac{\partial G_{p}}{\partial V_{G}} + {(\frac{\partial θ_{G}}{\partial P_{G}})}^{T} \frac{\partial G_{p}}{\partial θ_{G}} + {(\frac{\partial δ}{\partial P_{G}})}^{T} \frac{\partial G_{p}}{\partial δ} + {(\frac{\partial I_{d}}{\partial P_{G}})}^{T} \frac{\partial G_{p}}{\partial I_{d}} + {(\frac{\partial I_{q}}{\partial P_{G}})}^{T} \frac{\partial G_{p}}{\partial I_{q}} = 0

(44)

where $E$ is an identity matrix; $V_{G}$ and $θ_{G}$ are the vectors of the voltage magnitude and phase angle for generator buses, respectively; and $\partial G_{p} / \partial V_{G}$ , $\partial G_{p} / \partial θ_{G}$ , $\partial G_{p} / \partial δ$ , $\partial G_{p} / \partial I_{d}$ , and $\partial G_{p} / \partial I_{q}$ are all diagonal matrices, which can be directly derived by differentiation rule as:

\frac{\partial G_{p}}{\partial V_{G}} = - d i a g (I_{d} * s i n (δ - θ_{G}) + I_{q} * c o s (δ - θ_{G}))

(45)

\frac{\partial G_{p}}{\partial θ_{G}} = d i a g (V_{G} * I_{d} * c o s (δ - θ_{G}) - V_{G} * I_{q} * s i n (δ - θ_{G}))

(46)

\frac{\partial G_{p}}{\partial δ} = - d i a g (V_{G} * I_{d} * c o s (δ - θ_{G}) - V_{G} * I_{q} * s i n (δ - θ_{G}))

(47)

\frac{\partial G_{p}}{\partial I_{d}} = d i a g (- V_{G} * s i n (δ - θ_{G}))

(48)

\frac{\partial G_{p}}{\partial I_{q}} = d i a g (- V_{G} * c o s (δ - θ_{G}))

(49)

The differentiation of (13) with respect to $P_{G}$ is expressed as:

E + {(\frac{\partial V_{G}}{\partial P_{G}})}^{T} \frac{\partial G_{g}}{\partial V_{G}} + {(\frac{\partial θ_{G}}{\partial P_{G}})}^{T} \frac{\partial G_{g}}{\partial θ_{G}} + {(\frac{\partial δ}{\partial P_{G}})}^{T} \frac{\partial G_{g}}{\partial δ} + {(\frac{\partial I_{d}}{\partial P_{G}})}^{T} \frac{\partial G_{g}}{\partial I_{d}} + {(\frac{\partial I_{q}}{\partial P_{G}})}^{T} \frac{\partial G_{g}}{\partial I_{q}} = 0

(50)

where $\partial G_{g} / \partial V_{G}$ , $\partial G_{g} / \partial θ_{G}$ , $\partial G_{g} / \partial δ$ , $\partial G_{g} / \partial I_{d}$ , and $\partial G_{g} / \partial I_{q}$ are all diagonal matrices, which can also be directly derived by differentiation rule. $\partial G_{g} / \partial V_{G}$ can be obtained as:

\frac{\partial G_{g}}{\partial V_{G}} = d i a g (- I_{d} * s i n (δ - θ_{G}) - I_{q} * c o s (δ - θ_{G}))

(51)

The differentiation of (15), (21) with respect to $P_{G}^{}$ can be derived similarly. Then a group of differentiated initial condition equations can be obtained. Solving these linear equations with $\partial V / \partial P_{G}$ and $\partial θ / \partial P_{G}$ from (42) and (43), $\partial x / \partial P_{G}$ , $\partial I_{d} / \partial P_{G}$ , and $\partial I_{q} / \partial P_{G}$ can be obtained.

2)　Formulation of $\partial λ_{i} / \partial X$

Since $\partial λ_{i} / \partial θ$ is the most complicated part among $\partial λ_{i} / \partial X$ , this paper only shows its formulation. The block diagonal matrices with respect to $θ$ is easy to derive the following differentiation rule. For example, the $D_{21}$ , which is a submatrix of $\tilde{D}$ , can be expressed as:

D_{21} = d i a g (- V_{G} * s i n (δ - θ_{G}))

(52)

The differentiation of $D_{21}$ with respect to $θ$ is:

\frac{\partial (τ_{2}^{T} D_{21})}{\partial θ} = d i a g (τ_{2} * V_{G} * c o s (δ - θ_{G}))

(53)

The differentiations of $D_{43}, D_{44}, D_{53},$ and $D_{54}$ are more complicated since they are full matrices. In fact, they are the Hessian matrices of the power flow equations of (10) and (11) as follows:

\frac{\partial (τ_{4}^{T} D_{43})}{\partial θ} = τ_{4}^{T} \frac{\partial^{2} G_{P}}{\partial V \partial θ} = (d i a g (b V) - a b^{T}) d + (- d i a g (τ_{4}^{T} a b) + a d b)

(54)

\frac{\partial (τ_{4}^{T} D_{44})}{\partial θ} = τ_{4}^{T} \frac{\partial^{2} G_{P}}{\partial θ \partial θ} = (d i a g (c V) - a c^{T}) a d - (- d i a g (τ_{4}^{T} a c) + a d c) a

(55)

\frac{\partial (τ_{5}^{T} D_{53})}{\partial θ} = τ_{5}^{T} \frac{\partial^{2} G_{Q}}{\partial V \partial θ} = (- d i a g (c V) + a c^{T}) e - (- d i a g (τ_{5}^{T} a c) + a e c)

(56)

\frac{\partial (τ_{5}^{T} D_{54})}{\partial θ} = τ_{5}^{T} \frac{\partial^{2} G_{Q}}{\partial θ \partial θ} = (d i a g (b V) - a b^{T}) a e - (- d i a g (τ_{5}^{T} a b) + a e b) a

(57)

where $a$ is the $d i a g (V)$ ; $b$ is the $Y * s i n φ$ ; $c$ is the $Y * c o s φ$ ; $d$ is the $d i a g (τ_{4}^{T})$ ; and $e$ is the $d i a g (τ_{5}^{T})$ .

3)　Calculation of CFES

Above all, (39) can be rewritten in the following form:

\begin{array}{l} \frac{\partial λ_{η}}{\partial P_{G}} = {(\frac{\partial X}{\partial P_{G}})}^{T} \frac{\partial λ_{η}}{\partial X} = {(\frac{\partial x}{\partial P_{G}})}^{T} \frac{\partial λ_{η}}{\partial x} + {(\frac{\partial I_{d}}{\partial P_{G}})}^{T} \frac{\partial λ_{η}}{\partial I_{d}} + \\ {(\frac{\partial I_{q}}{\partial P_{G}})}^{T} \frac{\partial λ_{η}}{\partial I_{q}} + {(\frac{\partial θ}{\partial P_{G}})}^{T} \frac{\partial λ_{η}}{\partial θ} + {(\frac{\partial V}{\partial P_{G}})}^{T} \frac{\partial λ_{η}}{\partial V} \end{array}

(58)

The real part of $\partial λ_{η} / \partial P_{G}$ is the closed-form spectral abscissa sensitivity (CFSAS).

In summary, the flow chart to calculate the CFES with respect to active power generation is shown in Fig. 1. The formulations proposed in this paper are all combinations of matrices and vectors. Therefore, by directly using some low-level routines, the computation efficiency will be improved. Furthermore, the proposed formulation is more clear and is not error-prone. It should also be noted that the proposed formulation can calculate all of the CFES, not just the spectral abscissa sensitivity.

Fig.1 Flow chart of computing CFES with respect to active power outputs.

4)　Additional Processing for PV and Slack Buses in CFSAS

At an equilibrium point, the voltage magnitude for PV and slack buses are fixed. Also, phase angles for slack buses are fixed. Actually, they are constants for the small-signal analysis model at an equilibrium point. Thus $\partial V_{k} / \partial P_{G}$ for a slack bus or PV bus $k$ in Section IV should be set to be zero. $\partial θ_{k} / \partial P_{G}$ in Section IV should be set to be zero if the bus $k$ is a slack bus.

However, since the operation point in each iteration is not basically an equilibrium point, there is no PV bus model in some power system optimization models with only one reference bus. To calculate the eigenvalue sensitivity with respect to $P_{G}$ , $\partial θ_{k} / \partial P_{G}$ in Section IV should be set to be zero if bus $k$ is a reference bus.

V. Case Studies

The proposed formulation and NES are both implemented in ANSI C language with CSparse of SuiteSparse [

24] and BLAS [23], which support matrix and vector operation. The code for element-based formulation uses ANSI C language. Since the cases are relatively small, LAPACK [25] is employed to calculate the eigenvalues and eigenvectors via QR decomposition.

A. Spectral Abscissa Sensitivities with Respect to Controller Parameters of Wind Turbines

A slightly modified version of the New England 10-machine 39-bus system [

26] is considered to assess the impact of the controller parameters of wind turbines on small-signal stability of the conventional power system. Two synchronous generators at buses 32 and 38 in the New England 10-machine 39-bus system are replaced by two doubly-fed induction generators. The model and the required data can be found in the PSAT MATLAB toolbox [27]. The perturbed quantity

ε

is set to be 0.1, which is small enough for calculating the NES in this case. The spectral abscissa sensitivity with respect to the parameters of pitch and voltage controller are presented in Table I, where NSAS stands for numerical spectral abscissa sensitivity. From Table I, both of NSAS and CFSAS are zero for pitch control gain

K_{p}

and pitch control time constant

T_{p}

. It implies that the two parameters of the pitch controller do not affect spectral abscissa. Besides, negative values of NSAS and CFSAS indicate that increasing the voltage control gain

K_{V}

will enhance the small-signal stability. This further demonstrates the same tendency and small difference of NSAS and CFSAS in the small system.

TABLE I Spectral Abscissa Sensitivity with Respect to Controller Parameters of Wind Turbines for New England 10-machine 39-bus System

Bus	$K_{p}$		$T_{p}$		$K_{V}$
Bus	NSAS	CFSAS	NSAS	CFSAS	NSAS	CFSAS
32	0	0	0	0	-0.0149	-0.0350
38	0	0	0	0	-0.0148	-0.0294

B. Spectral Abscissa Sensitivities with Respect to Power Generation

The proposed formulation is applied to three benchmark systems, namely the WSCC 3-machine 9-bus system [

18], New England 10-machine 39-bus system, and the modified IEEE 54-machine 118-bus system [28], and is compared with the numerical differentiation method. The generators are described by a two-axis model with an IEEE type DC-1 exciter. The loads are modeled as constant power. The required data for the three systems and detailed results in this paper can also be found in GitHub [17]. The perturbed quantity

ε

is set to be a enough small value 0.1 p.u. for calculating a numerically stable NES. For ease of comparison, only the spectral abscissa sensitivity is considered, which is most critical for power system small-signal stability. In these cases, the state matrix is the implicit function of power generation.

Detailed solutions are presented in Tables II-IV, where CFSAS-Y stands for the CFSAS with additional processing in Section IV and CFSAS-N stands for the CFSAS without additional processing.

1) Comparison of spectral abscissa sensitivities with respect to active power generation between NSAS, CFSAS-Y, and CFSAS-N

Note that the slack bus and the PV buses have no NES. The positive value in Tables II-IV means that if the power generation increases, the spectral abscissa sensitivity will become bigger and the system stability will be weakened. By contrast, the negative value in Tables II-IV means that if the power generation increases, the spectral abscissa sensitivity will become smaller and the system stability will be enhanced. It can be observed in Tables II-IV that the eigenvalue sensitivities of some generators obtained from NSAS, CFSAS-Y, and CFSAS-N can be significantly different, especially for IEEE 54-machine 118-bus system. Furthermore, some eigenvalue sensitivities have opposite signs such as those with respect to active power for the buses 25, 26, 59, 66, 80, 89, 100, and 103 of the IEEE 118-bus system, which indicates different adjustment directions.

TABLE II Spectral Abscissa Sensitivity for WSCC 3-machine 9-bus System

Bus	Active power			Reactive power
Bus	NSAS	CFSAS-Y	CFSAS-N	NSAS	CFSAS-Y	CFSAS-N
1		-0.0014	-0.2365		$- 0.0042$	-0.0673
2	0.1210	0.0794	0.0549		$- 0.0840$	0.0721
3	0.0369	-0.0114	-0.0364		$- 0.0187$	0.0671

To test the effectiveness of CFSAS, the following model (59) similar to the one in [

8] is used to simulate the redispatch of active power.

TABLE III Spectral Abscissa Sensitivity for New England 10-machine 39-bus System

Bus	Active power			Reactive power
Bus	NSAS	CFSAS-Y	CFSAS-N	NSAS	CFSAS-Y	CFSAS-N
30	$- 5.8 \times 10^{- 4}$	$- 1.4 \times 10^{- 4}$	$- 1.8 \times 10^{- 3}$		$- 2.1 \times 10^{- 5}$	$- 2.2 \times 10^{- 3}$
31		$- 9.3 \times 10^{- 6}$	$- 7.9 \times 10^{- 6}$		$- 5.9 \times 10^{- 6}$	$- 1.6 \times 10^{- 5}$
32	$1.0 \times 10^{- 5}$	$1.5 \times 10^{- 4}$	$- 1.9 \times 10^{- 4}$		$- 5.8 \times 10^{- 5}$	$- 3.2 \times 10^{- 5}$
33	$5.3 \times 10^{- 5}$	$8.1 \times 10^{- 5}$	$- 1.0 \times 10^{- 4}$		$- 1.7 \times 10^{- 4}$	$- 5.2 \times 10^{- 5}$
34	$- 1.8 \times 10^{- 5}$	$5.5 \times 10^{- 5}$	$7.2 \times 10^{- 5}$		$1.6 \times 10^{- 5}$	$7.3 \times 10^{- 6}$
35	$7.3 \times 10^{- 5}$	$1.9 \times 10^{- 4}$	$- 2.1 \times 10^{- 4}$		$- 9.6 \times 10^{- 5}$	$- 3.1 \times 10^{- 5}$
36	$- 1.1 \times 10^{- 4}$	$1.2 \times 10^{- 5}$	$4.5 \times 10^{- 5}$		$4.3 \times 10^{- 6}$	$2.5 \times 10^{- 5}$
37	$8.2 \times 10^{- 3}$	$6.3 \times 10^{- 3}$	$- 6.9 \times 10^{- 3}$		$- 3.1 \times 10^{- 2}$	$- 3.6 \times 10^{- 2}$
38	$- 6.3 \times 10^{- 6}$	$1.4 \times 10^{- 4}$	$- 1.1 \times 10^{- 4}$		$- 1.1 \times 10^{- 6}$	$- 8.6 \times 10^{- 5}$
39	$- 7.3 \times 10^{- 5}$	$3.2 \times 10^{- 4}$	$1.8 \times 10^{- 4}$		$- 9.1 \times 10^{- 8}$	$7.5 \times 10^{- 5}$

TABLE IV Spectral Abscissa Sensitivity for IEEE 54-machine 118-bus System

Bus	Active power			Reactive power
Bus	NSAS	CFSAS-Y	CFSAS-N	NSAS	CFSAS-Y	CFSAS-N
10	$3.1 \times 10^{- 6}$	$2.1 \times 10^{- 10}$	$1.5 \times 10^{- 10}$		$- 4.0 \times 10^{- 11}$	$- 1.0 \times 10^{- 9}$
12	$3.5 \times 10^{- 6}$	$1.1 \times 10^{- 10}$	$1.3 \times 10^{- 10}$		$- 1.1 \times 10^{- 11}$	$1.2 \times 10^{- 10}$
25	$2.6 \times 10^{- 6}$	$- 5.0 \times 10^{- 9}$	$- 1.1 \times 10^{- 8}$		$- 3.9 \times 10^{- 11}$	$- 1.9 \times 10^{- 9}$
26	$2.7 \times 10^{- 6}$	$- 6.9 \times 10^{- 10}$	$8.8 \times 10^{- 9}$		$- 1.3 \times 10^{- 11}$	$1.8 \times 10^{- 9}$
49	$1.3 \times 10^{- 5}$	$2.5 \times 10^{- 3}$	$- 6.3 \times 10^{- 2}$		$1.6 \times 10^{- 5}$	$- 2.3 \times 10^{- 2}$
54	$- 3.5 \times 10^{- 2}$	$- 3.3 \times 10^{- 2}$	$- 2.0 \times 10^{- 2}$		$5.9 \times 10^{- 3}$	$1.2 \times 10^{- 2}$
59	$- 8.3 \times 10^{- 5}$	$1.9 \times 10^{- 3}$	$- 8.2 \times 10^{- 2}$		$1.5 \times 10^{- 6}$	$- 2.6 \times 10^{- 2}$
61	$- 4.1 \times 10^{- 5}$	$- 8.3 \times 10^{- 5}$	$3.7 \times 10^{- 3}$		$- 3.5 \times 10^{- 11}$	$1.1 \times 10^{- 3}$
65	$- 1.1 \times 10^{- 5}$	$- 3.2 \times 10^{- 7}$	$- 7.5 \times 10^{- 5}$		$- 1.8 \times 10^{- 11}$	$- 3.8 \times 10^{- 7}$
66	$4.5 \times 10^{- 6}$	$- 4.8 \times 10^{- 5}$	$2.3 \times 10^{- 3}$		$3.8 \times 10^{- 9}$	$7.9 \times 10^{- 4}$
69		$4.7 \times 10^{- 11}$	$9.7 \times 10^{- 3}$		$- 1.6 \times 10^{- 10}$	$3.9 \times 10^{- 3}$
80	$- 4.2 \times 10^{- 6}$	$2.5 \times 10^{- 8}$	$2.0 \times 10^{- 8}$		$- 1.7 \times 10^{- 10}$	$6.9 \times 10^{- 8}$
89	$- 3.6 \times 10^{- 6}$	$2.4 \times 10^{- 9}$	$3.9 \times 10^{- 8}$		$- 1.8 \times 10^{- 11}$	$- 1.8 \times 10^{- 7}$
100	$- 4.0 \times 10^{- 6}$	$2.0 \times 10^{- 9}$	$7.5 \times 10^{- 10}$		$- 4.4 \times 10^{- 11}$	$- 5.0 \times 10^{- 9}$
103	$- 4.2 \times 10^{- 6}$	$7.6 \times 10^{- 10}$	$5.3 \times 10^{- 10}$		$- 2.1 \times 10^{- 11}$	$- 1.8 \times 10^{- 10}$
111	$- 4.3 \times 10^{- 6}$	$- 5.7 \times 10^{- 10}$	$- 3.0 \times 10^{- 10}$		$6.0 \times 10^{- 12}$	$2.2 \times 10^{- 9}$

\{\begin{array}{l} m i n \sum_{i \in S_{G}} Δ P_{G i}^{2} \\ s . t . \sum_{i \in S_{G}} Δ P_{G i}^{2} = 0 \\ \sum_{i \in S_{G}} σ_{P i} Δ P_{G i} \leq \bar{η} - η_{0} \\ {\underset{̲}{P}}_{G i} \leq P_{G i}^{0} + Δ P_{G i} \leq {\bar{P}}_{G i} i \in S_{G} \end{array}

(59)

where $σ_{P i}$ is the spectral abscissa sensitivity of the active power output of the $i$ ^th generator at the equilibrium point; $η_{0}$ and $P_{G i}^{0}$ are the spectral abscissa and active power output of the equilibrium point which are obtained from a conventional optimal power flow (OPF) with the object function minimizing the generation cost, respectively; $Δ P_{G i}$ is the change of the output power of the $i$ ^th generator from the equilibrium point; $\bar{η}$ is the expected spectral abscissa sensitivity; and ${\underset{̲}{P}}_{G i}$ and ${\bar{P}}_{G i}$ are the minimum and maximum power outputs of the $i$ ^th generator, respectively.

Tests on the New England 10-machine 39-bus system are performed with $η_{0}$ as $- 0.10$ and the expected spectral abscissa sensitivity as $- 0.12$ and $- 0.13$ . In each case, either CFSAS-Y or NSAS is used to provide descent directions. The results are summarized in Table V. It can be observed that when using CFSAS-Y, a smaller amount of active power output adjustment can obtain the same spectral abscissa, indicating that CFSAS-Y can provide better descent direction than NSAS for this optimization problem.

TABLE V Active Power Generation Adjustment for 10-machine 39-bus System

Bus	$P_{G}^{0}$ (MW)	$Δ P_{G}$ (MW)
		$\bar{η} = - 0.12$		$\bar{η} = - 0.13$
		NSAS	CFSAS-Y	NSAS	CFSAS-Y
30	0	262.33	107.51	352.69	122.50
31	754.57	-7.07	-7.07	-7.07	-7.07
32	920.00	0.00	0.00	0.00	0.00
33	0.00	0.00	62.20	0.00	70.82
34	747.50	-31.11	0.00	-16.80	0.00
35	862.50	-263.10	-129.83	-308.93	-148.24
36	0.00	0.00	357.86	0.00	408.09
37	805.00	-313.29	-373.08	-372.13	-425.74
38	883.26	151.74	58.80	151.74	66.94
39	1179.50	200.50	-76.39	200.50	-87.28
Total		1229.13	1172.73	1409.87	1336.67

CFSAS-N fails in all cases, which suggests that additional processing for PV and slack buses is essential to the redispatch problem in (59). This is because the model (59) needs a spectral abscissa sensitivity at an equilibrium point. However, it needs to be emphasized that CFSAS with only slack bus processing is very suitable to be used in the small-signal stability constrained optimal power flow (SSSC-OPF) model, as shown in [

9]. In the OPF model, CFSAS-Y does not exist because there are no PV buses. Besides, since the operation point in each iteration is basically not an equilibrium point, NSAS cannot be solved either. Therefore, CFSAS with slack bus processing is the only method that can be used to calculate the spectral abscissa sensitivity in the SSSC-OPF model. The results in [9] show that the method using CFSAS with slack bus processing has good convergence and optimality.

2) Comparison of spectral abscissa sensitivities with respect to reactive power generation between NSAS, CFSAS-Y, and CFSAS-N

As Hopf bifurcations are associated with eigenvalue conditions, the reactive power also plays an important role in small-signal stability [

29]. Besides, the cost for adjusting reactive power is lower and the response time for adjusting reactive power is shorter because the exciter has a smaller time constant than the governor. Thus, the reactive power re-dispatch is also important to improve small-signal stability. Because most of the generators in the test systems are set as PV buses or slack bus, NSAS is invalid for reactive power redispatch. The following modified model (60) is used to simulate the re-dispatch with CFSAS-Y and CFSAS-N in Section V.

\{\begin{array}{l} m i n \sum_{i \in S_{G}} Δ P_{G i}^{2} \\ s . t . \sum_{i \in S_{G}} Δ P_{G i} = 0 \\ \sum_{i \in S_{G}} Δ Q_{G i} = 0 \\ \sum_{i \in S_{G}} σ_{P i} Δ P_{G i} + \sum_{i \in S_{G}} σ_{Q i} Δ Q_{G i} \leq \bar{η} - η_{0} \\ {\underset{̲}{P}}_{G i} \leq P_{G i}^{0} + Δ P_{G i} \leq {\bar{P}}_{G i} i \in S_{G} \\ {\underset{̲}{Q}}_{G i} \leq Q_{G i}^{0} + Δ Q_{G i} \leq {\bar{Q}}_{G i} i \in S_{G} \end{array}

(60)

where $σ_{Q i}$ is the eigenvalue sensitivity of the reactive power output of the $i$ ^th generator; $Δ Q_{G i}$ is the change of the reactive power output of the $i$ ^th generator from the base case; and ${\underset{̲}{Q}}_{G i}$ and ${\bar{Q}}_{G i}$ are the minimum and maximum reactive power outputs of the $i$ ^th generator, respectively.

The changes of the active power and reactive power outputs for all the generators are shown in Table VI. It can be observed that considering the eigenvalue sensitivity with respect to the reactive power outputs of generator, the expected spectral abscissa sensitivity can be achieved without the need to re-dispatch the active power outputs. It further implies that the reactive power outputs of generator should not be ignored in small-signal stability analysis. Besides, CFSAS-N needs more adjustment of reactive power than CFSAS-Y. This shows that CFSAS-Y should also be used in the reactive power re-dispatching according to spectral abscissa sensitivity.

TABLE VI Active Power and Reactive Power Generation Adjustment for 10-machine 39-bus System

Bus	$Q_{G}^{0}$ (Mvar)	$Δ P_{G}$ (MW)		$Δ Q_{G}$ (Mvar)
Bus	$Q_{G}^{0}$ (Mvar)	CFSAS-Y	CFSAS-N	CFSAS-Y	CFSAS-N
30	182.67	1.41	1.41	-147.43	-116.19
31	271.94	-7.07	-7.07	-42.49	-47.46
32	294.12	0	0	9.01	-46.89
33	103.99	1.41	1.41	-32.72	-39.85
34	188.75	0	0	-40.61	-39.58
35	225.62	0	0	-103.52	-206.65
36	64.93	1.41	1.41	-52.83	147.82
37	36.49	0	0	461.89	460.32
38	45.04	1.41	1.41	-4.19	-17.07
39	98.66	1.41	1.41	-47.11	-94.44
Total		14.14	14.14	941.79	1216.27

C. Efficiency

The time efficiency of the proposed closed-form formulation is compared with the numerical differentiation method and element-based formulation. All simulations are performed on an HP EliteOne with Intel Core i5-6500 3.20 GHz CPU and 8 GB of RAM memory without GPU hardware. From Table VII, it can be observed that the proposed CFES calculation is more time-efficient than the element-based formulation and the numerical differentiation method, especially for large systems. The calculation time of the numerical differentiation method is about 12 times as much as the proposed formulation for IEEE 118-bus system.

TABLE VII Time of Calculating CFSAS-Y and NSAS

Test system	Calculation time (s)
Test system	Proposed formulation	Element-based formulation	Numerical differentiation
WSCC 9-bus system	0.008	0.008	0.010
New England 39-bus system	0.016	0.080	0.067
IEEE 118-bus system	0.189	2.250	2.355

VI. Conclusion

The accuracy of the numerical eigenvalue decreases as the size of the system increases compared with the numerical differentiation method. The proposed CFES has higher accuracy and can provide better descent direction for the optimization problem used for improving the system small-signal stability. Besides, the efficiency of the proposed formulation of CFES is substantially higher than that of the element-based formulation and the numerical differentiation method. Also, the proposed formulation will be helpful for coordinating controller tuning and taking remedial actions.

References

J. Quintero, V. Vittal, G. T. Heydt et al., “The impact of increased penetration of converter control-based generators on power system modes of oscillation,” IEEE Transactions on Power Systems, vol. 29, no. 5, pp. 2248-2256, Sept. 2014. [Baidu Scholar]

P. Belkin, “Event of 10-22-09,” in Proceedings of CREZ Technical Conference, Taylor, USA, Jan. 2010. [Baidu Scholar]

R. S. M. Wu and L. Xie, “Parameter sensitivity analysis for subsybchronous control interactions in wind-integrated power systems,” in Proceedings of CIGRE Grid of the Future Symposium, Houston, USA, Oct. 2014. [Baidu Scholar]

X. Chen, M. Shi, H. Sun et al., “Distributed cooperative control and stability analysis of multiple DC electric springs in a DC microgrid,” IEEE Transactions on Industrial Electronics, vol. 65, no. 7, pp. 5611-5622, Jul. 2018. [Baidu Scholar]

N. J. B. P. Kundur and M. G. Lauby, Power System Stability and Control. New York: McGraw-hill, 1994. [Baidu Scholar]

R. Jalayer and B. Ooi, “Co-ordinated pss tuning of large power systems by combining transfer function-eigenfunction analysis (TFEA), optimization, and eigenvalue sensitivity,” IEEE Transactions on Power Systems, vol. 29, no. 6, pp. 2672-2680, Nov. 2014. [Baidu Scholar]

A. T. Sarić and A. M. Stanković, “Rapid small-signal stability assessment and enhancement following changes in topology,” IEEE Transactions on Power Systems, vol. 30, no. 3, pp. 1155-1163, May 2015. [Baidu Scholar]

C. Y. Chung, L. Wang, F. Howell et al., “Generation rescheduling methods to improve power transfer capability constrained by smallsignal stability,” IEEE Transactions on Power Systems, vol. 19, no. 1, pp. 524-530, Feb. 2004. [Baidu Scholar]

P. Li, J. Qi, J. Wang et al., “An SQP method combined with gradient sampling for small-signal stability constrained OPF,”IEEE Transactions on Power Systems, vol. 32, no. 3, pp. 2372-2381, May 2017. [Baidu Scholar]

R. Zarate-Minano, F. Milano, and A. J. Conejo, “An OPF methodology to ensure small-signal stability,” IEEE Transactions on Power Systems, vol. 26, no. 3, pp. 1050-1061, Aug. 2011. [Baidu Scholar]

J. Cao, W. Du, H. Wang et al., “A novel emergency damping control to suppress power system inter-area oscillations,” IEEE Transactions on Power Systems, vol. 28, no. 3, pp. 3165-3173, Aug. 2013. [Baidu Scholar]

C. Li, H. Chiang, and Z. Du, “Network-preserving sensitivity-based generation rescheduling for suppressing power system oscillations,”IEEE Transactions on Power Systems, vol. 32, no. 5, pp. 3824-3832, Sept. 2017. [Baidu Scholar]

J. Condren and T. W. Gedra, “Expected-security-cost optimal power flow with small-signal stability constraints,” IEEE Transactions on Power Systems, vol. 21, no. 4, pp. 1736-1743, Nov. 2006. [Baidu Scholar]

H. K. Nam, Y. K. Kim, K. S. Shim et al., “A new eigen-sensitivity theory of augmented matrix and its applications to power system stability analysis,” IEEE Transactions on Power Systems, vol. 15, no. 1, pp. 363-369, Feb. 2000. [Baidu Scholar]

K. W. Wang, C. Y. Chung, C. T. Tse et al., “Multimachine eigenvalue sensitivities of power system parameters,” IEEE Transactions on Power Systems, vol. 15, no. 2, pp. 741-747, May 2000. [Baidu Scholar]

S. Liu, P. X. Liu, and X. Wang, “Stochastic small-signal stability analysis of grid-connected photovoltaic systems,” IEEE Transactions on Industrial Electronics, vol. 63, no. 2, pp. 1027-1038, Feb. 2016. [Baidu Scholar]

GitHub. Eigenvalue sensitivity in GitHub. (2018, Dec.).[Online]. Available: https://github.com/eigenlpj/Eigenvalue-Sensitivity [Baidu Scholar]

P. W. Sauer and M. A. Pai, Power System Dynamics and Stability. Bergen County: Prentice Hall, USA, 1998. [Baidu Scholar]

F. Milano, Power System Modelling and Scripting. Boston: Springer, 2010. [Baidu Scholar]

J. R. Magnus, “On differentiating eigenvalues and eigenvectors,” Econometric Theory, vol. 1, no. 2, pp. 179-191, Aug. 1985. [Baidu Scholar]

F. Alvarado and C. Demarco, “Avoiding and suppressing oscillations,” [Baidu Scholar]

Power Systems Engineering Research Center, Cornell University, Ithaca, New York, Dec., 1999. [Baidu Scholar]

Y. Li, Z. Xu, and K. P. Wong,“Advanced control strategies of PMSG-based wind turbines for system inertia support,” IEEE Transactions on Power Systems, vol. 32, no. 4, pp. 3027-3037, Oct. 2017. [Baidu Scholar]

C. L. Lawson, R. J. Hanson, D. R. Kincaid et al., “Basic linear algebra subprograms for fortran usage,” ACM Transactions on Mathematical Software, vol. 5, no. 3, pp. 308-323, Sept. 1979. [Baidu Scholar]

T. A. Davis, Fundamentals of Algorithms Direct Methods for Sparse Linear Systems. Philadelphia: SIAM, 2006. [Baidu Scholar]

C. B. E. Anderson and Z. Bai, Lapack Users’ Guide, Release 2.0. Philadelphia: SIAM, 1995. [Baidu Scholar]

R. H. Yeu, “Small signal analysis of power systems: eigenvalue tracking method and eigenvalue estimation contingency screening for DSA,” Ph.D. dissertation, Electrical and Computer Engineering College, University of Illinois at Urbana-Champaign, Urbana, 2010. [Baidu Scholar]

F. Milano. Power system analysis toolbox (PSAT). (2018, Apr.). [Online]. Available: http://faraday1.ucd.ie/software.html [Baidu Scholar]

KIOS Research Center. IEEE 118-bus modified test system. (2018, Apr.). [Online]. Available: http://www.kios.ucy.ac.cy/testsystems/index.php [Baidu Scholar]

A. I. Zecevic and D. M. Miljkovic, “The effects of generation redispatch on Hopf bifurcations in electric power systems,” IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, vol. 49, no. 8, pp. 1180-1186, Aug. 2002. [Baidu Scholar]

Address:No.19 Chengxin Avenue, Jiangning District, Nanjing 211106, China

E-mail: mpce@alljournals.cn

Tel:86-25-81093060

Fax:86-25-81093040

Home

Introduction

Editorial Board

For Author

Call For Papers

APC

Sponsor & Publisher

A Closed-form Formulation of Eigenvalue Sensitivity Based on Matrix Calculus for Small-signal Stability Analysis in Power System PDF

Abstract

Keywords

I. Introduction

II. Small-signal Stability Model

A. State-space Analysis Model for Converter-fed Power Systems

B. Linearization Technique and Initial Conditions

III. NES and CFES

A. NES

B. CFES

C. Comparison Between Numerical Spectral Abscissa Sensitivity and Closed-form Spectral Abscissa Sensitivity

IV. Mathematical Formulation of CFES Based on Matrix Calculus

A. General Formulation of CFES

B. Formulation of $\partial λ_{i} / \partial P_{G}$

V. Case Studies

A. Spectral Abscissa Sensitivities with Respect to Controller Parameters of Wind Turbines

B. Spectral Abscissa Sensitivities with Respect to Power Generation

C. Efficiency

VI. Conclusion

References

Home

Introduction

Editorial Board

For Author

Call For Papers

APC

Sponsor & Publisher

A Closed-form Formulation of Eigenvalue Sensitivity Based on Matrix Calculus for Small-signal Stability Analysis in Power System PDF

Abstract

Keywords

I. Introduction

II. Small-signal Stability Model

A. State-space Analysis Model for Converter-fed Power Systems

B. Linearization Technique and Initial Conditions

III. NES and CFES

A. NES

B. CFES

C. Comparison Between Numerical Spectral Abscissa Sensitivity and Closed-form Spectral Abscissa Sensitivity

IV. Mathematical Formulation of CFES Based on Matrix Calculus

A. General Formulation of CFES

B. Formulation of ∂λi/∂PG

V. Case Studies

A. Spectral Abscissa Sensitivities with Respect to Controller Parameters of Wind Turbines

B. Spectral Abscissa Sensitivities with Respect to Power Generation

C. Efficiency

VI. Conclusion

References

B. Formulation of $\partial λ_{i} / \partial P_{G}$