Small-signal Stability Constrained Optimal Power Flow Considering Eigenvalue Distribution

Zheng Huang; Kewen Wang; Yi Wang; Jing Han; Jun Liang

网刊加载中。。。

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

确定继续浏览么?

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

Small-signal Stability Constrained Optimal Power Flow Considering Eigenvalue Distribution PDF

- ORCID：
Zheng Huang ¹
✉
- ORCID：
Kewen Wang ¹
✉
- ORCID：
Yi Wang ¹
✉
- ORCID：
Jing Han ¹
✉
- ORCID：
Jun Liang ²
✉

1. School of Electrical and Information Engineering, Zhengzhou University, Zhengzhou 450001, China； 2. Cardiff University, Cardiff CF24 3AA, U.K

Updated：2024-07-26

DOI：10.35833/MPCE.2023.000135

OUTLINE

Abstract

In the existing small-signal stability constrained optimal power flow (SSSC-OPF) algorithms, only the rightmost eigenvalue or eigenvalues that do not satisfy a given threshold, e.g., damping ratio threshold and real-part threshold of eigenvalue, are considered in the small-signal stability constraints. The effect of steady-state, i.e., operating point, changes on eigenvalues is not fully taken into account. In this paper, the small-signal stability constraint that can fully reflect the eigenvalue change and system dynamic performance requirement is formed by analyzing the eigenvalue distribution on the complex plane. The small-signal stability constraint is embedded into the standard optimal power flow model for generation rescheduling. The simultaneous solution formula of the SSSC-OPF is established and solved by the quasi-Newton approach, while penalty factors corresponding to the eigenvalue constraints are determined by the stabilization degree of constrained eigenvalues. To improve the computation speed, a hybrid algorithm for eigenvalue computation in the optimization process is proposed, which includes variable selection for eigenvalue estimation and strategy selection for eigenvalue computation. The effectiveness of the proposed algorithm is tested and validated on the New England 10-machine 39-bus system and a modified practical 68-machine 2395-bus system.

Keywords

Damping ratio; eigenvalue; optimal power flow; quasi-Newton approach; sensitivity; small-signal stability

A. Indices and Sets

$Ω_{C}^{(l)}$	Set of eigenvalues to be considered in the $l^{t h}$ unconstrained optimization
$Ω_{D, k}$	Set of dominant variable serial numbers for electromechanical oscillation mode k
$Ω_{D}^{(l)}$	Set of dominant variable serial numbers in the $l^{t h}$ unconstrained optimization
$Ω_{S}$	Set of serial numbers of all state variables
$j 1$ , $j 2$	Iteration numbers of one-dimensional search in the $l^{t h}$ unconstrained optimization, and $j 1 < j 2$
l	Iteration number of quasi-Newton approach
B. Parameters
$α, β$	Penalty factors for real-part constraints of eigenvalue and damping ratio constraints
$ε$	Convergence tolerance
$κ$	Scaling factor
$ϖ$	Threshold for determining N_R in Algorithm 1
$σ_{C}$	Real-part threshold of eigenvalue, and $σ_{C} < 0$
$ζ_{T}, ζ_{C}$	Damping ratio thresholds, and $ζ_{C} > ζ_{T}$
$Δ t$	Time interval
$a_{i}, b_{i},$ $c_{i}$	Cost coefficients of the $i^{t h}$ generator
$A$ , $I$	State matrix and identity matrix
$I_{m a x}$	The maximum iteration number in Algorithm 1
$m, r$	Numbers of power flow equations and technical limits
$N_{M}$	Number of considered contingency scenarios
$N_{G}, N_{S}$	Numbers of generators and state variables
$R_{G}$	Vector of ramping limits
S	Threshold for selecting eigenvalue computation strategy
$u_{i}, v_{j}$	Penalty factors for the $i^{t h}$ technical limit and the $j^{t h}$ small-signal stability constraint
C. Variables
$σ_{k}, ω_{k}$	Real and imaginary parts of the k^th eigenvalue
$σ_{r}$	The largest real part of eigenvalues
$λ_{k}, ζ_{k}$	The k^th eigenvalue and the corresponding damping ratio
$λ_{m}$	Eigenvalue corresponding to the minimum damping ratio
$ζ_{m}$	The minimum damping ratio of all modes
$Δ ζ_{m}$	Variation of the minimum damping ratio between optimal power flow (OPF) and small-signal stability constrained optimal power flow (SSSC-OPF) solutions
$Δ σ_{r}$	Variation of the largest real part of eigenvalues between OPF and SSSC-OPF solutions
$Δ L$	Variation of augmented Lagrangian function value between OPF and SSSC-OPF solutions
$η^{(l)}$	Optimal step length for the $l^{t h}$ unconstrained optimization
$d^{(l)}$ , $s^{(l)}$	Deviations of parameter and gradient vectors associated with the i^th unconstrained optimization
$H$	Approximation to inverse of Hessian matrix of augmented Lagrangian function
$H^{(l)}$	H at operating point K⁽^l⁾
$k_{r}$ , $k_{t}$	Total numbers of H restarts and optimal step length search
$k_{Q R}$ , $k_{I R A}$	Total number of calls to QR algorithm and implicitly restarted Arnoldi approach
$k_{E S}$	Total number of calls to the first-order eigenvalue sensitivity computation subroutine without sensitivity computation for updating H
K	Column vector consisting of state variables and Lagrangian multipliers
L(K⁽^l⁾)	Augmented Lagrangian function value at K⁽^l⁾
N_λ	Number of modes considered in $Ω_{C}^{(l)}$
N_R	Number of selected dominant state variables
$P_{G i}$	Active power output of the $i^{t h}$ generator
R, Q	Intermediate variables in Algorithm 1
U_k, W_k	Right and left eigenvectors corresponding to the k^th eigenvalue
$V_{G i}$	Voltage magnitude of the $i^{t h}$ generator
$w_{i}^{(\cdot)}$	Sensitivity of $ζ_{k}$ with respect to state variable x_i in the $l^{t h}$ unconstrained optimization
x, y	Column vectors of state variables and Lagrangian multipliers

I. Introduction

POWER flow distribution affects both the steady-state and stability performance of power systems. For computational reason, the standard optimal power flow (OPF) model usually does not consider stability constraints or simply takes the form of a branch phase angle difference [

1], [2]. With the significant improvement in computing power, it is feasible to consider stability constraints in the OPF model [3]-[5]. Small-signal stability constrained optimal power flow (SSSC-OPF) has gained considerable popularity in power system operating as economic efficiency and small-signal stability are both considered. However, solving the SSSC-OPF problem is very time-consuming compared with the OPF. There is always a problem of how to reduce the impact of the small-signal stability constraint (SSSC) on the computation speed [6].

