# Available Transfer Capability Calculation Constrained with Small-Signal Stability Based on Adaptive Gradient Sampling

Authors: Peijie Li; Ling Zhu; Xiaoqing Bai; Hua Wei
Journal: Complexity
 (2020)

---

## Abstract

Due to the nonsmoothness of the small-signal stability constraint, calculating the available transfer capability (ATC) limited by small-signal stability rigorously through the nonlinear programming is quite difficult. To tackle this challenge, this paper proposes a sequential quadratic programming (SQP) method combined with gradient sampling (GS) in a dual formulation. The highlighted feature is the sample size of the gradient changes dynamically in every iteration, yielding an adaptive gradient sampling (AGS) process. Thus, the computing efficiency is greatly improved owing to the decrease and the parallelization of gradient evaluation, which dominates the computing time of the whole algorithm. Simulations on an IEEE 10-machine 39-bus system and an IEEE 54-machine 118-bus system prove the effectiveness and high efficiency of the proposed method.


---

## Body

## 1. Introduction

ATC is determined by a number of operating limits such as thermal limit, voltage limit, and transient stability limit [1, 2]. As the power system interconnects sustainably and the renewable energy sources penetrate increasingly, small-signal stability becomes a more significant limiting factor of determining power transfer capabilities, which may cause the occurrence of low-frequency oscillation before reaching the power transmission limit determined by traditional constraints. Therefore, the calculation of ATC considering small-signal stability is a problem to be solved urgently. Due to the nonsmoothness of small-signal stability constraint, the calculation of ATC with small-signal stability constraint (SSSC-ATC) has become a challenge. For there is no common solution to this kind of problem, a few people have made some exploration of it.Chung et al. [3] propose a numerical-sensitivity-based rescheduling method to calculate the ATC subject to the small-signal stability constraint under a set of contingencies. The method is easy to handle; however, it cannot guarantee the convergency because it ignores the nonsmoothness of small-signal stability constraint calculation essentially.Othman and Busan [4] introduce a straightforward approach based on a normalized participation factor to identify critical generators for rescheduling so that a secure power transfer during the outage of critical is calculated. Compared with [3], its major contribution is that it formulates a closed-form of eigenvalue sensitivity with respect to power generation in terms of participation factors. However, the formulation makes a significant approximation as the derivative for power generation with respect to rotor angle is formulated based on a simplified two-order machine model. Besides, it does not address the nonsmoothness of small-signal stability in ATC calculation.Recently, Li et al. [5] propose an SQP method combined with GS to solve the optimal power flow problem with a nonsmooth small-signal stability constraint (SSSC-OPF). The GS procedure samples the gradients of the nonsmooth constraint within the neighborhood to yield a subgradient, thereby ensuring that good search directions are produced in nonsmooth regions. It can guarantee the SSSC-OPF is globally and efficiently convergent to stationary points.Carrying on its idea, this paper makes a step forward in SSSC-ATC calculation field with general applicability and practical significance. The main contributions of our work includes three aspects. Firstly, an AGS process is proposed, in which the sample size of the gradient changes dynamically during every iteration. Thus, the quantity of gradient evaluation is greatly decreased, which dominates the computing time of the whole algorithm, showing the superiority in sampling efficiency compared with SQP-GS [5]. Secondly, the use of a dual formulation of the QP subproblem clearly reflects the idea of gradient sampling when solving the nonsmooth constraints. Thirdly, as the sampling process is independent and random, a process of parallel sampling is applied to the simulation, greatly reducing the time-consumed by the algorithm which is shown in the test of IEEE 54-machine 118-bus system.Section2 describes a model of the ATC calculation with a small-signal stability constraint. The SQP-AGS algorithm applied to the proposed SSSC-ATC model is discussed in Section 3 in detail. Section 4 is case studies, and the conclusion is drawn in Section 5.


## 2. The Model of ATC Calculation with a Small-Signal Stability Constraint

Mathematically, ATC is defined as(1)ATC=TTC−TRM−CBM−ETC,where TTC is the total transfer capacity over the interconnected transmission network. TRM is the transmission reliability margin. CBM is the capacity benefit margin, and ETC is the existing transmission commitments. TRM and CBM are usually considered as constants or a percentage of TTC, and ETC can be easily obtained as the sum of the existing commitments. Thus, the calculation of TTC is the key to obtain ATC. An optimal power flow (OPF) model is formulated to the conventional ATC calculation [6]. Based on the model, the SSSC-ATC includes the following.(1)
The objective value:(2)minf•=−∑i,j∈SlinkPij,whereSlink is the set of tie lines between the source side and the receiving side. Pij indicates the active power of the tie line i,j, which is defined as Pij=Vi2Yijcosαij−ViYijVjcosδij. Vi represents the voltage amplitude of the ith bus. Yij and αij are the admittance amplitude and the phase angle of line i,j, respectively. δij=θi−θj−αij, in which θi is the phase angle of the ith bus.(2)
For all buses in system denoted as collectionSB, the power flow equations for bus i∈SB is(3)PGi−PLi−∑j∈iViYijVjcosδij=0,QGi−QLi−∑j∈iViYijVjsinδij=0,wherePGi is the active power output of the ith generator, QGi is the reactive power output of the ith generator correspondingly, and PLi and QLi are the active power load and the reactive power load of the ith bus, respectively.(3)
Technical constraints include(4)V¯i≤Vi≤V¯i,i∈SB,P¯Gi≤PGi≤P¯Gi,i∈SG,Q¯Gi≤QGi≤Q¯Gi,i∈SG,Sij2≤S¯ij2,i,j∈Sline,whereSG is the collection of all the generator buses and Sline is the collection of all lines. Sij represents the apparent power of line i,j, and ⋅¯ and ⋅¯, respectively, denote the upper and lower limits.(4)
Initial equations which link the variables of the conventional ATC with the state and algebraic variables of the small-signal stability model:(5)SIPG,QG,V,θ,x=0,wherex is the state variable vector including the rotor angle vector δ, the d-axis and q-axis component vectors of the internal voltage are Ed′ and Eq′, the d-axis and q-axis component vectors of the internal current are Id′ and Iq′, and the excitation output voltage vector is Efd.(5)
Small-signal stability constraint:(6)ηPG,QG,V,θ,x≤η¯,where η is the spectral abscissa of the system, which is the largest one of the real parts of the system’s eigenvalues λ. η¯ is an upper limit of spectral abscissa to represent a stability margin, which can be determined based on offline stability studies.The spectral abscissaη in small-signal stability constraint is a nonsmooth function [7]. Thus, it failed to solve it by the interior point method (IPM), which is the workhorse in solving smooth constrained optimization problems. Fortunately, the spectral abscissa function has been proved to be locally Lipschitz and continuously differentiable on open dense subsets D of Rn, which means that it is continuously almost everywhere and its gradient can be easily obtained where it is defined by calculating spectral abscissa derivative.


