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

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.)e highlighted feature is the sample size of the gradient changes dynamically in every iteration, yielding an adaptive gradient sampling (AGS) process. )us, 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 54machine 118-bus system prove the effectiveness and high efficiency of the proposed method.


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, smallsignal 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. erefore, 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. e 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). e 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. e 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. us, 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.
irdly, 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.
Section 2 describes a model of the ATC calculation with a small-signal stability constraint. e 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.

The Model of ATC Calculation with a Small-Signal Stability Constraint
Mathematically, ATC is defined as 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. us, 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) e objective value: where S link is the set of tie lines between the source side and the receiving side. P ij indicates the active power of the tie line (i, j), which is defined as represents the voltage amplitude of the ith bus. Y ij 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 collection S B , the power flow equations for bus i ∈ S B is where P Gi is the active power output of the ith generator, Q Gi is the reactive power output of the ith generator correspondingly, and P Li and Q Li are the active power load and the reactive power load of the ith bus, respectively. (3) Technical constraints include where S G is the collection of all the generator buses and S line is the collection of all lines. S ij 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: where x is the state variable vector including the rotor angle vector δ, the d-axis and q-axis component vectors of the internal voltage are E d ′ and E q ′ , the d-axis and q-axis component vectors of the internal current are I d ′ and I q ′ , and the excitation output voltage vector is E fd . (5) Small-signal stability constraint: 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. e spectral abscissa η in small-signal stability constraint is a nonsmooth function [7]. us, 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 R n , which means that it is continuously almost everywhere and its gradient can be easily obtained where it is defined by calculating spectral abscissa derivative.