The SSSC in the SSSC-OPF model needs not only to adequately describe the stability of the system, but also to consider the computation burden. Therefore, for small-scale power systems, the SSSC includes all or some of the eigenvalues; for large-scale power systems, only several critical eigenvalues can be considered due to the computation burden. In [

7] and [8], the SSSCs are formulated in terms of the sensitivity of the minimum damping ratio and the sensitivity of the maximum real part of the eigenvalue, respectively. In [9]-[11], the SSSC is that the largest real part of the eigenvalues is less than a threshold. An important issue, as reported in [7]-[11], is that the minimum damping ratio or the largest real part of the eigenvalues can be chosen as the small-signal stability index. However, as the system operating point changes, the eigenvalue corresponding to the stability index may change from one to another, which may lead to oscillations among several critical modes in the solution process [7], [10]. In [12]-[14], the constrained eigenvalues in the SSSC include the real part of all eigenvalues. To improve the overall dynamic performance of the system, the SSSC, which is described by all eigenvalues with damping ratio no less than the given requirement, is incorporated into the OPF model [15]. However, since the effect of steady-state changes on eigenvalues is not fully considered in the SSSC, eigenvalue oscillations may still occur during the solution process. To alleviate the oscillation of eigenvalues, an active set strategy is utilized to update the constrained eigenvalues [16]. However, it should be pointed out that in this approach, the minimum damping ratio threshold is changed successively, so that the optimal operating point is approached step by step. Therefore, in this paper, the SSSC is constructed considering the effect of steady-state changes on eigenvalues, and the eigenvalue oscillation in the simultaneous solution of the SSSC-OPF problem is also addressed.

At present, the SSSC-OPF problem is mainly solved by gradient-based optimization algorithms [

7]-[13], [15]-[19] and artificial intelligence algorithms such as differential evolution algorithm [14], self-adaptive evolutionary programming [20], and back propagation neural network algorithm [21]. Most of these algorithms are based on eigenvalue analysis, but it is time-consuming to compute eigenpairs repeatedly in the solution process. To reduce the computation requirement, partial eigenvalue approaches are utilized to compute the eigenvalues in the SSSC-OPF problem, e.g., Rayleigh’s iteration [13] and implicitly restarted Arnoldi (IRA) approach [7], [12], [15]. Alternatively, critical eigenvalues can be estimated by using the eigenvalue sensitivity approach, which avoids the formation of the state matrix and the computation of eigenpairs. However, it may result in the small-signal stability index being smaller or larger than the given threshold. To deal with this problem, a scaling factor is introduced to correct for eigenvalue sensitivity in eigenvalue estimation [7], [13]. This problem can also be solved by a hybrid algorithm using approximate (eigenvalue sensitivity approach) and accurate computations, e.g., QR algorithm [22] or partial eigenvalue approaches.

In this paper, a practical algorithm for the simultaneous solution of the SSSC-OPF problem is proposed. The SSSC is formed here by analyzing the eigenvalue distribution on the complex plane, which is composed of the real part of eigenvalue and the damping ratio of some critical modes. These critical modes include all unstable and underdamped modes as well as those that may become unstable and underdamped as the system operating point changes. The SSSC is integrated with the standard OPF model, and the penalty factors for these constrained eigenvalues are determined by their corresponding stabilization degree. The quasi-Newton approach is employed to simultaneously solve the proposed SSSC-OPF model, and a hybrid algorithm using approximate and accurate computation is proposed to compute eigenvalues in the optimization process.

The remainder of this paper is organized as follows. Section II presents the SSSC. The model and algorithm of the SSSC-OPF are given in Section III. In Section IV, the effectiveness of the proposed algorithm is tested and validated on the New England 10-machine 39-bus system and a modified practical 68-machine 2395-bus system. Finally, Section V outlines the conclusions.

II. SSSC

For small-signal stability, the state matrix A is formed by the linearized model of power systems at an equilibrium point [

22]. The stability of a system is determined by the eigenvalues of A. Let

λ_{k}

be a complex eigenvalue, we can obtain:

λ_{k} = σ_{k} \pm j ω_{k}

(1)

Each pair of the complex eigenvalues corresponds to an oscillatory mode. The damping ratio $ζ_{k}$ of an oscillatory mode is given by:

ζ_{k} = - \frac{σ_{k}}{\sqrt[]{σ_{k}^{2} + ω_{k}^{2}}}

(2)

For a power system with N_G generators, there are usually $N_{G} - 1$ electromechanical oscillation modes. To ensure the system dynamic performance, the eigenvalues of all electromechanical oscillation modes must have negative real parts, and the damping ratio should not be less than a given threshold $ζ_{T}$ , e.g., 5% [

16]. Thus, the traditional SSSC1 can be expressed as:

\{\begin{matrix} σ_{k} < 0 \\ ζ_{k} \geq ζ_{T} \end{matrix} k = 1,2, \dots, N_{G} - 1

(3)

For SSSC-OPF problems, if the SSSC1 of (3) is not satisfied, the system operating point will be adjusted to improve these modes. However, as the system operating point changes, some stable modes may become unstable and vice versa, i.e., eigenvalue oscillations may occur during the solution process. This oscillation may lead to poor convergence of the algorithm, or even make the algorithm difficult to converge.

To alleviate the oscillation of constrained eigenvalues during the solution process, the complex plane is divided into three parts in this paper, as shown in Fig. 1.

Fig. 1 Distribution range of eigenvalues.

Region 1 indicates the modes that do not satisfy SSSC1. All eigenvalues in region 1 are included in $Ω_{C}$ . Region 2 represents the modes with insufficient stability or security margin, i.e., the eigenvalues are close to the imaginary axis or $ζ_{T}$ . Since the eigenvalues move during the optimization process, the eigenvalues in region 2 may cross the boundary to region 1. Hence, the eigenvalues in region 2 should be considered in $Ω_{C}$ . The modes in region 3 are relatively stable and can be ignored. However, if the eigenvalues in region 3 move to region 2, they should also be considered in $Ω_{C}$ . From the above analysis, it can be observed that $Ω_{C}$ not only contains multiple critical modes, but also changes dynamically. Mathematically, $Ω_{C}$ can be described by:

Ω_{C} = {λ_{k} : σ_{k} \geq σ_{C}} ⋃ {λ_{k} : ζ_{k} < ζ_{C}} k = 1,2, \dots, N_{G} - 1

(4)

And the proposed SSSC2 can be expressed as:

\{\begin{matrix} σ_{k} < σ_{C} \\ ζ_{k} \geq ζ_{C} \end{matrix} k = 1,2, \dots, N_{G} - 1

(5)

If $σ_{C} = 0$ and $ζ_{C} = ζ_{T}$ , the SSSC2 of (5) becomes the SSSC1. The SSSC2 not only reflects the dynamic performance requirement of the system, but also fully considers the effect of the operating point changes on eigenvalues. It should be noted that the purpose of SSSC is not to require all modes to satisfy the SSSC2, but to prevent the eigenvalues in region 2 from moving to region 1 when the modes in region 1 are improved during the optimization process. To achieve this purpose, different penalty factors are set for eigenvalues in different regions.

III. Model and Algorithm of SSSC-OPF

A. Model of SSSC-OPF

SSSC-OPF problems are proposed by adding the SSSC to the standard OPF model [

12], [16]. The SSSC-OPF formulation can be written as:

m i n f (x)

(6)

s . t .

h (x) = 0

(7)

g (x) \leq 0

(8)

q (x) \leq 0

(9)

where $f (x)$ indicates the objective function; $h (x) = [h_{1} (x),$ $h_{2} (x), \dots, h_{m} {(x)]}^{T}$ is the vector corresponding to the equality constraints, e.g., power flow equations; $g (x) = [g_{1} (x), g_{2} (x), \dots,$ $g_{r} {(x)]}^{T}$ is the vector corresponding to the inequality constraints, including technical limits for active power, reactive power, voltage magnitudes, and branch flow; and $q (x) = [q_{1} (x), q_{2} (x), \dots, q_{2 N_{G} - 2} {(x)]}^{T}$ is the vector corresponding to the compact form of the SSSC, which is the implicit function of x.

B. Problem Reformulation

In general, the SSSC-OPF problem is a nonlinear programming problem in mathematics. There are many approaches to deal with this problem [

23], such as gradient-based optimization algorithms and artificial intelligence algorithms. In this paper, the Lagrangian multiplier approach [24] is utilized to deal with the equality constraint. Because of its simplicity and ease of implementation, the exterior penalty function approach [1] is employed to deal with the inequality constraints. Then, an augmented Lagrangian function can be formed as:

\begin{array}{l} L (x, y) = f (x) + y^{T} h (x) + \sum_{i = 1}^{r} (u_{i} \cdot m a x (0, g_{i} {(x)))}^{2} + \\ \sum_{j = 1}^{2 N_{G} - 2} (v_{j} \cdot m a x (0, q_{j} {(x)))}^{2} \end{array}

(10)

q_{2 k - 1} (x) = σ_{k} - σ_{C}

(11)

q_{2 k} (x) = ζ_{C} - ζ_{k}

(12)

Penalty factors for SSSCs are given as:

v_{2 k - 1} = \{\begin{array}{l} α & σ_{k} \geq 0 \\ κ α & σ_{C} \leq σ_{k} < 0 \\ 0 & σ_{k} < σ_{C} \end{array}

(13)

v_{2 k} = \{\begin{array}{l} β ζ_{k} < ζ_{T} \\ κ β ζ_{T} \leq ζ_{k} < ζ_{C} \\ 0 ζ_{k} \geq ζ_{C} \end{array}

(14)

If a mode is located in region 1 ( $σ_{k} \geq 0$ or $ζ_{k} < ζ_{T}$ ), the system is underdamped or unstable under small disturbances. To improve this mode, $α$ and $β$ need to be set relatively large. However, the penalty factors for the modes in region 2 ( $σ_{C} \leq σ_{k} < 0$ or $ζ_{T} \leq ζ_{k} < ζ_{C}$ ) should be relatively small. In other words, $κ$ is less than 1. This is because the penalty factors for the modes in region 2 are not to improve the modes, but to prevent the eigenvalues in region 2 from moving to region 1. It is effective to select the penalty factors by starting with a low value and then increasing it during the optimization process [

1].

C. Optimization Technique

The Karush-Kuhn-Tucker (KKT) conditions for the SSSC-OPF problem of (10) can be given as follows. $\nabla_{x} L$ and $\nabla_{y} L$ are gradient vectors of the Lagrangian with respect to state variables and Lagrangian multipliers, respectively.

\nabla_{x} L (x, y) = 0

(15)

\nabla_{y} L (x, y) = 0

(16)

Since the second-order eigenvalue sensitivity is very time-consuming [

25], the quasi-Newton approach is utilized to solve the KKT conditions (15) and (16). Let

K = {[x y]}^{T}

, an iterative formula is given by:

K^{(l + 1)} = K^{(l)} - η^{(l)} H^{(l)} \nabla L (K^{(l)})

(17)

where $\nabla L (K^{(l)}) = [\nabla_{x} L (K^{(l)})^{T} \nabla_{y} L (K^{(l)})^{T}]^{T}$ ; and $η^{(l)}$ can be obtained by using the one-dimensional search [

26]. In this paper, the Broyden-Fletcher-Goldfarb-Shanno (BFGS) approach [26] is used for updating matrix H.

H^{(l + 1)} = (I - \frac{s^{(l)} (d^{(l)})^{T}}{(d^{(l)})^{T} s^{(l)}}) H^{(l)} (I - \frac{d^{(l)} (s^{(l)})^{T}}{(d^{(l)})^{T} s^{(l)}}) + \frac{s^{(l)} (s^{(l)})^{T}}{(d^{(l)})^{T} s^{(l)}}

(18)

\{\begin{array}{l} s^{(l)} = - η^{(l)} H^{(l)} \nabla L (K^{(l)}) \\ d^{(l)} = \nabla L (K^{(l + 1)}) - \nabla L (K^{(l)}) \end{array}

(19)

After each iteration, the inequality constraints (8) and (9) are checked and $L (K^{(l)})$ is changed. Near the optimal solution of $L (K^{(l)})$ , the variation of $K^{(l)}$ is relatively small. However, due to the highly nonlinear relationship between eigenvalues and operating parameters [

25], any change in operating parameters may result in significant changes of eigenvalues. Although the range of eigenvalues considered is extended by the SSSC2, the modes considered in

L (K^{(l)})

may also be changed. If the modes considered in

L (K^{(l + 1)})

are different from those in

L (K^{(l)})

d^{(l)}

used to update the matrix

H^{(l + 1}

⁾ may be ineffective and may even provide a wrong direction. Therefore, in this case,

H^{(l + 1)}

has to be restarted from the identity matrix.

D. A Hybrid Algorithm for Eigenvalue Calculation

During the simultaneous solution of the SSSC-OPF using the quasi-Newton approach, $η^{(l)}$ is obtained by the one-dimensional search in the $l^{t h}$ unconstrained optimization. In the process of searching for $η^{(l)}$ , it is necessary to repeatedly compare the value of the augmented Lagrangian function L(K). This requires repeatedly forming the state matrix and calculating eigenvalues, which is time-consuming. For small-scale power systems, an accurate eigenvalue computation, i.e., QR algorithm, IRA approach, can be performed for each iteration. For large-scale power systems, the computation burden of eigenvalues is much greater than that of solving the equilibrium point. To reduce the computation cost, the eigenvalues considered in $L (K^{(l)})$ can also be computed by a hybrid algorithm using accurate eigenvalue computation and eigenvalue estimation.

1) Variable Selection for Eigenvalue Estimation

Since the number of state variables in large-scale power systems may be large, it is also time-consuming to compute the sensitivities of multiple eigenvalues with respect to all state variables in x. To reduce the sensitivity computation time, it is necessary to reduce the number of variables in x for eigenvalue estimation. This means that only dominant variables are used for eigenvalue estimation.

The small-signal stability index $ζ$ is utilized to explain how to select the dominant variables. For an electromechanical oscillation mode k, the first-order Taylor series expansion of $ζ_{k}^{(j 2)}$ at an operating point $x^{(j 1)}$ is given in (20), which can be rewritten as (21).

ζ_{k}^{(j 2)} = ζ_{k}^{(j 1)} + \sum_{i = 1}^{N_{S}} \frac{\partial ζ_{k}^{(j 1)}}{\partial x_{i}^{(j 1)}} (x_{i}^{(j 2)} - x_{i}^{(j 1)})

(20)

Δ ζ_{k} = ζ_{k}^{(j 2)} - ζ_{k}^{(j 1)} = \sum_{i = 1}^{N_{S}} w_{i}^{(j 1)} Δ x_{i}

(21)

In [

17], for an electromechanical oscillation mode, state variables with relatively large damping ratio sensitivity are selected for eigenvalue estimation. However, it can be observed from (21) that

Δ ζ_{k}

is related to both

w_{i}^{(j 1)}

and

Δ x_{i}

. Since

Δ x_{i}

is determined by (17), there is no guarantee that if

Δ ζ_{k}

is relatively large,

Δ x_{i}

will also be relatively large. If

w_{i}^{(j 1)}

is relatively small but

Δ x_{i}

is relatively large,

w_{i}^{(j 1)} Δ x_{i}

may also be relatively large. Compared with

w_{i}^{(j 1)}

w_{i}^{(j 1)} Δ x_{i}

reflects the impact on

Δ ζ_{k}

more accurately. Therefore, in this subsection, variable selection according to

w_{i} Δ x_{i}

is proposed for eigenvalue estimation of multiple electromechanical oscillation modes.

To reduce the computation time, if $w_{i}^{(j 1)} Δ x_{i}$ in (21) has a relatively small effect on $Δ ζ_{k}$ , it can be ignored. However, all $w_{i}^{(j 1)}$ are unknown beforehand, because the purpose of variable selection is to select $w_{i}^{(j 1)}$ that needs to be computed. Fortunately, all $w_{i}^{(l)}$ have been computed in each unconstrained optimization, and $w_{i}^{(l)} Δ x_{i}$ can be used for variable selection. For a mode k, the absolute values of all $w_{i}^{(l)} Δ x_{i}$ are sorted in decreasing order, i.e.,

|w_{1}^{(l)} Δ x_{1}| > |w_{2}^{(l)} Δ x_{2}| > \dots > |w_{N_{R}}^{(l)} Δ x_{N_{R}}| > \dots > |w_{N_{S}}^{(l)} Δ x_{N_{S}}|

(22)

where $|\cdot|$ represents the absolute value.

Assuming that the first N_R $|w_{i}^{(l)} Δ x_{i}|$ in (22) have a relatively large effect on $Δ ζ_{k}$ , and the first N_R $|w_{i}^{(l)} Δ x_{i}|$ corresponding to the state variables are considered as the dominant variables and will be selected. The value of N_R is determined by Algorithm 1.

Algorithm 1 : determine value of N_R

1: Initialization: set $i = 1$ , $I_{m a x} = N_{S}$ , $N_{R} = 0$ , $R = 0$ , and $Q = 0$ ; choose a threshold $ϖ$

2: while $i < I_{m a x}$ do

3: set $R \leftarrow R + w_{i}^{(l)} Δ x_{i}$ , $Q \leftarrow |w_{i + 1}^{(l)} Δ x_{i + 1} / R|$ , $N_{R} \leftarrow i$

4: if $Q < ϖ$

5: Output N_R and stop

6: end if

7: set $i \leftarrow i + 1$

8: end while

For the mode k, $Ω_{D, k}$ is given by:

Ω_{D, k} = {1,2, \dots, N_{R}} k \in Ω_{C}^{(l)}

(23)

Multiple modes are included in $Ω_{C}^{(l)}$ , and dominant variables for different modes are usually inconsistent. Therefore, for multiple modes, the dominant variables are selected as:

Ω_{D}^{(l)} = Ω_{D, 1} ⋃ Ω_{D, 2} ⋃ \dots ⋃ Ω_{D, k} k \in Ω_{C}^{(l)}

(24)

Note that the computation complexity of the proposed variable selection approach is very low. To reduce the impact of variable selection on the optimization process, dominant variables can be reselected when $w_{i}^{(j 1)}$ needs to be computed in each iteration.

2) Strategy Selection for Eigenvalue Computation

The hybrid algorithm can be used for eigenvalue computation in the one-dimensional search process, but it is still a problem to choose an appropriate eigenvalue computation strategy (accurate eigenvalue computation or eigenvalue estimation) in each iteration. An accurate eigenvalue computation can be done after several iterations, and the eigenvalue estimation will be used for the rest of the iterations. However, this may lead to poor convergence or even non-convergence of the algorithm because the variation of state variables may be large. Therefore, the following condition, as given in (25), is utilized to choose an appropriate eigenvalue computation strategy.

\sum_{i = 1}^{N_{S}} |x_{i}^{(j 1)} - x_{i}^{(j 2)}| \leq S

(25)

Assume that an accurate eigenvalue computation is performed at $x^{(j 1)}$ and the value of S is given. If (25) is not satisfied at operating point x⁽^j²⁾, the eigenvalues corresponding to $x^{(j 2)}$ must be computed accurately. If (25) is satisfied at $x^{(j 2)}$ , the eigenvalue sensitivity approach will be utilized to compute the eigenvalues until (25) is not satisfied at $x^{(j 2)}$ .

3) Eigenvalue Estimation with Sensitivity Approach

If (25) is satisfied at $x^{(j 2)}$ , it is assumed that the modes considered in $L (K^{(j 2)})$ are the same as those considered in $L (K^{(j 1)})$ . Based on the eigenvalue sensitivity approach, the eigenvalues at x⁽^j²⁾ can be estimated by:

λ_{k} (x^{(j 2)}) = λ_{k} (x^{(j 1)}) + \sum_{i} (x_{i}^{(j 2)} - x_{i}^{(j 1)}) \frac{\partial λ_{k} (x^{(j 1)})}{\partial x_{i}^{(j 1)}}

(26)

\frac{\partial λ_{k}}{\partial x_{i}} = W_{k}^{T} \frac{\partial A}{\partial x_{i}} U_{k} k \in Ω_{C}^{(j 1)}

(27)

\{\begin{array}{l} i \in Ω_{S} x^{(j 1)} = x^{(l)} \\ i \in Ω_{D}^{(l)} x^{(j 1)} \neq x^{(l)} \end{array}

(28)

The derivative of A with respect to state variable x_i can be obtained using plug-in modeling technology [

25].

E. Flowchart of Proposed Algorithm

For the sake of simplicity, a computation flowchart of the proposed algorithm for the SSSC-OPF problem is given in Fig. 2.

Fig. 2 Flowchart of proposed algorithm.

F. Extension of Proposed Algorithm

The extension of the proposed algorithm considering both the normal and contingency operating conditions will be discussed below. The OPF model considering SSSCs under the normal and contingency operating conditions of the power system can be formulated as follows. The subscripts 0 and $𝓁$ represent the normal operating condition and contingency scenario, respectively.

m i n f (x_{0})

(29)

s.t.

h_{0} (x_{0}) = 0

(30)

G_{0} (x_{0}) \leq 0

(31)

h_{𝓁} (x_{𝓁}) = 0 𝓁 = 1,2, \dots, N_{M}

(32)

G_{𝓁} (x_{𝓁}) \leq 0 𝓁 = 1,2, \dots, N_{M}

(33)

|x_{0} - x_{𝓁}| \leq R_{G} Δ t 𝓁 = 1,2, \dots, N_{M}

(34)

$G_{0} (x_{0}) \leq 0$ is the compact form of constraints (8) and (9). Formula (34) represents linking constraints relating state variables of the normal operating condition and state variables of any of the contingency scenarios.

The problem in (29)-(34) can be solved by several algorithms, e.g., Benders decomposition utilized in this paper. Using the Benders decomposition, the problem can be decomposed into a normal operating master problem and a set of contingency sub-problems. The whole problem is iteratively solved between the master problem and sub-problems. The proposed algorithm in this paper can be utilized to solve both the master problem and sub-problems. A detailed description of the Benders decomposition used in optimization problems can be found in [

27] and [28].

IV. Case Studies

In this section, the proposed algorithm is applied to the New England 10-machine 39-bus system and a modified practical 68-machine 2395-bus system to illustrate the effectiveness. For all systems, the loads are modeled as constant impedance. Except for the slack generator, P_G_i and V_G_i of all generators are included in the state variables. The objective function is to minimize the cost of power generation, which can be expressed as:

m i n f (x) = \sum_{i = 1}^{N_{G}} (a_{i} P_{G i}^{2} + b_{i} P_{G i} + c_{i})

(35)

Double precision computation is typically used in small-signal stability analysis, because solving eigenpairs is computationally intensive and introduces relatively large round-off errors. The scaling factor $κ$ is set to be 10^-2. The initial values of the penalty factors $α$ and $β$ are set to be 100 and the maximum value is set to be 10000. After each unconstrained optimization, the penalty factors corresponding to be the violated constraints will be increased proportionally [

1], [29]. The convergence tolerance

ε

is set to be 10^-6.

A. New England 10-machine 39-bus System

The New England 10-machine 39-bus system is often used for stability analysis. Details of network parameters, nodal power, and dynamic parameters can be found in [

30]. The generator cost coefficients can be found in [2]. The voltage magnitude limits are from [10]. The generator limits are from [13]. The generators are described by the fourth-order model. All the generators are equipped with an IEEE type-1 exciter, except for the generator 1 which represents an equivalent of the New York network. The number of eigenvalues is 65. The parameters of the SSSC are chosen as

σ_{C} = - 0.05

ζ_{T} = 3 %

, and

ζ_{C} = 3.5 %

1) Effectiveness of Proposed Algorithm

To analyze the effectiveness of the proposed algorithm, the following four cases are taken into account. All cases are solved by the proposed algorithm, where all eigenvalues are computed by the QR algorithm.

1) Case 1: OPF without SSSC.

2) Case 2: OPF with SSSC1 (OPF-SSSC1), without considering H restart.

3) Case 3: OPF-SSSC1, considering H restart.