## 3. A Sequential Quadratic Programming Algorithm Combined with Gradient Sampling

There is no general method to solve the optimization problem with nonsmooth nonlinear constraints these years. Recently, a gradient sampling (GS) method is proposed for minimizing an objective function that is locally Lipschitz continuous on an open dense subset ofRn, which makes a step forward in nonsmooth optimization field [8]. Curtis and Overton [9] further generalized this idea to the optimization problems with nonsmooth nonlinear constraints or objective function by introducing the SQP algorithm framework.


### 3.1. The Gradient Sampling Theory with Unconstrained Minimization

The GS is basically a stabilized steepest descent algorithm. Considering an unconstrained optimization problem of the form(7)minx∈Rnfx,defined in an open, dense subset D∈Rn, the objective function f is locally Lipschitz continuous which is nonsmooth and/or nonconvex.As Figure1 shows, the gradients of smooth points like x1 and x2 are easily obtained and unique. But it is unable to obtain the gradient at a nonsmooth point x0. Moreover, the subgradient of nonsmooth point exists and is not unique. According to the definition,(8)∂f≔ufx≥fx0+u⊤x−x0,∀x∈Rn⟶R,u∈∂f denotes the subgradient at the nonsmooth point x0. The set of subgradients u is called the subdifferential of f at x0.Figure 1
The subgradients of a nonsmooth point.The GS technique approximates the subgradients of the objective function through the gradients sampled randomly nearx0 as(9)∂fx0≈∑i∈Pyki∇fxki∀i:yki>0∧∑i∈Pyki=1,where P is the amount of the sampled gradients. Defining that Gk≔∇fxk0,∇fxk1,…,∇fxkp is a set of sampled gradients evaluating from a set of random and independent sampling points Xk≔xk0,xk1,…,xkp⊂Bϵx∩D near x0, where(10)Bϵx0≔xx−x0≤ϵ,is the closed ball with a sample radius ϵ centered at x0. The function f is locally Lipschitz and almost continuously differentiable.The expression of∂fx0 is actually a convex hull of sampled gradients, which is the set of all convex combinations of the gradients. In the convex combination, each sampled gradient ∇fxki in Gk is assigned a nonnegative weight or coefficient yki that all of them are summed to one. In such a way, the set of subgradients of x0 is derived by changing yki into different values.It is known that the conventional way to find the stationary point of a smooth functionf is to find the point whose gradient values are zero, and the gradient value is generally taken as a key criterion to decide the descent direction of a smooth optimization problem. Corresponding to it, the Clarke stationarity is a criterion for determining the stationary point in a nonsmooth function, which points out that a point x0 is Clarke stationary if(11)0∈∂fx.Figure1 marks out a zero-valued subgradient u1 determining that x0 is a Clarke stationary point. The Clarke stationarity is the key to the approach of finding a descent direction.Hence, the descent direction of nonsmooth optimization can be determined as the subgradient which has the shortest distance to origin among all the subgradients. To find the subgradientuk at the kth iteration, a following quadratic programming (QP) model is established:(12)miny∈Rp+1Gkyk22,s.t.eTyk=1,yk≥0,where yk is the coefficient vector at the kth iteration. The 2-norm in the objective function (12) represents the shortest distance from the set of subgradients to the origin. Solving this model, the coefficient vector yk can be obtained. Then, the subgradient uk can be calculated according to (9). The vector dk=−uk/uk is defined as an approximate steepest descent direction. If uk=0, the stationary point is found.


### 3.2. Sequential Quadratic Programming Algorithm Combined with Gradient Sampling for Nonsmoothly Constrained Optimization

