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.
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.
Section 2 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 aswhere 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: where is the set of tie lines between the source side and the receiving side. indicates the active power of the tie line , which is defined as . represents the voltage amplitude of the ith bus. and are the admittance amplitude and the phase angle of line , respectively. , in which is the phase angle of the ith bus.(2)For all buses in system denoted as collection , the power flow equations for bus is where is the active power output of the ith generator, is the reactive power output of the ith generator correspondingly, and and are the active power load and the reactive power load of the ith bus, respectively.(3)Technical constraints include where is the collection of all the generator buses and is the collection of all lines. represents the apparent power of line , 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 is the state variable vector including the rotor angle vector , the d-axis and q-axis component vectors of the internal voltage are and , the d-axis and q-axis component vectors of the internal current are and , and the excitation output voltage vector is .(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.
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 of , 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 of , 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 formdefined in an open, dense subset , the objective function f is locally Lipschitz continuous which is nonsmooth and/or nonconvex.
As Figure 1 shows, the gradients of smooth points like and are easily obtained and unique. But it is unable to obtain the gradient at a nonsmooth point . Moreover, the subgradient of nonsmooth point exists and is not unique. According to the definition, denotes the subgradient at the nonsmooth point . The set of subgradients is called the subdifferential of f at .
The GS technique approximates the subgradients of the objective function through the gradients sampled randomly near aswhere P is the amount of the sampled gradients. Defining that is a set of sampled gradients evaluating from a set of random and independent sampling points near , whereis the closed ball with a sample radius ϵ centered at . The function f is locally Lipschitz and almost continuously differentiable.
The expression of 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 in is assigned a nonnegative weight or coefficient that all of them are summed to one. In such a way, the set of subgradients of is derived by changing 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 is Clarke stationary if
Figure 1 marks out a zero-valued subgradient determining that 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 subgradient at the kth iteration, a following quadratic programming (QP) model is established:where 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 can be obtained. Then, the subgradient can be calculated according to (9). The vector is defined as an approximate steepest descent direction. If , 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:
The SQP-GS algorithm solves the model above in each iteration mainly by determining a search direction centering on solving a regularized QP dual subproblem with the following form:where ρ is a penalty parameter and is the approximate Hessian of the Lagrangian optimization model. Incorporating the ideas from the GS algorithm, , , and are sets of independent and identically distributed random points uniformly sampled from (10) with a sample size of p:
Consistent with the GS strategy, , , and 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 .
As a significant step to update during each iteration, we introduce the limited-memory Broyden–Fletcher–Goldfrab–Shanno (L-BFGS) method which is an attempt to achieve fast convergence. We defineas the approximation of the minimum-norm element of the differential of (13) Lagrange function, relating to the choice of made during each iteration. We initialize at the 1th iteration and perform the update asfor . Here, and are the displacements in 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 is an infeasibility vector at the kth iteration defined aswhich is also an indicator in the step of parameter update.
Furthermore, let us define the local model
For a computed , the model reduction is defined to be , which values
The model reduction is nonnegative because yields an objective value, and can be zero only if is ϵ-stationary. So a sufficiently small reduction in 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 points . Initialize a sampling size value p at first, which is the total number of sample points in every iteration. For (with for some i and , AGS setswhere the sample set is composed of points generated uniformly in . Denote the number of points selected by is . The actual sampling size of kth iteration changes to , equaling to the number of points in .
It is apparent to see that AGS needs only to compute the unknown gradients at the updating sample points in 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 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.
|
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 power and reactive power , the sample radiuses and are stipulated as of the capacity range. For voltage-type variables (per-unit value) like and , the sample radiuses and are chosen as 0.03 p.u. For angle variables like θ, the sample radius is chosen as .
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 function quadprog 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., and . The parameters of SQP-AGS are chosen as and . The selection of sample radius ϵ depends on the characteristics of the dependent variables as described in Subsection 3.3. Set , , , , and from [9]. The tolerances and are set to be and , 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 Figure 2. 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.
For the voltage magnitude limits, we choose p.u. and 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 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 (3)Case 2: SSSC-ATC with an abscissa spectrum constraint of
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 and detailedly listed in Table 2.
The 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 , 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 , 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 , 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.
In order to test the efficiency of SQP-AGS, the algorithm is going with a sample size of . 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.
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 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 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. 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 as , in which n is the number of parallel workers and is the computing time with only one worker. The efficiency of parallel computing is defined as . Table 3 separately lists the speedup ratios of 2, 4, and 8 workers, respectively, which are , , and . The efficiencies are , , and , 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.
Data Availability
The 10-machine 39-bus system and 54-machine 118-bus system data that support the findings of this study are available in the thesis [11] and at http://www.kios.ucy.ac.cy/testsystems/index.php/dynamicieee-test-systems/ieee-118-bus-modified-test-system [12].
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported in part by the National Natural Science Foundation of China under Grant nos. 51967001 and 51667003, in part by Natural Science Foundation of Guangxi under Grant no. 2018GXNSFAA281194, and in part by Guangxi Key Laboratory of Power System Optimization and Energy Technology Research Grant.