4) Case 4: OPF with SSSC2 (OPF-SSSC2), considering H restart.

Simulation results and electromechanical oscillation modes in four cases for New England 10-machine 39-bus are shown in Table I and Fig. 3, respectively.

TABLE I Simulation Results in Four Cases for New England 10-machine 39-bus System

Case	L	ζ_m(%)	l	$k_{t}$	$k_{r}$
Case 1	40308.18	1.07	9	106
Case 2	43310.49	2.96	9	120
Case 3	43210.34	2.99	14	187	4
Case 4	42836.08	3.01	9	117	1

Fig. 3 Electromechanical oscillation modes in four cases for New England 10-machine 39-bus system.

From Table I and Fig. 3, it can be observed that the minimum damping ratio $ζ_{m}$ of Case 1 and Case 2 are lower than the desired damping ratio $ζ_{T}$ . The minimum damping ratio of Case 3 is acceptable, and all electromechanical oscillation modes obtained by Case 4 reach $ζ_{T}$ . Comparing the results of Case 2 and Case 3, it can be observed that although the iteration number of Case 2 is less than that of Case 3, the augmented Lagrangian function value of Case 2 is higher than that of Case 3. This is because the eigenvalues considered in the optimization process have changed and H has not been restarted from the identity matrix, resulting in optimization results far from the optimal point. The simulation results of Case 2 and Case 3 show that H restart is necessary for solving the SSSC-OPF problem by the quasi-Newton approach. Comparing the results of Case 3 and Case 4, it can be observed that l and $k_{t}$ of Case 3 are more than those of Case 4. This is because the number of H restarts $k_{r}$ in Case 3 is four times that in Case 4. H restart means that the algorithm degenerates from the quasi-Newton approach to the steepest descent approach, which affects the convergence of the algorithm.