A Sequential Quadratic Programming
Algorithm Combined with Gradient Sampling ere 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 of R n , 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.
e GS is basically a stabilized steepest descent algorithm. Considering an unconstrained optimization problem of the form defined in an open, dense subset D ∈ R n , the objective function f is locally Lipschitz continuous which is nonsmooth and/or nonconvex. As Figure 1 shows, the gradients of smooth points like x 1 and x 2 are easily obtained and unique. But it is unable to obtain the gradient at a nonsmooth point x 0 . Moreover, the subgradient of nonsmooth point exists and is not unique. According to the definition, u ∈ zf denotes the subgradient at the nonsmooth point x 0 . e set of subgradients u is called the subdifferential of f at x 0 . e GS technique approximates the subgradients of the objective function through the gradients sampled randomly near x 0 as where P is the amount of the sampled gradients. Defining that G k ≔ ∇f(x k0 ), ∇f(x k1 ), . . . , ∇f(x kp ) is a set of sampled gradients evaluating from a set of random and independent sampling points X k ≔ x k0 , x k1 , . . . , is the closed ball with a sample radius ϵ centered at x 0 . e function f is locally Lipschitz and almost continuously differentiable. e expression of zf(x 0 ) 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 ∇f(x ki ) in G k is assigned a nonnegative weight or coefficient y ki that all of them are summed to one. In such a way, the set of subgradients of x 0 is derived by changing y ki into different values.
It is known that the conventional way to find the stationary point of a smooth function f 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 x 0 is Clarke stationary if 0 ∈ zf(x).
(11) Figure 1 marks out a zero-valued subgradient u 1 determining that x 0 is a Clarke stationary point. e 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 subgradient u k at the kth iteration, a following quadratic programming (QP) model is established: where y k is the coefficient vector at the kth iteration. e 2norm in the objective function (12) represents the shortest distance from the set of subgradients to the origin. Solving this model, the coefficient vector y k can be obtained. en, the subgradient u k can be calculated according to (9). e vector d k � − u k /‖u k ‖ is defined as an approximate steepest descent direction. If ‖u k ‖ � 0, the stationary point is found.

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) e SQP-GS algorithm solves the model above in each iteration mainly by determining a search direction d k centering on solving a regularized QP dual subproblem with the following form: where ρ is a penalty parameter and H k is the approximate Hessian of the Lagrangian optimization model. Incorporating the ideas from the GS algorithm, are sets of independent and identically distributed random points uniformly sampled from (10) with a sample size of p: Consistent with the GS strategy, y f � (y As a significant step to update H k 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 as the approximation of the minimum-norm element of the differential of (13) Lagrange function, relating to the choice of H k made during each iteration. We initialize H 1 ⟵I at the 1th iteration and perform the update as for l � 2, . . . , k. Here, s l ≔ x l − x l− 1 and y l ≔ G l − G l− 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: where σ(x) is an infeasibility vector at the kth iteration defined as which is also an indicator in the step of parameter update. Furthermore, let us define the local model For a computed d k , the model reduction is defined to be △q ρ (d k ; e model reduction is nonnegative because d k � 0 yields an objective value, and △q ρ (d k ; x k , B ϵ,k , H k ) can be zero only if x k is ϵ-stationary. So a sufficiently small 4 Complexity reduction in △q ρ (·; x k , B ϵ,k , H k ) 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. e difference between GS and AGS is mainly in the process of sample points x k,p ∈ B ϵ,k . Initialize a sampling size value p at first, which is the total number of sample points in every iteration. For X k ≔ x k,1 , . . . , x k,p (with where the sample set X k ≔ x k,1 , . . . , x k,p is composed of p points generated uniformly in B k . Denote the number of points selected by (X k− 1 ∩ B k ) is p re . e actual sampling size of kth iteration changes to p � p − p re , 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 in X k 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. e complete SQP-AGS algorithm is presented in Algorithm 1.
What is worth mentioning is that the sample size p 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.

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 smallsignal stability constraint for it is the only the nonsmooth constraint in the proposed SSSC-ATC model. erefore, 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 power P G and reactive power Q G , the sample radiuses ϵ P G and ϵ Q G 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 ∘ .
irdly, 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. e application of parallel computing technology will further improve the efficiency of SQP-AGS.

Case Studies
e proposed SQP-AGS method is applied to the IEEE 10machine 39-bus system and the modified IEEE 54-machine 118-bus system to illustrate the effectiveness in solving SSSC-ATC. e 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. e 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.
e SQP-AGS method is implemented in MATLAB by using the function quadprog as the QP solver for the subproblem. e 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., P G � (P G + P G )/2 and Q G � (Q G + Q G )/2. e parameters of SQP-AGS are chosen as ρ � 2 × 10 − 4 and μ ρ � 0.05. e 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 c � 0.8 from [9].

Comparison of Results under Different Small-Signal Stability Constraints Run from IEEE 10-Machine 39-Bus System.
e 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. e specific regional division is shown in Figure 2. e 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.
For the voltage magnitude limits, we choose V � 1.1 p.u. and V � 0.9 p.u. for all buses. e 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.192 (i) Case 0: the classic IPM is used to solve the ATC problem without any small-signal stability Complexity 5 constraint.
e results are listed in the second column of Table 1, 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. is 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. e 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. e 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. e 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. e 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.

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. e results run by the three methods are all restricted to a spectral abscissa constraint as η ≤ − 0.1 and detailedly listed in Table 2. e 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 smallsignal 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. is reveals the superiority of SQP-AGS in solving the SSSC-ATC problem compared with the other two mentioned methods.

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 118bus 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. e specific regional division of IEEE 118bus system is as follows: (1) Initialization: choose a sample radius ϵ > 0, penalty parameter ρ > 0, constraint violation tolerance τ > 0, sample size p, line search constant ϖ ∈ (0, 1), backtracking constant c ∈ (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 x 0 ∈ D and set k � 1 and k max � 200.  (d k , y o , y h , y g ). (5) L-BFGS update: update H k by (26) for next iteration. (6) Parameter update: if d k > ] s ϵ 2 , then go to step 9. Otherwise, if max(σ k ) ≤ τ, set τ ⟵ μ τ τ, but if max(σ k ) > τ, set ρ ⟵ μ ρ ρ. en, set ϵ ⟵ μ ϵ ϵ, β k ⟵ 0, and go to step 8. (7) Line search: set β k as the largest value in the sequence 1, c, c 2 , . . . such that  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. e 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.
In order to test the efficiency of SQP-AGS, the algorithm is going with a sample size of p � 50. Under the serial computing condition, we make a comparison between SQP-GS and SQP-AGS dealing with the same SSSC-ATC problem. e computing time of them is listed in second and third columns of Table 3.
It 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  Complexity 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. e intuitive sample point changes are shown in Figure 5.
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 the parpool 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. e 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. e speedup ratio of parallel computing is defined as S n � T 1 /T n , in which n is the number of parallel workers and T 1 is the computing time with only one worker. e efficiency of parallel computing is defined as E n � S n /n. Table 3 separately lists the speedup ratios of 2, 4, and 8 workers, respectively, which are S 2 � 1.71, S 4 � 2.86, and S 8 � 4.05. e efficiencies are E 2 � 0.86, E 4 � 0.72, and E 8 � 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  8 Complexity 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.

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. e 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. e 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.

Conflicts of Interest
e authors declare that they have no conflicts of interest.