An Asymptotic Numerical Continuation Power Flow to Cope with Non-Smooth Issue

: Continuation power ﬂow (CPF) calculation is very important for analyzing voltage stability of power system. CPF calculation needs to deal with non-smooth constraints such as the generator buses reactive power limits. It is still a technical challenge to determine the step size while dealing with above non-smooth constraints in CPF calculation. In this paper, an asymptotic numerical method (ANM) based on Fischer-Burmeister (FB) function, is proposed to calculate CPF. We ﬁrst used complementarity constraints to cope with non-smooth issues and introduced the FB function to formulate the complementarity constraints. Meanwhile, we introduced new variables for substitution to meet the quadratic function requirements of ANM. Compared with the conventional predictor-corrector method combining with heuristic PV-PQ (PV and PQ are used to describe bus types. PV means that the active power and voltage of the bus are known. PQ means that the active and reactive power of bus are known.) bus type switching, ANM can e ﬀ ectively solve the PV-PQ bus type switching problem in CPF calculation. Furthermore, to assure high e ﬃ ciency, ANM can rapidly approach the voltage collapse point by self-adaptive step size adjustment and constant Jacobian matrix used for power series expansion. However, conventional CPF needs proper step set in advance and calculates Jacobian matrix for each iteration. Numerical tests on a nine-bus network and a 182-bus network validate that the proposed method is more robust than existing methods.


Introduction
Power flow calculation is one of the most basic calculation in power system analysis, and it is also the basis of power system stability calculation and fault analysis [1]. Continuation power flow (CPF) calculation combines the continuous method with the static power flow of power system, which is widely used in static voltage stability analysis of power system [2][3][4][5][6][7][8]. However, it is still a technical challenge to deal with the non-smooth constraints such as the reactive power limits violation in CPF calculation [9].
The reactive power limit in CPF calculation has important effects on the voltage stability limit of power system. The traditional method of solving reactive power limit is the logical switch of PV-PQ (PV and PQ are used to describe bus types. PV means that the active power and voltage of the bus are known. PQ means that the active and reactive power of bus are known.) bus type. In [10], a new definition of bus types was presented. The solvability criteria and solution method of bus-type extended power flow are given. And the double switching logic of the new bus types is given to handle the reactive power limits of generators. It is complicated to deal with PV-PQ switch by traditional logic switch method. In [11], a CPF method which considers the reactive limits of generators as a part of the algorithmic procedure was proposed. Lagrange polynomial interpolation formula is used to

Semi-Smooth Quadratic Power Flow Equations
The power flow problem is formulated as a set of nonlinear mismatch equations for active and reactive power, as follows, For a network with N buses, active and reactive power satisfy that: (G ij f j + B ij e j ) = 0, where P i , Q i are active power and reactive power of bus i, e i , f i are real and imaginary part of bus i voltage, G ij , B ij are conductance and susceptance matrix elements respectively. For PV buses, Equation (2) can be replaced by the equation of voltage: where U set i is voltage target of PV buses. For voltage control buses PV and slack bus, when the reactive power exceeds the lower limit, the bus reactive power is set to the lower limit. When the reactive power exceeds the upper limit, the bus reactive power is set to the upper limit. When the bus reactive power does not exceed the limit maintain the original value.
Describe the above relationship with complementary constraints as follows: where Q max i and Q min i represent the maximum and minimum limits of reactive power of bus i respectively, U + i and U − i are slack variables for voltage regulation. There is a general rule of power system. When the reactive power of power system is greater than the maximum limit, the voltage regulation target will be lowered by −U − i , when the reactive power of power system is less than the minimum limit, the voltage regulation target will be enlarged by +U + i . To cope with the complementarity constraints u⊥v, the Fischer-Burmesiter (FB) function [19], is introduced: where a slack variable µ = 10 −20 is introduced to avoid non-differentiable problem of FB at (u = 0, v = 0). In order to turn Equation (5) into a quadratic equation to satisfy the calculation requirements of ANM, we defined a new variable again for variable substitution, the Equation (5) is converted as follows: where w = u 2 + v 2 + µ. FB function has nice properties, such as strong semi-smoothness. Moreover, the squared norm of FB function has a Lipschitz continuous gradient. The FB functions for the complementarity constraints in the power flow can be arranged as: where w + and w − are intermediate variables.
Finally, power mismatch Equations (1)-(3) and (7) constitute the power flow equations, which can be solved by using ANM.