The eigenvalue changes for Case 3 and Case 4 are provided by Table II. For Case 3, it can be observed from Table II that eigenvalue oscillations between two modes occur in the first 9 iterations, after which the minimum damping ratio varies above and below the given damping ratio threshold, and the constrained mode remains until convergence after 11 iterations. Although the modes considered in Case 4 also change, the damping ratios of all modes are above the given threshold after 6 iterations, and the constrained modes remain until convergence after 3 iterations. The simulation results show that the proposed algorithm can efficiently alleviate the eigenvalue oscillation and improve the convergence.

TABLE II Eigenvalue Changes for Case 3 and Case 4

Case 3			Case 4
λ_m	ζ_m (%)	$N_{λ}$	λ_m	ζ_m (%)	$N_{λ}$
$- 0.0667 \pm j 6.4528$	1.03	1	$- 0.0667 \pm j 6.4528$	1.03	2
$- 0.0667 \pm j 6.4528$	1.03	1	$- 0.1681 \pm j 6.5620$	2.56	2
$- 0.1475 \pm j 6.5588$	2.24	1	$- 0.1591 \pm j 6.3751$	2.49	3
$- 0.1475 \pm j 6.5588$	2.24	1	$- 0.1960 \pm j 6.5526$	2.99	3
$- 0.1475 \pm j 6.5588$	2.24	1	$- 0.1960 \pm j 6.5526$	2.99	3
$- 0.1807 \pm j 6.5451$	2.76	2	$- 0.1902 \pm j 6.3388$	3.00	3
$- 0.1933 \pm j 6.4959$	2.97	2	$- 0.1902 \pm j 6.3388$	3.00	3
$- 0.1884 \pm j 6.3118$	2.98	2	$- 0.1984 \pm j 6.5806$	3.01	3
$- 0.1905 \pm j 6.3954$	2.97	1	$- 0.1984 \pm j 6.5806$	3.01	3
$- 0.1928 \pm j 6.3957$	3.01	0
$- 0.1911 \pm j 6.3942$	2.99	1
$- 0.1911 \pm j 6.3942$	2.99	1
$- 0.1911 \pm j 6.3942$	2.99	1
$- 0.1911 \pm j 6.3942$	2.99	1

2) Feasibility of Hybrid Algorithm Without Variable Selection

Since the New England 10-machine 39-bus system is relatively small, there is no need to reduce the number of state variables in x. Taking Case 4 as an example, the SSSC-OPF results under different threshold S are provided in Table III. From Table III, it can be observed that minimum damping ratios all reach $ζ_{T}$ , and the variation of the augmented Lagriangan function value is within an acceptable range. Also, when the threshold S is 0.01, 0.05, 0.15, or 0.25, l is more than that of Case 4, but the computation time required is less than that of Case 4. This is because $k_{Q R}$ is dramatically reduced. When the variation of state variables is less than $S$ , eigenvalues are computed by the eigenvalue sensitivity approach. Note that the computation time may not decrease as the threshold S increases because eigenvalues are the nonlinear functions of state variables. Furthermore, the value of S must be within a reasonable range according to the simulations. If S takes a relatively small value, the improvement in computation speed will not be significant; if S takes a relatively large value, it may cause a significant increase in iterations or even non-convergence of the algorithm.

TABLE III Performance Comparison of Proposed Algorithm Under Different S for New England 10-machine 39-bus System

S	L	ζ_m (%)	l	$k_{t}$	$k_{r}$	$k_{Q R}$	$k_{E S}$	Computation time (s)
	42836.08	3.01	9	117	1	117	0	44.73
0.01	42590.68	3.01	12	170	1	67	4	31.14
0.05	42720.13	3.01	12	143	1	30	4	22.12
0.15	42780.83	3.00	12	153	1	24	4	21.13
0.25	42826.50	3.01	12	165	1	16	3	18.96

3) Comparison with Existing Approaches