Similar to the QP model used in the GS algorithm to calculate a search direction, a sequential quadratic programming algorithm (SQP) framework is built by centering on a QP subproblem to determine a search direction of nonsmooth constraint optimization problems. Generally, the limitation of the wide-used SQP algorithm lies in that it can only solve the smooth optimization problems but fail for nonsmooth optimization problems. However, combined with the GS algorithm, the SQP-GS algorithm can be used in in many nonsmooth problems even if the objective function is not locally Lipschitz continuous, which is basically in the following form:(13)minxfx,s.t.hx=0,gx≤0.The SQP-GS algorithm solves the model above in each iteration mainly by determining a search directiondk centering on solving a regularized QP dual subproblem with the following form:(14)mindk,yo,yh,yg−fxk−hxkhxk−gxk⊤yfyhyg+dk⊤,yf⊤,yh⊤,yg⊤Hkdkyfyhyg,(15)s.t.Hkdk+∇fxfTyf+∇hxh−∇hxh∇gxg⊤yhyg=0,∀xf∈ℬϵ,kf,∀xh∈ℬϵ,kh,∀xg∈ℬϵ,kg,(16)eTyf=ρ,(17)eTyh=1,(18)eTyg=1,(19)yf≥0,(20)yh≥0,(21)yg≥0,where ρ is a penalty parameter and Hk is the approximate Hessian of the Lagrangian optimization model. Incorporating the ideas from the GS algorithm, ℬϵ,kf, ℬϵ,kh, and ℬϵ,kg are sets of independent and identically distributed random points uniformly sampled from (10) with a sample size of p:(22)ℬϵ,kf≔xk,0f,xk,1f,…,xk,pf,wherexk,0f≔xk,(23)ℬϵ,khi≔xk,0hi,xk,1hi,…,xk,phi,wherexk,0hi≔xk,(24)ℬϵ,kgj≔xk,0gj,xk,1gj,…,xk,pgj,wherexk,0gj≔xk.Consistent with the GS strategy,yf=yk,0f,yk,1f,…,yk,pf, yh=yk,0hi,yk,1hi,…,yk,phi, and yg=yk,0gj,yk,1gj,…,yk,pgj are approximated minimum-norm elements of the subdifferential for objective function, equality and inequality constraints, respectively, which are obtained by solving the subproblem (14)–(21) at xk.As a significant step to updateHk during each iteration, we introduce the limited-memory Broyden–Fletcher–Goldfrab–Shanno (L-BFGS) method which is an attempt to achieve fast convergence. We define(25)Gk≔∑n=1pyk,nf∇fxk,nf+∑m=1j∑n=1pyk,nhm∇hxk,nhm+∑m=1j∑n=1pyk,ngm∇gxk,ngm,as the approximation of the minimum-norm element of the differential of (13) Lagrange function, relating to the choice of Hk made during each iteration. We initialize H1⟵I at the 1th iteration and perform the update as(26)Hl⟵Hl−1−Hl−1slsl⊤Hl−1⊤sl⊤Hl−1sl+ylyl⊤sl⊤yl,for l=2,…,k. Here, sl≔xl−xl−1 and yl≔Gl−Gl−1 are the displacements in x and in the approximations of the minimum-norm elements of the differential of (13) Lagrange function.In order to reflect the iteration process, consider the penalty function with the penalty parameterρ as the following form:(27)ϕρx≔ρfx+σx,where σx is an infeasibility vector at the kth iteration defined as(28)σkxk≔hxkmaxgxk,0,which is also an indicator in the step of parameter update.Furthermore, let us define the local model(29)qρdk;xk,ℬϵ,k,Hk≔ρmaxx∈ℬϵ,kffxk+∇fx⊤dk+∑i=1mhmaxx∈ℬϵ,khihixk+∇hix⊤dk+∑i=1mgmaxx∈ℬϵ,kgimaxgixk+∇gix⊤dk,0+12dk⊤Hkdk.For a computeddk, the model reduction is defined to be △qρdk;xk,ℬϵ,k,Hk≔qρ0;xk,ℬϵ,k,Hk−qρdk;xk,ℬϵ,k,Hk, which values(30)△qρdk;xk,ℬϵ,k,Hk=ϕρx−qρdk;xk,ℬϵ,k,Hk.The model reduction is nonnegative becausedk=0 yields an objective value, and △qρdk;xk,ℬϵ,k,Hk can be zero only if xk is ϵ-stationary. So a sufficiently small reduction in △qρ⋅;xk,ℬϵ,k,Hk indicates that a decrease in ϵ is appropriate.However, GS produces the approximateϵ-steepest descent directions by evaluating gradients of the objective function and constraint functions at p randomly generated points during each iteration, resulting in a high computational cost. As the QP subproblem data are computed afresh for every iteration, the computational effort required to solve each of these subproblems can be significant for the large-scale problem [10]. As a more time-saving way to reduce the computation, a generic adaptive gradient sampling (AGS) algorithm is introduced here in which GS is a special case.The difference between GS and AGS is mainly in the process of sample pointsxk,p∈ℬϵ,k. Initialize a sampling size value p at first, which is the total number of sample points in every iteration. For Xk≔xk,1,…,xk,p (with xk,i=xk for some i and Bk≔Bϵxk∩D, AGS sets(31)Xk⟵Xk−1∩Bk∪xk∪X¯k,where the sample set X¯k≔x¯k,1,…,x¯k,p¯ is composed of p¯ points generated uniformly in Bk. Denote the number of points selected by Xk−1∩Bk is pre. The actual sampling size of kth iteration changes to p¯=p−pre, equaling to the number of points in X¯k.It is apparent to see that AGS needs only to compute the unknown gradients at the updating sample points inXk and change the sample size p adaptively, which cut down some of the workloads in computing the gradients and renewing the QP data at the kth iteration. The complete SQP-AGS algorithm is presented in Algorithm 1.What is worth mentioning is that the sample sizep can be chosen as zero when solving a smooth function for it is unnecessary to sample the function’s gradients at nearby points, in which case the SQP-AGS algorithm reduces to the general SQP algorithm. Especially in a model with smooth and nonsmooth functions, the algorithm samples points of the nonsmooth part only, cutting down a lot of work in the process of sample points and evaluating gradients of the smooth part.Algorithm 1: SQP-AGS algorithm.
(1)
Initialization: choose a sample radius ϵ>0, penalty parameter ρ>0, constraint violation tolerance τ>0, sample size p, line search constant ϖ∈0,1, backtracking constant γ∈0,1, sample radius reduction factor μϵ∈0,1, penalty parameter reduction factor μρ∈0,1, constraint violation tolerance reduction factor μτ∈0,1, infeasibility tolerance νin>0, and stationarity tolerance parameter νs>0. Choose an initial iterate x0∈D and set k=1 and kmax=200.
(2)
Termination conditions check: while k<kmax, if dk<νs and maxσk<νin, output solution and stop.
(3)
Adaptive gradient sampling: generate ℬϵ,kf, ℬϵ,kh, and ℬϵ,kg by (22)–(24) and (31).
(4)
Search direction calculation: choose Hk≻0, solve (14)–(21) to get dk,yo,yh,yg.
(5)
L-BFGS update: update Hk by (26) for next iteration.
(6)
Parameter update: if dk>νsϵ2, then go to step 9. Otherwise, if maxσk≤τ, set τ⟵μττ, but if maxσk>τ, set ρ⟵μρρ. Then, set ϵ⟵μϵϵ, βk⟵0, and go to step 8.
(7)
Line search: set βk as the largest value in the sequence 1,γ,γ2,… such that xk+1⟵xk+βkdk satisfies:ϕρxk+1≤ϕρxk−ϖβkdk and xk+βkdk−xk+1≤minβk,ϵdk
(8)
Iteration increment: set k⟵k+1 and go to step 3.
(9)
end do


### 3.3. Solving SSSC-ATC by SQP-AGS

Combined the characteristic of the proposed SSSC-ATC model, we discuss the SQP-AGS in the specific solution process.Firstly, the AGS technique is only applied to the small-signal stability constraint for it is the only the nonsmooth constraint in the proposed SSSC-ATC model. Therefore, the sample size of the SQP-AGS algorithm reduces to zero in solving the objective function and other smooth constraints.Secondly, considering the different characteristics and capacity range of variables in the model, the sampling radiuses of variables are also specified. In detail, for active powerPG and reactive power QG, the sample radiuses ϵPG and ϵQG are stipulated as 20% of the capacity range. For voltage-type variables (per-unit value) like V and E, the sample radiuses ϵV and ϵE are chosen as 0.03 p.u. For angle variables like θ, the sample radius ϵθ is chosen as 3∘.Thirdly, owing to the sample points are independent and randomly sampled, a multicore parallel computing technology is applied to solve the gradient calculation which occupies the majority of computing work of the whole process. The application of parallel computing technology will further improve the efficiency of SQP-AGS.


## 3.1. The Gradient Sampling Theory with Unconstrained Minimization

The GS is basically a stabilized steepest descent algorithm. Considering an unconstrained optimization problem of the form(7)minx∈Rnfx,defined in an open, dense subset D∈Rn, the objective function f is locally Lipschitz continuous which is nonsmooth and/or nonconvex.As Figure1 shows, the gradients of smooth points like x1 and x2 are easily obtained and unique. But it is unable to obtain the gradient at a nonsmooth point x0. Moreover, the subgradient of nonsmooth point exists and is not unique. According to the definition,(8)∂f≔ufx≥fx0+u⊤x−x0,∀x∈Rn⟶R,u∈∂f denotes the subgradient at the nonsmooth point x0. The set of subgradients u is called the subdifferential of f at x0.Figure 1
The subgradients of a nonsmooth point.The GS technique approximates the subgradients of the objective function through the gradients sampled randomly nearx0 as(9)∂fx0≈∑i∈Pyki∇fxki∀i:yki>0∧∑i∈Pyki=1,where P is the amount of the sampled gradients. Defining that Gk≔∇fxk0,∇fxk1,…,∇fxkp is a set of sampled gradients evaluating from a set of random and independent sampling points Xk≔xk0,xk1,…,xkp⊂Bϵx∩D near x0, where(10)Bϵx0≔xx−x0≤ϵ,is the closed ball with a sample radius ϵ centered at x0. The function f is locally Lipschitz and almost continuously differentiable.The expression of∂fx0 is actually a convex hull of sampled gradients, which is the set of all convex combinations of the gradients. In the convex combination, each sampled gradient ∇fxki in Gk is assigned a nonnegative weight or coefficient yki that all of them are summed to one. In such a way, the set of subgradients of x0 is derived by changing yki into different values.It is known that the conventional way to find the stationary point of a smooth functionf is to find the point whose gradient values are zero, and the gradient value is generally taken as a key criterion to decide the descent direction of a smooth optimization problem. Corresponding to it, the Clarke stationarity is a criterion for determining the stationary point in a nonsmooth function, which points out that a point x0 is Clarke stationary if(11)0∈∂fx.Figure1 marks out a zero-valued subgradient u1 determining that x0 is a Clarke stationary point. The Clarke stationarity is the key to the approach of finding a descent direction.Hence, the descent direction of nonsmooth optimization can be determined as the subgradient which has the shortest distance to origin among all the subgradients. To find the subgradientuk at the kth iteration, a following quadratic programming (QP) model is established:(12)miny∈Rp+1Gkyk22,s.t.eTyk=1,yk≥0,where yk is the coefficient vector at the kth iteration. The 2-norm in the objective function (12) represents the shortest distance from the set of subgradients to the origin. Solving this model, the coefficient vector yk can be obtained. Then, the subgradient uk can be calculated according to (9). The vector dk=−uk/uk is defined as an approximate steepest descent direction. If uk=0, the stationary point is found.


## 3.2. Sequential Quadratic Programming Algorithm Combined with Gradient Sampling for Nonsmoothly Constrained Optimization

Similar to the QP model used in the GS algorithm to calculate a search direction, a sequential quadratic programming algorithm (SQP) framework is built by centering on a QP subproblem to determine a search direction of nonsmooth constraint optimization problems. Generally, the limitation of the wide-used SQP algorithm lies in that it can only solve the smooth optimization problems but fail for nonsmooth optimization problems. However, combined with the GS algorithm, the SQP-GS algorithm can be used in in many nonsmooth problems even if the objective function is not locally Lipschitz continuous, which is basically in the following form:(13)minxfx,s.t.hx=0,gx≤0.The SQP-GS algorithm solves the model above in each iteration mainly by determining a search directiondk centering on solving a regularized QP dual subproblem with the following form:(14)mindk,yo,yh,yg−fxk−hxkhxk−gxk⊤yfyhyg+dk⊤,yf⊤,yh⊤,yg⊤Hkdkyfyhyg,(15)s.t.Hkdk+∇fxfTyf+∇hxh−∇hxh∇gxg⊤yhyg=0,∀xf∈ℬϵ,kf,∀xh∈ℬϵ,kh,∀xg∈ℬϵ,kg,(16)eTyf=ρ,(17)eTyh=1,(18)eTyg=1,(19)yf≥0,(20)yh≥0,(21)yg≥0,where ρ is a penalty parameter and Hk is the approximate Hessian of the Lagrangian optimization model. Incorporating the ideas from the GS algorithm, ℬϵ,kf, ℬϵ,kh, and ℬϵ,kg are sets of independent and identically distributed random points uniformly sampled from (10) with a sample size of p:(22)ℬϵ,kf≔xk,0f,xk,1f,…,xk,pf,wherexk,0f≔xk,(23)ℬϵ,khi≔xk,0hi,xk,1hi,…,xk,phi,wherexk,0hi≔xk,(24)ℬϵ,kgj≔xk,0gj,xk,1gj,…,xk,pgj,wherexk,0gj≔xk.Consistent with the GS strategy,yf=yk,0f,yk,1f,…,yk,pf, yh=yk,0hi,yk,1hi,…,yk,phi, and yg=yk,0gj,yk,1gj,…,yk,pgj are approximated minimum-norm elements of the subdifferential for objective function, equality and inequality constraints, respectively, which are obtained by solving the subproblem (14)–(21) at xk.As a significant step to updateHk during each iteration, we introduce the limited-memory Broyden–Fletcher–Goldfrab–Shanno (L-BFGS) method which is an attempt to achieve fast convergence. We define(25)Gk≔∑n=1pyk,nf∇fxk,nf+∑m=1j∑n=1pyk,nhm∇hxk,nhm+∑m=1j∑n=1pyk,ngm∇gxk,ngm,as the approximation of the minimum-norm element of the differential of (13) Lagrange function, relating to the choice of Hk made during each iteration. We initialize H1⟵I at the 1th iteration and perform the update as(26)Hl⟵Hl−1−Hl−1slsl⊤Hl−1⊤sl⊤Hl−1sl+ylyl⊤sl⊤yl,for l=2,…,k. Here, sl≔xl−xl−1 and yl≔Gl−Gl−1 are the displacements in x and in the approximations of the minimum-norm elements of the differential of (13) Lagrange function.In order to reflect the iteration process, consider the penalty function with the penalty parameterρ as the following form:(27)ϕρx≔ρfx+σx,where σx is an infeasibility vector at the kth iteration defined as(28)σkxk≔hxkmaxgxk,0,which is also an indicator in the step of parameter update.Furthermore, let us define the local model(29)qρdk;xk,ℬϵ,k,Hk≔ρmaxx∈ℬϵ,kffxk+∇fx⊤dk+∑i=1mhmaxx∈ℬϵ,khihixk+∇hix⊤dk+∑i=1mgmaxx∈ℬϵ,kgimaxgixk+∇gix⊤dk,0+12dk⊤Hkdk.For a computeddk, the model reduction is defined to be △qρdk;xk,ℬϵ,k,Hk≔qρ0;xk,ℬϵ,k,Hk−qρdk;xk,ℬϵ,k,Hk, which values(30)△qρdk;xk,ℬϵ,k,Hk=ϕρx−qρdk;xk,ℬϵ,k,Hk.The model reduction is nonnegative becausedk=0 yields an objective value, and △qρdk;xk,ℬϵ,k,Hk can be zero only if xk is ϵ-stationary. So a sufficiently small reduction in △qρ⋅;xk,ℬϵ,k,Hk indicates that a decrease in ϵ is appropriate.However, GS produces the approximateϵ-steepest descent directions by evaluating gradients of the objective function and constraint functions at p randomly generated points during each iteration, resulting in a high computational cost. As the QP subproblem data are computed afresh for every iteration, the computational effort required to solve each of these subproblems can be significant for the large-scale problem [10]. As a more time-saving way to reduce the computation, a generic adaptive gradient sampling (AGS) algorithm is introduced here in which GS is a special case.The difference between GS and AGS is mainly in the process of sample pointsxk,p∈ℬϵ,k. Initialize a sampling size value p at first, which is the total number of sample points in every iteration. For Xk≔xk,1,…,xk,p (with xk,i=xk for some i and Bk≔Bϵxk∩D, AGS sets(31)Xk⟵Xk−1∩Bk∪xk∪X¯k,where the sample set X¯k≔x¯k,1,…,x¯k,p¯ is composed of p¯ points generated uniformly in Bk. Denote the number of points selected by Xk−1∩Bk is pre. The actual sampling size of kth iteration changes to p¯=p−pre, equaling to the number of points in X¯k.It is apparent to see that AGS needs only to compute the unknown gradients at the updating sample points inXk and change the sample size p adaptively, which cut down some of the workloads in computing the gradients and renewing the QP data at the kth iteration. The complete SQP-AGS algorithm is presented in Algorithm 1.What is worth mentioning is that the sample sizep can be chosen as zero when solving a smooth function for it is unnecessary to sample the function’s gradients at nearby points, in which case the SQP-AGS algorithm reduces to the general SQP algorithm. Especially in a model with smooth and nonsmooth functions, the algorithm samples points of the nonsmooth part only, cutting down a lot of work in the process of sample points and evaluating gradients of the smooth part.Algorithm 1: SQP-AGS algorithm.
(1)
Initialization: choose a sample radius ϵ>0, penalty parameter ρ>0, constraint violation tolerance τ>0, sample size p, line search constant ϖ∈0,1, backtracking constant γ∈0,1, sample radius reduction factor μϵ∈0,1, penalty parameter reduction factor μρ∈0,1, constraint violation tolerance reduction factor μτ∈0,1, infeasibility tolerance νin>0, and stationarity tolerance parameter νs>0. Choose an initial iterate x0∈D and set k=1 and kmax=200.
(2)
Termination conditions check: while k<kmax, if dk<νs and maxσk<νin, output solution and stop.
(3)
Adaptive gradient sampling: generate ℬϵ,kf, ℬϵ,kh, and ℬϵ,kg by (22)–(24) and (31).
(4)
Search direction calculation: choose Hk≻0, solve (14)–(21) to get dk,yo,yh,yg.
(5)
L-BFGS update: update Hk by (26) for next iteration.
(6)
Parameter update: if dk>νsϵ2, then go to step 9. Otherwise, if maxσk≤τ, set τ⟵μττ, but if maxσk>τ, set ρ⟵μρρ. Then, set ϵ⟵μϵϵ, βk⟵0, and go to step 8.
(7)
Line search: set βk as the largest value in the sequence 1,γ,γ2,… such that xk+1⟵xk+βkdk satisfies:ϕρxk+1≤ϕρxk−ϖβkdk and xk+βkdk−xk+1≤minβk,ϵdk
(8)
Iteration increment: set k⟵k+1 and go to step 3.
(9)
end do


## 3.3. Solving SSSC-ATC by SQP-AGS

Combined the characteristic of the proposed SSSC-ATC model, we discuss the SQP-AGS in the specific solution process.Firstly, the AGS technique is only applied to the small-signal stability constraint for it is the only the nonsmooth constraint in the proposed SSSC-ATC model. Therefore, the sample size of the SQP-AGS algorithm reduces to zero in solving the objective function and other smooth constraints.Secondly, considering the different characteristics and capacity range of variables in the model, the sampling radiuses of variables are also specified. In detail, for active powerPG and reactive power QG, the sample radiuses ϵPG and ϵQG are stipulated as 20% of the capacity range. For voltage-type variables (per-unit value) like V and E, the sample radiuses ϵV and ϵE are chosen as 0.03 p.u. For angle variables like θ, the sample radius ϵθ is chosen as 3∘.Thirdly, owing to the sample points are independent and randomly sampled, a multicore parallel computing technology is applied to solve the gradient calculation which occupies the majority of computing work of the whole process. The application of parallel computing technology will further improve the efficiency of SQP-AGS.


## 4. Case Studies

The proposed SQP-AGS method is applied to the IEEE 10-machine 39-bus system and the modified IEEE 54-machine 118-bus system to illustrate the effectiveness in solving SSSC-ATC. The full dynamic data of these two systems is extracted from [11, 12], respectively. For both systems, the generators are described by the two-axis model with an IEEE type-I exciter. The loads are modeled as constant power. All the simulations are performed on a workstation with the Intel(R) Xeon(R) E5-2667 v2 processor, which includes 16 cores with the base frequency of 3.30 GHz.The SQP-AGS method is implemented in MATLAB by using the functionquadprog as the QP solver for the subproblem. The eigenvalues and eigenvectors are computed by QR decomposition using the MATLAB function eig. A flat start is used for which all voltage angles are set to be zero, and all the voltage magnitudes are set to be 1.0 p.u., PG=P¯G+P¯G/2 and QG=Q¯G+Q¯G/2. The parameters of SQP-AGS are chosen as ρ=2×10−4 and μρ=0.05. The selection of sample radius ϵ depends on the characteristics of the dependent variables as described in Subsection 3.3. Set μϵ=1×10−3, τ=0.1, μτ=0.8, ϖ=1, and γ=0.8 from [9]. The tolerances νin and νs are set to be 5×10−3 and 2×10−3, respectively.


### 4.1. Comparison of Results under Different Small-Signal Stability Constraints Run from IEEE 10-Machine 39-Bus System

The IEEE 10-machine 39-bus system is used for stability analysis with different small-signal stability constraints, which are three different spectral abscissa constrained values in the SSSC-ATC calculation model above. The specific regional division is shown in Figure2. The generators G2, G3, and G10 which, respectively, correspond to bus 31, 32, and 39 are selected to locate at the receiving side, while the other generators are located at the source side.Figure 2
The regional division of the IEEE 39-bus system.For the voltage magnitude limits, we chooseV¯=1.1 p.u. and V¯=0.9 p.u. for all buses. The results under different spectral abscissa constraints are listed in Table 1. In order to reflect the effect on the SQP-AGS algorithm applied on SSSC-ATC, we set the sample size as p=30 and make a comparison of the following cases:(1)
Case 0: base case, which is a standard ATC calculation without any small-signal stability constraint(2)
Case 1: SSSC-ATC with an abscissa spectrum constraint ofη≤−0.10(3)
Case 2: SSSC-ATC with an abscissa spectrum constraint ofη≤−0.192Table 1
Active output of generators and TTC in IEEE 39-bus system.
Active output (MW)
Active output range (MW)
Case 0
Case 1
Case 2
G1
590.13
1102.59
1106.74
100∼1380
G2
50.00
105.35
136.58
50∼747.5
G3
50.00
677.18
676.12
50∼920.0
G4
862.50
665.66
639.15
50∼862.5
G5
747.50
529.28
564.18
50∼747.5
G6
862.50
671.67
667.62
50∼862.5
G7
862.50
581.69
603.44
50∼862.5
G8
805.00
597.35
513.94
100∼805.0
G9
1035.00
899.08
898.57
100∼1035
G10
402.50
318.83
346.63
50∼402.5
TTC (MW)
1704.99
498.62
465.73
—(i)
Case 0: the classic IPM is used to solve the ATC problem without any small-signal stability constraint. The results are listed in the second column of Table1, which shows the TTC is obtained as 1704.99 MW and the spectral abscissa is calculated to be 2.15, indicating the TTC is overoptimistic due to the fact that the system is not small-signal stable in this operation state. This further indicates that the influence of small-signal stability on ATC calculation needs to be considered. Extracting information from Table 1, it can be seen that the active generation of generators G2 and G3 which are located in the receiving side all reach the lower limit of the active power output in order to ensure the maximum of TTC.
(ii)
Case 1: a spectral abscissa constraintη≤−0.10 is added to the ATC model as the small-signal stability constraint and SQP-AGS algorithm is run to solve it. The results are listed in the third column of Table 1, which shows the TTC is obtained as 498.62 MW, and the system is small-signal stable in this operation state. The result also shows that the generation in the receiving side is greatly increased compared with Case 0, especially as generator G3, possibly due to the larger time inertia coefficient of it, resulting in the reduction of TTC. It reflects that the small-signal stability does constrain the regional transmission directly.
(iii)
Case 2: the spectral abscissa constraint is further reduced toη≤−0.192 in the SSSC-ATC model. The results run from the SQP-AGS algorithm shows that the calculation is still convergent and the TTC is reduced to 465.73 MW.Comparing Case 0, Case 1, and Case 2 will find out a rule that TTC value reduces as the spectral abscissa decreases, which means the stricter restriction effect on ATC calculation when a spectral abscissa constraint gets smaller. The results under spectral abscissa constraints are more reasonable than the one without any small-signal stability constraint so that the ATC in this situation can be a guide to avoid the power oscillation caused by small-signal stability in the system.


### 4.2. Comparison of Different Methods in Calculating SSSC-ATC

Here, we compare the SQP-AGS’s calculation effect on the same problem with a sensitivity method [3] and a participation factor method [4] for IEEE 39-bus system. The results run by the three methods are all restricted to a spectral abscissa constraint as η≤−0.1 and detailedly listed in Table 2.Table 2
Comparison of different methods for the 39-bus system.
SQP-AGS (MW)
Sensitivity method (MW)
Participation factor method (MW)
G1
1102.59
1210.27
1357.84
G2
105.35
556.43
150.46
G3
677.18
407.27
704.98
G4
665.66
732.91
680.85
G5
529.28
599.08
673.28
G6
671.67
416.40
862.50
G7
581.69
764.68
50.00
G8
597.35
715.46
733.05
G9
899.08
790.18
627.21
G10
318.83
50.00
402.50
TTC (MW)
498.60
289.94
262.78The TTC calculated by SQP-AGS is 498.60 MW, while the TTC calculated by the sensitivity method and the participation factor method is 289.94 MW, and the one calculated by the participation factor method is 262.78 MW. All the three methods obtain a TTC value much smaller than which is calculated from the ATC model without small-signal stability constraint. In addition, compared with the result calculated by SQP-AGS, the other two methods’ results are relatively conservative. In other words, it is to say that SQP-AGS can evaluate the SSSC-ATC problem in an optimization way. It further means this ATC result may bring more economic benefits that are more acceptable to the power industry in practice.On the other hand, it is worth mentioning that when the spectral abscissa constraint decreases asη≤−0.15, the sensitivity method and the participation method both cannot converge to obtain a result, while the SQP-AGS can still effectively reduce the spectral abscissa of the system and achieve a good convergent result, as is shown in Table 1. This reveals the superiority of SQP-AGS in solving the SSSC-ATC problem compared with the other two mentioned methods.


### 4.3. Efficiency Tested in IEEE 54-Machine 118-Bus System

To test the efficiency of SQP-AGS, we make some comparisons on a modified version of the IEEE 54-machine 118-bus benchmark system which has 54 synchronous machines with IEEE type-1 exciters, 20 of which are synchronous compensators used only for reactive power support and 15 of which are motors. The specific regional division of IEEE 118-bus system is as follows:(i)
Receiving side: the region includes buses of 1, 4, 6, 8, 10, 12, 15, 18, 19, 24, 25, 26, 27, 31, 32, 34, 36, 40, 42, 46, 49, 54, 55, 56, 59, 61, 62, 65, 66, 69, 70, 72, 73, 113, and 116.(ii)
Source side: the region includes buses of 74, 76, 77, 80, 85, 87, 89, 90, 91, 92, 99, 100, 103, 104, 105, 107, 110, 111, and 112.Lines 74–70, 75–70, 75–69, 77–69, and 81–68 are the tie lines between the two regions. For an ATC model without any small-signal stability constraint, the TTC is obtained as 1410.20 MW and the spectral abscissa isη=1.1, which also means the calculated ATC is overoptimistic, and power oscillation will occur before the transmission between regions reaches this value. Applying SQP-AGS to solve an SSSC-ATC problem with a spectral abscissa constraint η≤−0.1, the spectral abscissa is decreased meeting the constraint requirement after 17 iterations and the result is convergent, obtaining the TTC as 446.08 MW. The spectral abscissa changing with iterations is shown in Figure 3. It is easy to find that the descending trend of spectral abscissa is stable. Moreover, Figure 4 reflects the changes in the termination conditions with the iteration increases. As the trend shows, the infeasibility defined in (28) and the model reduction calculated by (30) both decrease to an acceptable tolerance.Figure 3
Spectral abscissa changes with iterations for the 118-bus system.Figure 4
The termination conditions change with iterations for the 118-bus system.In order to test the efficiency of SQP-AGS, the algorithm is going with a sample size ofp=50. Under the serial computing condition, we make a comparison between SQP-GS and SQP-AGS dealing with the same SSSC-ATC problem. The computing time of them is listed in second and third columns of Table 3.Table 3
Efficiency comparison of the IEEE 118-bus system.
Serial computing
Parallel computing for SQP-AGS
SQP-GS method
SQP-AGS method
Two cores
Four cores
Eight cores
Total time (s)
173.48
66.29
39.88
25.15
16.36
Gradient sampling time (s)
170.40
63.32
36.90
20.24
13.37
QP solver time (s)
2.78
2.67
2.68
2.62
2.69
Others’ time (s)
0.30
0.30
0.30
0.29
0.30
Speedup ratio,Sn
—
—
1.71
2.86
4.05
Efficiency,En
—
—
0.86
0.72
0.51
Iteration
18
17
17
17
17It is clear to see that, for the 118-bus system, applying the serial SQP-GS method to solve an SSSC-ATC problem consumes 173.48s and 18 iterations while applying the serial SQP-AGS needs only 66.29s and 17 iterations, cutting down a lot of time in solving the same model. Particularly, the time-consuming of gradient calculating accounts for the vast majority of the total time because of the repeating calculating processes of small-signal stability calculation and the sampling gradient calculation of nonsmooth constraint. Compared with SQP-GS, SQP-AGS spends much less time on gradient calculation due to the adaptive identification of sampling points which are within the sampling neighborhood of the current iteration point. In each iteration, SQP-AGS effectively utilizes the historical sampling data in the previous iteration and saves a lot of time for gradient calculation at current iteration, which greatly improves the efficiency of calculation. The intuitive sample point changes are shown in Figure5.Figure 5
Number change of new sampling points.Moreover, notice that the sample points for gradient calculation in each iteration are random and independent, we adopt the multicore parallel computing technology of MATLAB to carry out the gradient calculation when applying SQP-AGS. Parallel computing opens the parallel computing pool of specified number of workers through theparpool instruction in MATLAB and uses parfor statement to execute the operation process of the loop body. Testing it on the IEEE 118-bus system and comparing the efficiency with serial computing, we can obtain the result of calculation time, as shown in the fourth to the sixth columns of Table 3.It shows that the time-consuming of SQP-AGS using parallel computing is much less than the one under the serial computing process, while the iterations of both of them are the same. The time-consuming is reduced with the number of increasing workers as well. On this basis, we calculate the speedup ratio and the efficiency of parallel computing with different numbers of workers running. The speedup ratio of parallel computing is defined asSn=T1/Tn, in which n is the number of parallel workers and T1 is the computing time with only one worker. The efficiency of parallel computing is defined as En=Sn/n. Table 3 separately lists the speedup ratios of 2, 4, and 8 workers, respectively, which are S2=1.71, S4=2.86, and S8=4.05. The efficiencies are E2=0.86, E4=0.72, and E8=0.51, correspondingly. Due to the characteristic of sampled points, multicore parallel computing plays an important role in shortening the calculation time. Compared with the SQP-GS which is used in [5], the SQP-AGS with the multicore parallel computing technique greatly improves the computing efficiency of the SSSC-ATC problem of a large-scale power system, which has absolute advantages and more practical application prospects in practical engineering applications.


## 4.1. Comparison of Results under Different Small-Signal Stability Constraints Run from IEEE 10-Machine 39-Bus System

The IEEE 10-machine 39-bus system is used for stability analysis with different small-signal stability constraints, which are three different spectral abscissa constrained values in the SSSC-ATC calculation model above. The specific regional division is shown in Figure2. The generators G2, G3, and G10 which, respectively, correspond to bus 31, 32, and 39 are selected to locate at the receiving side, while the other generators are located at the source side.Figure 2
The regional division of the IEEE 39-bus system.For the voltage magnitude limits, we chooseV¯=1.1 p.u. and V¯=0.9 p.u. for all buses. The results under different spectral abscissa constraints are listed in Table 1. In order to reflect the effect on the SQP-AGS algorithm applied on SSSC-ATC, we set the sample size as p=30 and make a comparison of the following cases:(1)
Case 0: base case, which is a standard ATC calculation without any small-signal stability constraint(2)
Case 1: SSSC-ATC with an abscissa spectrum constraint ofη≤−0.10(3)
Case 2: SSSC-ATC with an abscissa spectrum constraint ofη≤−0.192Table 1
Active output of generators and TTC in IEEE 39-bus system.
Active output (MW)
Active output range (MW)
Case 0
Case 1
Case 2
G1
590.13
1102.59
1106.74
100∼1380
G2
50.00
105.35
136.58
50∼747.5
G3
50.00
677.18
676.12
50∼920.0
G4
862.50
665.66
639.15
50∼862.5
G5
747.50
529.28
564.18
50∼747.5
G6
862.50
671.67
667.62
50∼862.5
G7
862.50
581.69
603.44
50∼862.5
G8
805.00
597.35
513.94
100∼805.0
G9
1035.00
899.08
898.57
100∼1035
G10
402.50
318.83
346.63
50∼402.5
TTC (MW)
1704.99
498.62
465.73
—(i)
Case 0: the classic IPM is used to solve the ATC problem without any small-signal stability constraint. The results are listed in the second column of Table1, which shows the TTC is obtained as 1704.99 MW and the spectral abscissa is calculated to be 2.15, indicating the TTC is overoptimistic due to the fact that the system is not small-signal stable in this operation state. This further indicates that the influence of small-signal stability on ATC calculation needs to be considered. Extracting information from Table 1, it can be seen that the active generation of generators G2 and G3 which are located in the receiving side all reach the lower limit of the active power output in order to ensure the maximum of TTC.
(ii)
Case 1: a spectral abscissa constraintη≤−0.10 is added to the ATC model as the small-signal stability constraint and SQP-AGS algorithm is run to solve it. The results are listed in the third column of Table 1, which shows the TTC is obtained as 498.62 MW, and the system is small-signal stable in this operation state. The result also shows that the generation in the receiving side is greatly increased compared with Case 0, especially as generator G3, possibly due to the larger time inertia coefficient of it, resulting in the reduction of TTC. It reflects that the small-signal stability does constrain the regional transmission directly.
(iii)
Case 2: the spectral abscissa constraint is further reduced toη≤−0.192 in the SSSC-ATC model. The results run from the SQP-AGS algorithm shows that the calculation is still convergent and the TTC is reduced to 465.73 MW.Comparing Case 0, Case 1, and Case 2 will find out a rule that TTC value reduces as the spectral abscissa decreases, which means the stricter restriction effect on ATC calculation when a spectral abscissa constraint gets smaller. The results under spectral abscissa constraints are more reasonable than the one without any small-signal stability constraint so that the ATC in this situation can be a guide to avoid the power oscillation caused by small-signal stability in the system.


## 4.2. Comparison of Different Methods in Calculating SSSC-ATC

Here, we compare the SQP-AGS’s calculation effect on the same problem with a sensitivity method [3] and a participation factor method [4] for IEEE 39-bus system. The results run by the three methods are all restricted to a spectral abscissa constraint as η≤−0.1 and detailedly listed in Table 2.Table 2
Comparison of different methods for the 39-bus system.
SQP-AGS (MW)
Sensitivity method (MW)
Participation factor method (MW)
G1
1102.59
1210.27
1357.84
G2
105.35
556.43
150.46
G3
677.18
407.27
704.98
G4
665.66
732.91
680.85
G5
529.28
599.08
673.28
G6
671.67
416.40
862.50
G7
581.69
764.68
50.00
G8
597.35
715.46
733.05
G9
899.08
790.18
627.21
G10
318.83
50.00
402.50
TTC (MW)
498.60
289.94
262.78The TTC calculated by SQP-AGS is 498.60 MW, while the TTC calculated by the sensitivity method and the participation factor method is 289.94 MW, and the one calculated by the participation factor method is 262.78 MW. All the three methods obtain a TTC value much smaller than which is calculated from the ATC model without small-signal stability constraint. In addition, compared with the result calculated by SQP-AGS, the other two methods’ results are relatively conservative. In other words, it is to say that SQP-AGS can evaluate the SSSC-ATC problem in an optimization way. It further means this ATC result may bring more economic benefits that are more acceptable to the power industry in practice.On the other hand, it is worth mentioning that when the spectral abscissa constraint decreases asη≤−0.15, the sensitivity method and the participation method both cannot converge to obtain a result, while the SQP-AGS can still effectively reduce the spectral abscissa of the system and achieve a good convergent result, as is shown in Table 1. This reveals the superiority of SQP-AGS in solving the SSSC-ATC problem compared with the other two mentioned methods.


## 4.3. Efficiency Tested in IEEE 54-Machine 118-Bus System

To test the efficiency of SQP-AGS, we make some comparisons on a modified version of the IEEE 54-machine 118-bus benchmark system which has 54 synchronous machines with IEEE type-1 exciters, 20 of which are synchronous compensators used only for reactive power support and 15 of which are motors. The specific regional division of IEEE 118-bus system is as follows:(i)
Receiving side: the region includes buses of 1, 4, 6, 8, 10, 12, 15, 18, 19, 24, 25, 26, 27, 31, 32, 34, 36, 40, 42, 46, 49, 54, 55, 56, 59, 61, 62, 65, 66, 69, 70, 72, 73, 113, and 116.(ii)
Source side: the region includes buses of 74, 76, 77, 80, 85, 87, 89, 90, 91, 92, 99, 100, 103, 104, 105, 107, 110, 111, and 112.Lines 74–70, 75–70, 75–69, 77–69, and 81–68 are the tie lines between the two regions. For an ATC model without any small-signal stability constraint, the TTC is obtained as 1410.20 MW and the spectral abscissa isη=1.1, which also means the calculated ATC is overoptimistic, and power oscillation will occur before the transmission between regions reaches this value. Applying SQP-AGS to solve an SSSC-ATC problem with a spectral abscissa constraint η≤−0.1, the spectral abscissa is decreased meeting the constraint requirement after 17 iterations and the result is convergent, obtaining the TTC as 446.08 MW. The spectral abscissa changing with iterations is shown in Figure 3. It is easy to find that the descending trend of spectral abscissa is stable. Moreover, Figure 4 reflects the changes in the termination conditions with the iteration increases. As the trend shows, the infeasibility defined in (28) and the model reduction calculated by (30) both decrease to an acceptable tolerance.Figure 3
Spectral abscissa changes with iterations for the 118-bus system.Figure 4
The termination conditions change with iterations for the 118-bus system.In order to test the efficiency of SQP-AGS, the algorithm is going with a sample size ofp=50. Under the serial computing condition, we make a comparison between SQP-GS and SQP-AGS dealing with the same SSSC-ATC problem. The computing time of them is listed in second and third columns of Table 3.Table 3
Efficiency comparison of the IEEE 118-bus system.
Serial computing
Parallel computing for SQP-AGS
SQP-GS method
SQP-AGS method
Two cores
Four cores
Eight cores
Total time (s)
173.48
66.29
39.88
25.15
16.36
Gradient sampling time (s)
170.40
63.32
36.90
20.24
13.37
QP solver time (s)
2.78
2.67
2.68
2.62
2.69
Others’ time (s)
0.30
0.30
0.30
0.29
0.30
Speedup ratio,Sn
—
—
1.71
2.86
4.05
Efficiency,En
—
—
0.86
0.72
0.51
Iteration
18
17
17
17
17It is clear to see that, for the 118-bus system, applying the serial SQP-GS method to solve an SSSC-ATC problem consumes 173.48s and 18 iterations while applying the serial SQP-AGS needs only 66.29s and 17 iterations, cutting down a lot of time in solving the same model. Particularly, the time-consuming of gradient calculating accounts for the vast majority of the total time because of the repeating calculating processes of small-signal stability calculation and the sampling gradient calculation of nonsmooth constraint. Compared with SQP-GS, SQP-AGS spends much less time on gradient calculation due to the adaptive identification of sampling points which are within the sampling neighborhood of the current iteration point. In each iteration, SQP-AGS effectively utilizes the historical sampling data in the previous iteration and saves a lot of time for gradient calculation at current iteration, which greatly improves the efficiency of calculation. The intuitive sample point changes are shown in Figure5.Figure 5
Number change of new sampling points.Moreover, notice that the sample points for gradient calculation in each iteration are random and independent, we adopt the multicore parallel computing technology of MATLAB to carry out the gradient calculation when applying SQP-AGS. Parallel computing opens the parallel computing pool of specified number of workers through theparpool instruction in MATLAB and uses parfor statement to execute the operation process of the loop body. Testing it on the IEEE 118-bus system and comparing the efficiency with serial computing, we can obtain the result of calculation time, as shown in the fourth to the sixth columns of Table 3.It shows that the time-consuming of SQP-AGS using parallel computing is much less than the one under the serial computing process, while the iterations of both of them are the same. The time-consuming is reduced with the number of increasing workers as well. On this basis, we calculate the speedup ratio and the efficiency of parallel computing with different numbers of workers running. The speedup ratio of parallel computing is defined asSn=T1/Tn, in which n is the number of parallel workers and T1 is the computing time with only one worker. The efficiency of parallel computing is defined as En=Sn/n. Table 3 separately lists the speedup ratios of 2, 4, and 8 workers, respectively, which are S2=1.71, S4=2.86, and S8=4.05. The efficiencies are E2=0.86, E4=0.72, and E8=0.51, correspondingly. Due to the characteristic of sampled points, multicore parallel computing plays an important role in shortening the calculation time. Compared with the SQP-GS which is used in [5], the SQP-AGS with the multicore parallel computing technique greatly improves the computing efficiency of the SSSC-ATC problem of a large-scale power system, which has absolute advantages and more practical application prospects in practical engineering applications.


## 5. Conclusion

In this paper, we propose an SSSC-ATC calculating model with a spectral abscissa constraint which directly and fully considers the impact of small-signal stability on ATC calculation. An SQP-AGS method is adopted to solve the proposed SSSC-ATC problem, which is a nonsmooth optimization problem due to the property of the spectral abscissa function. The closed-form eigenvalue sensitivity is used to calculate the gradient of the spectral abscissa. In SQP-AGS, an AGS theory is adopted to evaluate the gradients around the current iteration points adaptively to make the search direction computation effective in the nonsmooth regions and practically improve the efficiency of calculation. The multicore computing technology further improves the effect on calculation efficiency. Simulation results on two test systems show that the proposed method solves the SSSC-ATC problem effectively without any convergence problem. Comparison of other two methods proves that the proposed method is more effective in calculating the SSSC-ATC problem.



---

*Source: 3912717-2020-02-11.xml*