Algorithm of ANM-CPF
The ANM algorithm is a continuous method for high-order prediction using power series expansion [20]. The solution curve of the nonlinear equation with parameters is segmented, and the solution curve of each segment can be analytically expressed in the form of closed power series. The nonlinear problem is transformed into an infinite amount of linear subproblems, and the sum of the solutions of the first few nonlinear subproblems is used to approximate the real solution. Since the predicted solution is almost the same as the real solution, the ANM-CPF continuous process does not require a correction step in general, and the calculation step size can also be adaptively adjusted. In order to clearly explain the calculation principle of the ANM algorithm, a two-bus example is introduced in the Appendix A. The detailed calculation process and results are listed there.
Further decompose the above Equations (1)-(3) and Equation (7) by a homotopy transformation: where x represents variables, including parameters such as the real and imaginary part of bus voltage e i , f i , reactive power Q i , and intermediate variables w + and w − . λ is the scalar embedding parameter used to describe the increase of active power P i and reactive power Q i . L and Q represent linear and bilinear operators respectively in Equations (1)-(3) and (7). F is the increase direction of active and reactive power at each bus. H is a constant vector consisting of the reactive power limits Q min i , Q max i , slack variable µ, etc. The detailed decomposition process of above power flow equations for the two-bus network corresponds to (A1) and (A4)-(A7) in Appendix A.
Suppose that (x t , λ t ) is the t th (t = 0, 1, 2 . . .) computed point, the q th point between t th and (t + 1) th step can be expressed as the series like: where K is the truncation order of the series, (∆s) q is the step size of q th point between t th and (t is the Taylor coefficients. The maximum step size is divided into M equal parts. And the step size of q th point between t th and (t + 1) th satisfies: where ∆s max is the maximum step size of the t th computed point and when q is equal to M, x q t = x t+1 . The relationship between step size ∆s and the maximum step size ∆s max for a practical two-bus network can be seen in the Table A1 in Appendix A.
The Taylor coefficients (x (p) t , λ (p) t ) at different orders are determined by a series of equations: Order p ≥ 2: where f t x is the value of Jacobian matrix at x t , and v 1 , v p are intermediate variables. Constant Jacobian matrix at x t is used for each point between x t and x t+1 when using power series expansion. The calculation of Jacobian matrix and Taylor coefficient in a practical two-bus network corresponds to (A8) in Appendix A.
) is a matrix composed of functions that construct quadratic terms of the equations: where A h Q (1 ≤ h ≤ n) is a sparse coefficient matrix, the value of n is equal to the dimension of variable x. The maximum step size ∆s max satisfies: where ε is the calculation accuracy control parameter. In view of (A9) and (A10) in Appendix A, the calculation process of the maximum step size ∆s max can be seen clearly. Given (12) and (14), ∆s max can be rewritten as follows:

Case Study
In this section, case studies are conducted on a small IEEE nine-bus network and a large 182-bus network. The detailed data is referred to [21,22]. The comparison between the ANM and the conventional predictor-corrector method is as follows. Both methods start from the same initial points.

Without the Consideration of Reactive Power Limit in CPF
In order to verify the correctness of the ANM, we set the limit of the system reactive power constraint to infinity. This ensures the PV buses or slack bus do not violate the limit in CPF, so the reactive power constraints do not work, as shown in Figure 1. When the scalar embedding parameter λ increases to λ = 1.641, the system load power increases to 2.641 times the initial value. The voltage stability limit will be reached, causing the system voltage to collapse.
where ε is the calculation accuracy control parameter. In view of (A9) and (A10) in Appendix A, the calculation process of the maximum step size max s Δ can be seen clearly.