In [

10], the SSSC is described by the largest real part of the eigenvalues. Due to the eigenvalue oscillations, a sequential quadratic programming approach with gradient sampling (SQP-GS) is introduced to solve the SSSC-OPF problem. In [16], the minimum damping ratio is selected as the small-signal stability index, and a sequential approach (SA) is proposed to alleviate the eigenvalue oscillations. The minimum damping ratio threshold is changed successively, so that the optimal operating point is approached step by step. Therefore, the proposed algorithm is compared with SQP-GS and SA. The oscillation of constrained eigenvalues has been considered in different ways by SQP-GS, SA, and the proposed algorithm. SA is an alternating solution algorithm, while SQP-GS and the proposed algorithm are simultaneous solution ones, which have advantages in convergence and computation speed. Furthermore, to improve the computation speed, a hybrid algorithm for eigenvalue computation in the optimization process is proposed in this paper, which includes variable selection for eigenvalue estimation and strategy selection for eigenvalue computation. Therefore, the proposed algorithm has an advantage in terms of computation speed.

The results of different algorithms for New England 10-machine 39-bus system are provided in Table IV. Computer configurations and runtime environments are listed in Table V. In Table IV, the results of this paper are obtained by using the proposed algorithm without variable selection under $S = 0.01$ , the results of the SQP-GS are from [

10], and the results of the SA are from [16]. Both SQP-GS and SA have been tested in the New England 10-machine 39-bus system, but the system parameters used, i.e., steady-state parameters and dynamic parameters, are from different references. Considering the differences in system parameters and computing resources, the results presented in Table IV are for reference only. As shown in Table IV, all algorithms can improve the small-signal stability of the system while ensuring its economy. The difference in the economic loss

Δ L

required by different algorithms may be due to the inconsistency of steady-state parameters and dynamic parameters.

TABLE IV Results of Different Algorithms for New England 10-machine 39-bus System

Algorithm	σ_r	ζ_m (%)	Δσ_r	Δζ_m (%)	ΔL (%)	Computation time (s)
Proposed	0.19	3.01	0.12	1.94	5.66	31.14
SQP-GS [10]	0.20		0.16		3.17	40.60
SA [16]	0.32	5.00	0.24	4.22	4.46	86.39

TABLE V Computing Resources Used in [10], [16], and Proposed Algorithm

Algorithm

Machine

Runtime environment

Proposed

Dell Precision T7920 with a

2.4 GHz CPU and 16 GB RAM

Visual Fortran 11.0

SQP-GS [10]

Dell Precision T5810 with a

3.5 GHz CPU and 64 GB RAM

CPLEX 12.60

SA [16]

1.7 GHz Intel Xeon CPU and

32 GB RAM

MATLAB 9.2

B. A Modified Practical 68-machine 2395-bus System

A modified practical 68-machine 2395-bus system shown in Fig. 4 is employed to test the scalability of the proposed algorithm.

Fig. 4 A modified practical 68-machine 2395-bus system.

The system consists of 66 synchronous machines and 2 synchronous compensators, where 57 generators are rescheduled to improve the small-signal stability. There are 2395 buses, 5488 transmission lines, 1496 transformers, and 458 constant impedance loads. The generators are described by the sixth-order model. Except for the equivalent generator, all the generators are equipped with excitation system, turbine governor, and power system stabilizer [

31]. The number of eigenvalues is 1220. The parameters of SSSC are chosen as

σ_{C} = - 0.03

ζ_{T} = 4 %

, and

ζ_{C} = 6 %

1) Comparison of OPF-SSSC1 and OPF-SSSC2

Results of OPF-SSSC1 and OPF-SSSC2 are provided by Table VI, where the eigenvalues are computed by the QR algorithm. Although the augmented Lagrangian function value of the OPF-SSSC2 is 955.29 higher than that of the OPF-SSSC1, the minimum damping ratio increases from 3.98% to 4.17%. In addition, the OPF-SSSC2 significantly reduces the number of iterations and the computation time compared with the OPF-SSSC1.

TABLE VI Results of OPF-SSSC1 and OPF-SSSC2 for Modified Practical 68-machine 2395-bus System

OPF-SSSC	L	ζ_m (%)	l	$k_{t}$	$k_{r}$	Computation time (s)
OPF-SSSC1	115677.92	3.98	23	233	7	8821.81
OPF-SSSC2	116633.21	4.17	10	106	1	4137.84

2) Variable Selection for Eigenvalue Estimation

Since the modified practical 68-machine 2395-bus system is relatively large, it is necessary to reduce the number of state variables in x. The thresholds are set to be $ϖ = 10 %$ and $S = 0.25$ . Taking the second iteration of the OPF-SSSC2 as an example, the logic of variable selection can be explained. According to the results of the second iteration, there are four modes with damping ratios smaller than $ζ_{C}$ . The eigenvalues and damping ratios of the four modes, $λ_{1108}$ , $λ_{1004}$ , $λ_{1011}$ , and $λ_{976}$ , are shown in Table VII. To alleviate the eigenvalue oscillations, these four modes must be considered in the optimization process. Therefore, the dominant variables are selected by $w_{i}^{(2)} Δ x_{i}$ of damping ratios ( $ζ_{1108}$ , $ζ_{1004}$ , $ζ_{1011}$ , and $ζ_{976}$ ) to eigenvalue estimation.

TABLE VII Dominant Variable Serial Numbers of Each Mode and All Considered Modes After the Second Iteration for Modified Practical 68-machine 2395-bus System

No.	λ	ζ (%)	Dominant variable serial number
No.	λ	ζ (%)	Each mode	All considered mode
1107, 1108	$- 0.1432 \pm$ j3.4817	4.11	40, 84, 52, 38, 28, 18	40, 84, 52, 88, 62, 60, 28, 18, 38, 2, 68, 10
1003, 1004	$- 0.2772 \pm$ j5.4111	5.11	84, 38, 52, 60, 40, 10, 28, 18
1010, 1011	$- 0.2235 \pm$ j4.7853	4.66	40, 84, 52, 88, 62, 60, 28, 18, 38, 2, 68, 10
975, 976	$- 0.3954 \pm$ j7.5428	5.23	40, 60, 84, 88, 28, 52

In the proposed algorithm, all state variables of the modified practical 68-machine 2395-bus system are $[P_{G 1}, V_{G 1},$ $P_{G 2}, V_{G 2}, \dots, P_{G 57}, V_{G 57}]$ and the corresponding serial numbers are [ $1, 2, 3, 4, \dots, 113, 114$ ]. Figure 5 shows all $w_{i}^{(2)} Δ x_{i}$ of $ζ_{1108}$ , $ζ_{1004}$ , $ζ_{1011}$ , and $ζ_{976}$ . It can be observed that the w_i⁽²⁾ $Δ x_{i}$ of $ζ_{1108}$ and $ζ_{1011}$ are on the same order of magnitude, while the $w_{i}^{(2)} Δ x_{i}$ of $ζ_{1004}$ and $ζ_{976}$ are on a smaller order of magnitude. For different modes, only a small number of $w_{i}^{(2)} Δ x_{i}$ has a significant impact. According to the Algorithm 1 in Section III-D, dominant variable serial numbers for each mode and all modes are shown in columns 4 and 5 of Table VII. It can be observed that different modes have some same dominant variables. For the four modes, 12 variables are selected for eigenvalue estimation.

Fig. 5 All w_i⁽²⁾Δx_i of $ζ_{1108}$ , $ζ_{1004}$ , $ζ_{1011}$ , and $ζ_{976}$ . (a) $ζ_{1108}$ . (b) $ζ_{1004}$ . (c) $ζ_{1011}$ .(d) $ζ_{976}$ .

3) Effectiveness of Hybrid Algorithm with Variable Selection

The performance comparison of the hybrid algorithm with variable selection is shown in Table VIII. As shown in row 3 of Table VI and rows 2 and 6 of Table VIII, when $S = 0.15$ , the hybrid algorithm with or without variable selection can reduce the computation time. It can also be observed from rows 2 and 6 of Table VIII that the computation time of the hybrid algorithm with variable selection is reduced by 476.87 s compared with the hybrid algorithm without variable selection. This is because eigenvalue estimation by the sensitivity approach only needs to compute the sensitivities with respect to the dominant variables. When S takes different values, it can be observed from rows 2 to 5 of Table VIII that although the hybrid algorithm with variable selection increases the number of iterations, it can effectively reduce the computation time. In addition, when S takes the value of 0.35 or 0.45, the number of iterations and the optimization results are the same, but the computation time difference is 62.46 s. This may be due to the different numbers of dominant variables used for eigenvalue estimation.

TABLE VIII Performance Comparison of Hybrid Algorithm with Variable Selection for Modified Practical 68-machine 2395-bus System

Variable	S	L	ζ_m (%)	l	$k_{t}$	$k_{r}$	$k_{Q R}$	$k_{E S}$	Computation time (s)
Dominant variables	0.15	116813.51	4.25	11	124	1	34	5	3115.23
	0.25	116847.25	4.33	10	116	1	27	5	2720.37
	0.35	116473.37	4.24	13	142	1	26	3	3101.72
	0.45	116473.37	4.24	13	142	1	26	3	3039.26
All variables	0.15	116724.37	4.24	11	106	1	36	7	3592.10

4) Comparison with Partial Eigenvalue Approach

In the proposed algorithm, the eigenvalues are computed using the hybrid algorithm with variable selection that combines accurate eigenvalue computation and eigenvalue estimation. In addition to the QR algorithm, partial eigenvalue approaches, e.g., Krylov-Schur approach [

32] and IRA approach [33], can be used for accurate eigenvalue computation. Taking the IRA approach as an example, simulation results are listed in Table IX.

TABLE IX Simulation Results of Accurate Eigenvalues Calculated by IRA Approach

S	L	ζ_m (%)	l	$k_{t}$	$k_{r}$	$k_{I R A}$	$k_{E S}$	Computation time (s)
0.15	116812.71	4.25	11	124	1	34	5	2837.57
0.45	116473.97	4.24	13	142	1	26	3	2826.66

From Tables VIII and IX, it can be observed that whether the QR algorithm or the IRA approach is used to compute eigenvalues, the influence on the optimization results is small and can even be ignored. This is because the residual norms ${‖A x_{i} - λ_{i} x_{i}‖}_{2}$ of the QR algorithm and the IRA approach are about $10^{- 6}$ and $10^{- 10}$ , respectively. It can also be observed from Tables VIII and IX that when S is taken as 0.15 and 0.45, the computation time is reduced by 277.66 s and 212.60 s, respectively. This is due to the fact that the QR algorithm takes 12.23 s to compute eigenvalues and left and right eigenvectors, while the IRA approach takes only 4.06 s. The simulation results show that both the QR algorithm and the IRA approach can be used for accurate eigenvalue computation in the case study of this paper. Compared with the QR algorithm, the required computation time of the proposed algorithm can be reduced by the IRA approach. For the systems with more generators, partial eigenvalue approaches are more advantageous and necessary. Furthermore, in our experiments, the computation of sensitivities with respect to dominant variables is still the most computationally expensive routine after using the variable selection approach, which can be accelerated by using parallel calculation techniques.

C. Effectiveness of Security Constrained OPF with SSSC2

Taking the New England 10-machine 39-bus system as an example, the OPF-SSSC2 under the normal operating condition is extended to the security constrained OPF with SSSC2 (SCOPF-SSSC2). The following 14 scenarios are taken into account. 1) Scenario 1: base case. 2) Scenario 2: lines 3-18 and 25-26 are out of service. 3) Scenario 3: lines 16-17 and 4-14 are out of service. 4) Scenario 4: line 6-11 is out of service. 5) Scenario 5: 360 MW load increase. 6) Scenario 6: lines 16-17, 4-14, and 25-26 are out of service. 7) Scenario 7: lines 16-17, 4-14, 25-26, and 1-39 are out of service. 8) Scenario 8: line 21-22 is out of service. 9) Scenario 9: lines 9-39 is out of service. 10) Scenario 10: $- 30 %$ loading. 11) Scenario 11: $+ 15 %$ loading. 12) Scenario 12: $+ 20 %$ loading. 13) Scenario 13: $- 20 %$ loading. 14) Scenario 14: $+ 50 %$ loading at buses 16 and 21 and lines 21-22 are out of service. Ramping limits $R_{G, i}^{u p} = R_{G, i}^{d o w n} = (P_{G, i}^{m a x} - P_{G, i}^{m i n}) / 60$ p.u./min [

13]. Critical modes of OPF-SSSC2 and SCOPF-SSSC2 are provided by Table X.

TABLE X Critical Modes of OPF-SSSC2 and SCOPF-SSSC2

Scenario	Critical mode
Scenario	OPF-SSSC2	SCOPF-SSSC2
1	$- 0.0700 \pm j 6.5416$	$- 0.0677 \pm j 6.5235$
2	$- 0.0739 \pm j 2.6772$	$- 0.0314 \pm j 2.5473$
3	$0.3691 \pm j 1.9953$	$- 0.0188 \pm j 2.3533$
4	$- 0.0708 \pm j 6.5452$	$- 0.0695 \pm j 6.5290$
5	$- 0.0699 \pm j 6.5409$	$- 0.0674 \pm j 6.5213$
6	$0.3622 \pm j 1.9962$	$- 0.0205 \pm j 2.3546$
7	$0.4069 \pm j 1.7608$	$- 0.0054 \pm j 2.1619$
8	$- 0.0737 \pm j 6.5280$	$- 0.0027 \pm j 2.5332$
9	$- 0.0710 \pm j 6.5410$	$- 0.0692 \pm j 6.5175$
10	$0.6800 \pm j 1.7742$	$- 0.0754 \pm j 6.4808$
11	$- 0.0708 \pm j 6.5538$	$- 0.0655 \pm j 6.5285$
12	$- 0.0717 \pm j 6.5623$	$- 0.0629 \pm j 2.5556$
13	$0.2695 \pm j 2.1696$	$- 0.0785 \pm j 6.5391$
14	$- 0.0734 \pm j 6.5323$	$- 0.0716 \pm j 6.5139$

The critical mode for Scenario 1 in column 2 of Table X is obtained from the OPF-SSSC2 in Section III-A. Based on the results of Scenario 1, the small-signal stability analysis of other scenarios is performed sequentially to obtain the critical mode. Critical modes in column 3 of Table X are obtained from the SCOPF-SSSC2 in Section III-F. As shown in column 2 of Table X, five scenarios are unstable under small disturbances, and the most unstable scenario is Scenario 10. Compared with the OPF-SSSC2, the value of augmented Lagrangian function of the SCOPF-SSSC2 increases from 40308.18 to 40469.41, but all considered scenarios are stable under small disturbances. From columns 2 and 3 of Table X, it can also be observed that the critical mode of Scenario 10 has changed from one mode to another, and the critical modes of Scenarios 8 and 13 have also changed. To alleviate the oscillation of constrained eigenvalues, it is useful to extend the proposed algorithm considering both the normal and contingency operating conditions.

V. Conclusion

A practical algorithm for the SSSC-OPF problem is proposed in this paper. The effectiveness of the proposed algorithm is validated on the New England 10-machine 39-bus system and a modified practical system. The unique features of the proposed algorithm are provided as follows.