Case Study
In this section, case studies are conducted on a small IEEE nine-bus network and a large 182-bus network. The detailed data is referred to [21,22]. The comparison between the ANM and the conventional predictor-corrector method is as follows. Both methods start from the same initial points.

Without the Consideration of Reactive Power Limit in CPF
In order to verify the correctness of the ANM, we set the limit of the system reactive power constraint to infinity. This ensures the PV buses or slack bus do not violate the limit in CPF, so the reactive power constraints do not work, as shown in Figure 1. When the scalar embedding parameter λ increases to 1.641 λ = , the system load power increases to 2.641 times the initial value. The voltage stability limit will be reached, causing the system voltage to collapse. From Figure 1, we can see that the solutions by using ANM is completely consistent with those by using the conventional predictor-corrector method in PF iteration. It shows that ANM inherits the excellent characteristics of high precision, verifying the correctness of the model and algorithm. In addition, it is easy to conclude that the step size can be adaptively adjusted to improve the calculation   From Figure 1, we can see that the solutions by using ANM is completely consistent with those by using the conventional predictor-corrector method in PF iteration. It shows that ANM inherits the excellent characteristics of high precision, verifying the correctness of the model and algorithm. In addition, it is easy to conclude that the step size can be adaptively adjusted to improve the calculation speed in PF iteration by using ANM.
At present, CPF is the most commonly used voltage stability analysis method of large power systems [23,24]. In order to validate the results of the IEEE nine-bus network, we carried out transient stability studies in DlgSILENT. The voltage of slack bus 1, PV buses 2 and 3 remain constant until the scalar embedding parameter λ increases to λ = 1.641. When time = 100 s, increase the load power to 2.641 times the initial value, we can see that the voltage of bus 1, 2, and 3 constantly oscillate during the subsequent iterative process, as shown in Figure 2. It proves that the system voltage will be unstable when it reaches the voltage stability limit, which is consistent with the calculation of CPF. . When time = 100 s, increase the load power to 2.641 times the initial value, we can see that the voltage of bus 1, 2, and 3 constantly oscillate during the subsequent iterative process, as shown in Figure 2. It proves that the system voltage will be unstable when it reaches the voltage stability limit, which is consistent with the calculation of CPF. In view of (15), ( can be understood as the vertical increment of function ( , ) t t x λ . When the PV curve gently changes, the vertical increment is small, the maximum step size is large, whereas, when the curve steeply changes, the vertical increment is large, the maximum step size is small. So the step size adaptive adjustment can be realized, as shown in Figure 3. The maximum step size at different points on IEEE nine-bus network can be seen in Figure 4. From Figure  4, we can see that the maximum step size near the voltage bifurcation point is smaller than other points. It can be further verified that the step size of ANM can be adjusted adaptively.  In view of (15), (− f t x v K ) can be understood as the vertical increment of function f (x, λ) at point (x t , λ t ). When the PV curve gently changes, the vertical increment is small, the maximum step size is large, whereas, when the curve steeply changes, the vertical increment is large, the maximum step size is small. So the step size adaptive adjustment can be realized, as shown in Figure 3. The maximum step size at different points on IEEE nine-bus network can be seen in Figure 4. From Figure 4, we can see that the maximum step size near the voltage bifurcation point is smaller than other points. It can be further verified that the step size of ANM can be adjusted adaptively.

Consideration of Reactive Power Limit in CPF
Considering reactive power limit in CPF, the reactive power limit data of IEEE nine-bus network can be referred to [21]. When scalar embedding parameter λ increases to 1.533 λ = , the system load power increases to 2.533 times the initial value. The reactive power of bus 1 will violate its reactive power limit, causing the system voltage to collapse, as shown in Figure 5.

Consideration of Reactive Power Limit in CPF
Considering reactive power limit in CPF, the reactive power limit data of IEEE nine-bus network can be referred to [21]. When scalar embedding parameter  increases to 1.533  = , the system load power increases to 2.533 times the initial value. The reactive power of bus 1 will violate its reactive power limit, causing the system voltage to collapse, as shown in Figure 5.

Consideration of Reactive Power Limit in CPF
Considering reactive power limit in CPF, the reactive power limit data of IEEE nine-bus network can be referred to [21]. When scalar embedding parameter λ increases to λ = 1.533, the system load power increases to 2.533 times the initial value. The reactive power of bus 1 will violate its reactive power limit, causing the system voltage to collapse, as shown in Figure 5.
The comparison on the central processing unit (CPU) time for both methods is listed in Table 1. In view of Equations (11) and (12), ANM can reduce the time of calculating Taylor coefficient by using constant Jacobian matrix at each point between x t and x t+1 when using power series expansion. The constant Jacobian matrix can lessen calculative burden on the premise of ensuring the accuracy of calculation. However, conventional CPF needs to calculate the Jacobian matrix for each iteration. It can be seen that ANM-CPF using constant Jacobian matrix for power series expansion assure the proposed method more efficient. The comparison on the central processing unit (CPU) time for both methods is listed in Table 1. In view of Equations (11) and (12), ANM can reduce the time of calculating Taylor coefficient by using constant Jacobian matrix at each point between t x and 1 t x + when using power series expansion.
The constant Jacobian matrix can lessen calculative burden on the premise of ensuring the accuracy of calculation. However, conventional CPF needs to calculate the Jacobian matrix for each iteration. It can be seen that ANM-CPF using constant Jacobian matrix for power series expansion assure the proposed method more efficient. When the reactive power of some PV buses and slack bus violate their limits, they are converted to PQ buses or Qθ (Qθ is used to describe bus types, which means that the reactive power and phase angle of the bus are known) bus respectively. In the nine-bus system, the conventional predictor-corrector method cannot smoothly realize the PV-PQ bus switching. The step size is 0.01 s = . The PV curves of bus 9 for both methods are as follows. The points of PV curve are sparse near the voltage bifurcation point by using the conventional predictor-corrector method, which is not suitable for obtaining accurate maximum load point. But the ANM can smoothly realize the PV-PQ bus switching and the step size is automatically small near the voltage bifurcation point because of the self-adaptive adjustment. The PV curve can be smoothly drawn, as shown in Figure 6. Correspondingly, no matter how small the step size of conventional predictor-corrector CPF is, the points are still sparse near the voltage bifurcation point, as shown in Figure 7.  When the reactive power of some PV buses and slack bus violate their limits, they are converted to PQ buses or Qθ (Qθ is used to describe bus types, which means that the reactive power and phase angle of the bus are known) bus respectively. In the nine-bus system, the conventional predictor-corrector method cannot smoothly realize the PV-PQ bus switching. The step size is s = 0.01. The PV curves of bus 9 for both methods are as follows. The points of PV curve are sparse near the voltage bifurcation point by using the conventional predictor-corrector method, which is not suitable for obtaining accurate maximum load point. But the ANM can smoothly realize the PV-PQ bus switching and the step size is automatically small near the voltage bifurcation point because of the self-adaptive adjustment. The PV curve can be smoothly drawn, as shown in Figure 6. Correspondingly, no matter how small the step size of conventional predictor-corrector CPF is, the points are still sparse near the voltage bifurcation point, as shown in Figure 7.    Step size 0.001 s = for conventional predictor-corrector method.

Impact of Slack Variable μ on CPF
The above calculation results are calculated when the slack variable is 0 μ = in Equation (5). In order to further analyze the influence of slack variable value on the calculation result, Step size s = 0.001 for conventional predictor-corrector method.

Impact of Slack Variable µ on CPF
The above calculation results are calculated when the slack variable is µ = 0 in Equation (5). In order to further analyze the influence of slack variable value on the calculation result, µ = 10 −7 and µ = 10 −20 are taken successively. We selected the voltage magnitude of bus 1 for both methods to compare. The voltage bifurcation point calculated by ANM is marked in Figures 8 and 9. The comparison results are as follows:  In view of Figures 8 and 9, we can conclude that when the slack variable μ takes different values, it has different effects on results. The PV curve of conventional CPF is the real solution curve.  In view of Figures 8 and 9, we can conclude that when the slack variable µ takes different values, it has different effects on results. The PV curve of conventional CPF is the real solution curve. When µ = 10 −7 , ANM will calculate the wrong voltage bifurcation point; when µ = 10 −20 , calculated solutions are almost the same as the real solutions. The higher the accuracy of variables for ANM, the less influences on the results. At meanwhile, the calculated solutions are closer to the real solutions.
The comparison on the CPU time for different µ is listed in Table 2. We can see that when µ takes different values, there is little difference in CPU time consumed. But ANM still consumes less time and computes more efficiently than conventional predictor-corrector method.

Computation of ANM-CPF
For a 182-bus network, the results of conventional predictor-corrector method are not convergent with the consideration of reactive power limit when the scalar embedding parameter λ increases to λ = 0.01. The PV buses or slack bus can be converted to PQ buses or a Qθ () bus when the reactive power violates their limits by using conventional predictor-corrector method. However, the switch cannot be recovered if the violation temporarily disappears in the sequential iterations. Once the PV buses or slack bus are converted to PQ buses or a Qθ bus, they cannot be converted back to PV buses or slack bus when the reactive power do not violate the limit, which is not consistent with the actual situation. But ANM can calculate the correct solutions when the scalar embedding parameter λ increases to λ = 1.340, the reactive power of bus 70 will violate its reactive power limit causing the system voltage to collapse. We calculated the PV curve of PQ buses, and selected two buses to display, as shown in Figure 10. It further proved that ANM is more robust than conventional predictor-corrector method from the above discussion. , the reactive power of bus 70 will violate its reactive power limit causing the system voltage to collapse. We calculated the PV curve of PQ buses, and selected two buses to display, as shown in Figure 10. It further proved that ANM is more robust than conventional predictor-corrector method from the above discussion.

Computation Failures of Conventional CPF
Conventional buses type switch logic methods in CPF, as the number of PV buses changing to PQ increases, the use of logical judgment inevitably increases the burden of calculation, even causing the power flow calculation failing or converging to the wrong solution. When 0.01 λ = , from Figure  11, we can see the PV-PQ bus switching strategy leads the conventional predictor-corrector power flow diverging in the large system, and the conventional predictor-corrector method first changed PV buses 59, 65, 103, 19, 34, 92, and 105 to PQ buses at position ①, then converged in five iterations.

Computation Failures of Conventional CPF
Conventional buses type switch logic methods in CPF, as the number of PV buses changing to PQ increases, the use of logical judgment inevitably increases the burden of calculation, even causing the power flow calculation failing or converging to the wrong solution. When λ = 0.01, from Figure 11, we can see the PV-PQ bus switching strategy leads the conventional predictor-corrector power flow diverging in the large system, and the conventional predictor-corrector method first changed PV buses 59, 65, 103, 19, 34, 92, and 105 to PQ buses at position 1 , then converged in five iterations. On sequential iteration, PV buses 49, 55, 56, 61, and 62 changed to PQ buses at position 2 and converged in four iterations; PV bus 54 became a PQ bus at position 3 ; and, finally, PV bus 66 became a PQ bus at position 4 . However, due to the number constantly oscillated during the iterative process, the conventional predictor-corrector power flow solution did not converge after 50 iterations from 4 . When λ = 0.01, the correct solutions are that PV buses 59, 65, 56, 55, 61, 62, 54, 49, 103, and 66 are changed to PQ buses. Compared to the real solutions, we can see that the conventional predictor-corrector method mistakenly changes the PV buses 19, 34, 92, and 105 to PQ buses because of false logic judgment, which leads to the wrong results during the iteration. Now let us explain why our proposed method is feasible and effective in the CPF calculation process considering non-smooth constraints exchange issue.
There is a sharp corner at the constraints exchange stage from

The Reason for the Proposed Method Working
Now let us explain why our proposed method is feasible and effective in the CPF calculation process considering non-smooth constraints exchange issue.
There is a sharp corner at the constraints exchange stage from Q i − Q min i to U + i = 0, where the conventional CPF method cannot obtain an accurate direction, which can be mistakenly fixed to PQ buses and cannot go back to PV buses. However, the proposed ANM-CPF based on FB function can solve the problem of the constraints exchange positions at the sharp corner better. It can obtain its local maximum successfully near the sharp point since the direction becomes more smoothly.

Conclusions
In this paper, we proposed a CPF calculation method based on the asymptotic numerical method. The algorithm step size can be adaptively adjusted, in addition we introduced a new variable for substitution to meet the calculation requirements of ANM, thus effectively dealing with the problem of reactive power limit. It has certain practicability in solving the problem of CPF in power systems. The small and large power system verify the robustness and efficiency of the proposed method.

Abbreviations
The following abbreviations are used in this manuscript: CPF continuation power flow ANM asymptotic numerical method PF power flow FB Fischer-Burmeister CPU Central processing unit PV PV is used to describe bus types, which means that the active power and voltage of the bus are known. PQ PQ is used to describe bus types, which means that the active and reactive power of bus are known. Qθ Qθ is used to describe bus types, which means that the reactive power and phase angle of the bus are known

Appendix A
Here we use a two-bus example to add explanation of ANM calculation principles, as shown in Figure A1.

Abbreviations
The following abbreviations are used in this manuscript: CPF continuation power flow ANM asymptotic numerical method PF power flow FB Fischer-Burmeister CPU Central processing unit PV PV is used to describe bus types, which means that the active power and voltage of the bus are known. PQ PQ is used to describe bus types, which means that the active and reactive power of bus are known. Qθ Qθ is used to describe bus types, which means that the reactive power and phase angle of the bus are known

Appendix A
Here we use a two-bus example to add explanation of ANM calculation principles, as shown in Figure A1.  Figure A1. Two-bus power system. In the two-bus power system, bus 1 is a slack bus and bus 2 is a PV bus. Set the parameter values: Node Admittance Matrix: Define variable x: where x 1∼4 represent the real and imaginary parts of bus 1,2, x 5∼6 represent the reactive power of bus 1,2, x 7∼10 are slack variables for voltage regulation, x 11∼12 are the voltage amplitude of bus 1,2, and x 13∼16 are slack variables for variable substitution.
The power flow equations are as follows: Further decompose the above Equation (A2) by a homotopy transformation: In the two-bus example, the matrix L is composed of one term in the PF constraint equation. The matrix Q is composed of quadratic terms in the PF equation. The matrix F represents the initial values of bus active and reactive power, and the matrix H represents constants, such as reactive power limit, slack variable µ, initial voltage value, etc. The detailed values are shown below.
where A L is a sparse coefficient matrix: The storage format of the coefficient matrix is defined: A[a, b, y] means that the a th row and b th column of the matrix is y. The values of the rest of elements in the matrix are zero.
The values of the sparse coefficient matrix A L are as follows: Taylor coefficients are calculated by using Equations (11) and (12). Then we calculate the maximum step size ∆s max . We can get the value of Q(x (r) , x (p−r) ) after the calculation of Taylor coefficients. So we can get the value of In view of (14), we can get the maximum step size ∆s max : The detailed step size for each point between 1th and 2th step is listed in Table A1. Figure out every point between 1 th and 2 th step by using (9).

Table A1.
Step size for each point between 1th and 2th.

Iteration Number q
Step Size (∆s)