1) The SSSC adopted in this paper not only reflects the dynamic performance requirement of the system, but also fully considers the effect of steady-state changes on eigenvalues. The simultaneous solution formula of the SSSC-OPF model is also established and solved by the quasi-Newton approach.

2) Different penalty factors are set for SSSCs according to the stabilization degree of constrained eigenvalues, which effectively alleviates the eigenvalue oscillations when simultaneously solving the SSSC-OPF problem.

3) The solution time of the SSSC-OPF problem is significantly reduced by the proposed algorithm, which includes variable selection for eigenvalue estimation and strategy selection for eigenvalue computation.

The SSSC-OPF model and algorithm for power systems with high penetration of renewable energy will be developed in future work.

References

H. W. Dommel and W. F. Tinney, “Optimal power flow solutions,” IEEE Transactions on Power Apparatus and Systems, vol. PAS-87, no. 10, pp. 1866-1876, Oct. 1968. [Baidu Scholar]

R. D. Zimmerman, C. E. Murillo-Sánchez, and R. J. Thomas, “Matpower: steady-state operations, planning, and analysis tools for power systems research and education,” IEEE Transactions on Power Systems, vol. 26, no. 1, pp. 12-19, Feb. 2011. [Baidu Scholar]

E. Vaahedi, Y. Mansour, J. Tamby et al., “Large scale voltage stability constrained optimal var planning and voltage stability applications using existing OPF/optimal var planning tools,” IEEE Transactions on Power Systems, vol. 14, no. 1, pp. 65-74, Feb. 1999. [Baidu Scholar]

Y. Fu, X. Zhang, L. Chen et al., “Analytical representation of data-driven transient stability constraint and its application in preventive control,” Journal of Modern Power Systems and Clean Energy, vol. 10, no. 5, pp. 1085-1097, Sept. 2022. [Baidu Scholar]

Y. Lin, X. Zhang, J. Wang et al., “Voltage stability constrained optimal power flow for unbalanced distribution system based on semidefinite programming,” Journal of Modern Power Systems and Clean Energy, vol. 10, no. 6, pp. 1614-1624, Nov. 2022. [Baidu Scholar]

Y. Yu, Y. Liu, C. Qin et al., “Theory and method of power system integrated security region irrelevant to operation states: an introduction,” Engineering, vol. 6, no. 7, pp. 754-777, Jul. 2020. [Baidu Scholar]

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

C. Li, H. D. 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]

P. Li, H. Wei, B. Li et al., “Eigenvalue-optimisation-based optimal power flow with small-signal stability constraints,” IET Generation, Transmission & Distribution, vol. 7, no. 5, pp. 440-450, May 2013. [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]

H. A. Seyed, R. Abdorreza, and T. B. Samad, “Optimal re-dispatch of generating units ensuring small signal stability,” IET Generation, Transmission & Distribution, vol. 14, no. 18, pp. 3692-3701, Aug. 2020. [Baidu Scholar]

J. Condren and T. 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]

R. Zarate-Minano, F. Milano, and A. 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]

S. Kim, A. Yokoyama, Y. Takaguchi et al., “Small-signal stability-constrained optimal power flow analysis of multiterminal VSC-HVDC systems with large-scale wind farms,” IEEJ Transactions on Electrical and Electronic Engineering, vol. 14, no. 7, pp. 1033-1046, Jul. 2019. [Baidu Scholar]

J. Xing, C. Chen, and P. Wu, “Optimal active power dispatch with small-signal stability constraints,” Electric Power Components and Systems, vol. 38, no. 9, pp. 1097-1110, Sept. 2010. [Baidu Scholar]

Y. Li, G. Geng, Q. Jiang et al., “A sequential approach for small signal stability enhancement with optimizing generation cost,” IEEE Transactions on Power Systems, vol. 34, no. 6, pp. 4828-4836, Nov. 2019. [Baidu Scholar]

J. Liu, Z. Yang, J. Zhao et al., “Explicit data-driven small-signal stability constrained optimal power flow,” IEEE Transactions on Power Systems, vol. 37, no. 5, pp. 3726-3737, Nov. 2022. [Baidu Scholar]

P. Pareek and H. D. Nguyen, “A convexification approach for small-signal stability constrained optimal power flow,” IEEE Transactions on Control of Network Systems, vol. 8, no. 4, pp. 1930-1941, Nov. 2021. [Baidu Scholar]

D. Pullaguram, R. Madani, T. Altun et al., “Small-signal stability-constrained optimal power flow for inverter dominant autonomous microgrids,” IEEE Transactions on Industrial Electronics, vol. 69, no. 7, pp. 7318-7328, Jul. 2022. [Baidu Scholar]

L. Shi, C. Wang, L. Yao et al., “Optimal power flow solution incorporating wind power,” IEEE Systems Journal, vol. 6, no. 2, pp. 233-240, Jun. 2012. [Baidu Scholar]

Y. Yang, Y. Luo, and L. Yang, “Small-signal stability-constrained optimal power flow model based on BP neural network algorithm,” Sustainability, vol. 14, no. 20, p. 13386, Nov. 2022. [Baidu Scholar]

X. Wang, Y. Song, and M. Irving, Modern Power Systems Analysis. New York: Springer, 2010. [Baidu Scholar]

S. Frank, I. Steponavice, and S. Rebennack, “Optimal power flow: a bibliographic survey II,” Energy Systems, vol. 3, pp. 259-289, Mar. 2012. [Baidu Scholar]

A. M. H. Rashed and D. H. Kelly, “Optimal load flow solution using Lagrangian multipliers and the Hessian matrix,” IEEE Transactions on Power Apparatus and Systems, vol. PAS-93, no. 5, pp. 1292-1297, Sept. 1974. [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]

J. Nocedal and S. J. Wright, Numerical Optimization. New York: Springer, 2006. [Baidu Scholar]

W. M. Lebow, R. Rouhani, R. Nadira et al., “A hierarchical approach to reactive volt ampere (VAR) optimization in system planning,” IEEE Transactions on Power Apparatus and Systems, vol. PAS-104, no. 8, pp. 2051-2057, Aug. 1985. [Baidu Scholar]

A. Monticelli, M. V. F. Pereira, and S. Granville, “Security-constrained optimal power flow with post-contingency corrective rescheduling,” IEEE Transactions on Power Systems, vol. 2, no. 1, pp. 175-180, Feb. 1987. [Baidu Scholar]

A. M. Sasson, F. Viloria, and F. Aboytes, “Optimal load flow solution using the hessian matrix,” IEEE Transactions on Power Apparatus and Systems, vol. PAS-92, no. 1, pp. 31-41, Jan. 1973. [Baidu Scholar]

M. A. Pai, Energy Function Analysis for Power System Stability. New York: Springer, 1989. [Baidu Scholar]

Z. Wu and X. Zhou, “Power system analysis software package (PSASP) - an integrated power system analysis tool,” in Proceedings of International Conference on Power System Technology, Beijing, China, Aug. 1998, pp. 7-11. [Baidu Scholar]

Y. Li, G. Geng, and Q. Jiang, “An efficient parallel Krylov-Schur method for eigen-analysis of large-scale power systems,” IEEE Transactions on Power Systems, vol. 31, no. 2, pp. 920-930, Mar. 2016. [Baidu Scholar]

R. B. Lehoucq, D. C. Sorensen, and C. Yang, ARPACK Users’ Guide: Solution of Large Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods. Philadelphia: SIAM, 1998. [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