Next Article in Journal
Hour-Ahead Energy Trading Management with Demand Forecasting in Microgrid Considering Power Flow Constraints
Previous Article in Journal
Flexibility Services to Minimize the Electricity Production from Fossil Fuels. A Case Study in a Mediterranean Small Island
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

College of Information and Electrical Engineering, China Agricultural University, Haidian District, Beijing 100083, China
*
Author to whom correspondence should be addressed.
Energies 2019, 12(18), 3493; https://doi.org/10.3390/en12183493
Submission received: 19 August 2019 / Revised: 4 September 2019 / Accepted: 5 September 2019 / Published: 10 September 2019

Abstract

:
Continuation power flow (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 first 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 effectively solve the PV-PQ bus type switching problem in CPF calculation. Furthermore, to assure high efficiency, 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.

1. 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 find the Q-limit breaking point indices. Then the algorithmic continuation steps would be guided by skipping to the closest Q-limit breaking point, consequently reducing the number of continuation steps and saving computational time. With the development of power system research, complementary constraints have been introduced into CPF calculation to deal with the problem of reactive power limit, which greatly improves the solving efficiency.
Reactive power limits can be formulated as complementarity constraints. By using complementary constraints to describe the switching relationship of PV-PQ bus type, the relationship between reactive power and voltage of the bus can be effectively considered without making logical judgment in the process of power flow (PF) iteration. References [12,13] stated that limit-induced constraints on exchange may induce PF divergence. An optimization based model was proposed in [14] to cope with complementarity constraints and was solved using interior point method (IPM). This method is relatively time-consuming since the PF equations are replaced by a constrained optimization problem. To guarantee PF global convergence, in [15,16], trust region based method was proposed to cope with complementarity constraints, whose computation burden even more heavy.
Predictor-corrector method is a conventional CPF calculation method. It involves the prediction of next solution point and correcting the prediction to get the next point on the curve. It includes four parts: prediction, correction, step size control, and parameterization. The values of variables and parameters along the solution curve can be parameterized in a number of ways [17,18].
The contributions of this paper are as follows:
(1)
In this paper, a new CPF calculation method based on asymptotic numerical method (ANM) is proposed, which can only cope with quadratic functions. PV-PQ bus type switch is reformulated by Fischer-Burmeister (FB) functions. FB function is also transformed into quadratic equation to meet the requirements of ANM. Then, ANM is used to calculate the continuous power flow and more robust than predictor-corrector based method to cope with constraints exchange issue. 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 failed.
(2)
Compared with the conventional predictor-corrector method, ANM can quickly approach the voltage collapse point by self-adaptive step size adjustment and constant Jacobian matrix used for power series expansion, thus greatly improving the CPF calculation efficiency.

2. Establishment and Solution of ANM-CPF

2.1. 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:
P i e i j = 1 N ( G i j e j B i j f j ) f i j = 1 N ( G i j f j + B i j e j ) = 0 ,
Q i f i j = 1 N ( G i j e j B i j f j ) + e i j = 1 N ( G i j f j + B i j 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 i j , B i j are conductance and susceptance matrix elements respectively.
For PV buses, Equation (2) can be replaced by the equation of voltage:
( U i s e t ) 2 e i 2 f i 2 = 0 ,
where U i s e t 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:
{ e i 2 + f i 2 = U i s e t + U i + U i , U i + 0 , U i 0 ( Q i Q i min ) U i + = 0 , i { P V , S l a c k } ( Q i max Q i ) U i = 0 , i { P V , S l a c k } Q i min Q i Q i max ,
where Q i max and Q i min 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:
ϕ ( u , v ) = u + v u 2 + v 2 + μ , u 0 , v 0 ,
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:
ϕ ( u , v ) = u + v w , u 0 , v 0 , w 0 ,
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:
{ U i + + ( Q i Q i min ) w + = 0 ( w + ) 2 ( ( U i + ) 2 + ( Q i Q i min ) 2 + μ ) = 0 U i + ( Q i max Q i ) w = 0 ( w ) 2 ( ( U i ) 2 + ( Q i max Q i ) 2 + μ ) = 0 ,
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.

2.2. 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:
f ( x , λ ) = L ( x ) + Q ( x , x ) + λ F + H = 0 ,
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 i min , Q i max , 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 t h   ( t = 0 , 1 , 2 ) computed point, the q t h point between t t h and ( t + 1 ) t h step can be expressed as the series like:
x t q = x t + p = 1 K x t ( p ) ( Δ s ) q λ t q = λ t + p = 1 K λ t ( p ) ( Δ s ) q ,
where K is the truncation order of the series, ( Δ s ) q is the step size of q t h point between t t h and ( t + 1 ) t h , ( x t ( p ) , λ t ( p ) ) is the Taylor coefficients.
The maximum step size is divided into M equal parts. And the step size of q t h point between t t h and ( t + 1 ) t h satisfies:
( Δ s ) q = q Δ s max M , 1 q M ,
where Δ s max is the maximum step size of the t t h computed point and when q is equal to M , x t q = 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 t ( p ) , λ t ( p ) ) at different orders are determined by a series of equations:
Order p = 1 :
f x t v 1 = F ( x t ( 1 ) , λ t ( 1 ) ) = ± ( v 1 , 1 ) / 1 + v 1 2 2 ,
Order p 2 :
f x t v p = r = 1 p 1 Q ( x t ( r ) , x t ( p r ) ) λ t ( p ) = λ t ( 1 ) v p , x t ( 1 ) x t ( p ) = v p + λ t ( p ) v 1 ,
where f x t 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.
Notice that Q ( x t ( r ) , x t ( p r ) ) is a matrix composed of functions that construct quadratic terms of the equations:
Q ( x t ( r ) , x t ( p r ) ) = [ ( x t ( r ) ) T * A Q 1 * ( x t ( p r ) ) ( x t ( r ) ) T * A Q n * ( x t ( p r ) ) ] ,
where A Q h   ( 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:
Δ s max = ( ε / ( r = 1 K 1 Q ( x t ( r ) , x t ( K r ) ) ) 1 K ,
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:
Δ s max = ( ε / ( f x t v K ) ) 1 K .

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

3.1. Example 1: IEEE Nine-Bus Network

3.1.1. 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 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.
In view of (15), ( f x t 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.

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

3.1.3. 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 Figure 8 and Figure 9. The comparison results are as follows:
In view of Figure 8 and Figure 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.

3.2. Example 2: A Large 182-Bus Network

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

3.2.2. 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. On sequential iteration, PV buses 49, 55, 56, 61, and 62 changed to PQ buses at position ② and converged in four iterations; PV bus 54 became a PQ bus at position ③; and, finally, PV bus 66 became a PQ bus at position ④. 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 ④. 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.

3.2.3. 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 i min 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.

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

Author Contributions

Y.J. and Y.H.; methodology, Y.J.; software, Z.Z.; validation, Y.H., Y.J., Z.Z.; formal analysis, Y.J.; investigation, Y.H.; resources, Y.J.; data curation, Z.Z.; writing—original draft preparation, Y.H.; writing—review and editing, Y.H.; supervision, Y.J.; project administration, Y.J.; funding acquisition, Y.J.

Funding

This paper is supported by the National Natural Science Foundation of China (Grant No. 51707196) and the State Grid Science and Technology Program of China “On-line Risk Analysis and Risk Pre-control Key Technologies of Transport and Distribution Collaboration for Large-scale New Energy Security Disposal”.

Acknowledgments

The calculation resources for this research were provided by CAU Smart Energy System. Yuntao Ju also belongs to Beijing Key Laboratory of Demand Side Multi-Energy Carriers Optimization and Interaction Technique, Beijing 100192, China. The authors would like to thank the unknown reviewers for their work and suggestion regarding the paper.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
CPFcontinuation power flow
ANMasymptotic numerical method
PFpower flow
FBFischer-Burmeister
CPUCentral processing unit
PVPV is used to describe bus types, which means that the active power and voltage of the bus are known.
PQPQ 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.
Figure A1. Two-bus power system.
Energies 12 03493 g0a1
In the two-bus power system, bus 1 is a slack bus and bus 2 is a PV bus. Set the parameter values:
P = 1.5 , V 2 = 0.95 θ , Q 1 max = Q 2 max = 15 , Q 1 min = Q 2 min = 15
Node Admittance Matrix:
G + j B = [ j 10 j 10 j 10 j 10 ]
Define variable x :
x = [ e 1 , f 1 , e 2 , f 2 , Q 1 , Q 2 , V 1 + , V 1 , V 2 + , V 2 , | V 1 | , | V 2 | , w 1 + , w 1 , w 2 + , w 2 ] T
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:
{ f 1 = 0 B 21 f 2 + 5 ( λ + 1 ) P = 0 Q 1 ( B 11 e 1 2 B 11 f 1 2 B 12 e 1 e 2 B 12 f 1 f 2 ) = 0 Q 2 ( B 21 e 2 2 B 21 f 2 2 B 22 e 1 e 2 B 22 f 1 f 2 ) = 0 e 1 2 + f 1 2 | V 1 | 2 = 0 e 2 2 + f 2 2 | V 2 | 2 = 0 V 1 + V 1 + V 1 | V 1 | = 0 V 2 + V 2 + V 2 | V 2 | = 0 V 1 + + ( Q 1 Q 1 min ) w 1 + = 0 ( w 1 + ) 2 ( ( V 1 + ) 2 + ( Q 1 Q 1 min ) 2 + μ ) = 0 V 1 + ( Q 1 max Q 1 ) w 1 = 0 ( w 1 ) 2 ( ( V 1 ) 2 + ( Q 1 max Q 1 ) 2 + μ ) = 0 V 2 + + ( Q 2 Q 2 min ) w 2 + = 0 ( w 2 + ) 2 ( ( V 2 + ) 2 + ( Q 2 Q 2 min ) 2 + μ ) = 0 V 2 + ( Q 2 max Q 2 ) w 2 = 0 ( w 2 ) 2 ( ( V 2 ) 2 + ( Q 2 max Q 2 ) 2 + μ ) = 0 .
Further decompose the above Equation (A2) by a homotopy transformation:
f ( x , λ ) = L ( x ) + Q ( x , x ) + λ F + H = 0 .
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.
L ( x ) = [ f 1 B 21 f 2 0 0 0 0 V 1 + V 1 | V 1 | V 2 + V 2 | V 2 | V 1 + + Q 1 W 1 + 2 Q 1 min Q 1 V 1 Q 1 W 1 2 Q 1 max Q 1 V 2 + + Q 2 W 2 + 2 Q 2 min Q 2 V 2 Q 2 W 2 2 Q 2 max Q 2 ] = A L * x ,
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 t h row and b t h 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:
A L : [ 1 , 2 , 1 ] , [ 2 , 4 , B 21 ] , [ 7 , 7 , 1 ] , [ 7 , 8 , 1 ] , [ 7 , 11 , 1 ] , [ 8 , 8 , 1 ] , [ 8 , 9 , 1 ] , [ 9 , 5 , 1 ] , [ 9 , 7 , 1 ] , [ 9 , 13 , 1 ] , ( 10 , 5 , 2 Q 1 min ] , [ 11 , 5 , 1 ] , [ 11 , 8 , 1 ] , [ 11 , 14 , 1 ] , [ 12 , 5 , 2 Q 1 max ] , [ 13 , 6 , 1 ] , [ 13 , 9 , 1 ] , [ 13 , 15 , 1 ] , [ 14 , 6 , 2 Q 2 min ] , [ 15 , 6 , 1 ] , [ 15 , 10 , 1 ] , [ 15 , 16 , 1 ] , [ 16 , 6 , 2 Q 2 max ]
Q ( x , x ) = [ 0 0 B 11 e 1 2 + B 11 f 1 2 + B 12 e 1 e 2 + B 12 f 1 f 2 B 21 e 2 2 + B 21 f 2 2 + B 22 e 1 e 2 + B 22 f 1 f 2 e 1 2 + f 1 2 | V 1 | 2 e 2 2 + f 2 2 | V 2 | 2 0 0 0 ( W 1 + ) 2 ( V 1 + 2 + Q 1 2 ) 0 ( W 1 ) 2 ( V 1 2 + Q 1 2 ) 0 ( W 2 + ) 2 ( V 2 + 2 + Q 2 2 ) 0 ( W 2 ) 2 ( V 2 2 + Q 2 2 ) ] = [ x T * A Q 1 * x x T * A Q 16 * x ] ,
where A Q h   ( 1 h 16 ) are sparse coefficient matrices.
The values of the sparse coefficient matrix A Q h are as follows:
  • A Q 3 : [ 1 , 1 , B 11 ] , [ 2 , 1 , B 12 ] , [ 2 , 2 , B 11 ] , [ 4 , 2 , B 12 ] ;
  • A Q 4 : [ 1 , 3 , B 22 ] , [ 2 , 4 , B 22 ] , [ 3 , 2 , B 21 ] , [ 4 , 4 , B 21 ] ;
  • A Q 5 : [ 1 , 1 , 1 ] , [ 2 , 2 , 1 ] , [ 11 , 11 , 1 ] ;
  • A Q 6 : [ 3 , 3 , 1 ] , [ 4 , 4 , 1 ] , [ 12 , 12 , 1 ] ;
  • A Q 10 : [ 5 , 5 , 1 ] , [ 7 , 7 , 1 ] , [ 13 , 13 , 1 ] ;
  • A Q 12 : [ 5 , 5 , 1 ] , [ 8 , 8 , 1 ] , [ 14 , 14 , 1 ] ;
  • A Q 14 : [ 6 , 6 , 1 ] , [ 9 , 9 , 1 ] , [ 15 , 15 , 1 ] ;
  • A Q 16 : [ 6 , 6 , 1 ] , [ 10 , 10 , 1 ] , [ 16 , 16 , 1 ] ;
    F = [ 0 , 5 P , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ] T = [ 0 , 7.5 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ] T ,
    H = [ 0 , 0 , 0 , 0 , 0 , 0 , V 1 , V 2 , Q 1 min , ( Q 1 min ) 2 μ , Q 1 max , ( Q 1 max ) 2 μ , Q 2 min , ( Q 2 min ) 2 μ , Q 2 max , ( Q 2 max ) 2 μ ] .
Set the parameter values: ε = 1.0 × 10 6 , K = 20 , M = 15 . Take the calculation of 1 t h point as an example.
The value of Jacobian matrix f x 1 at 1 t h point is as follows:
f x 1 : [ 3 , 1 , 14.169048 ] , [ 4 , 1 , 5.8309519 ] , [ 5 , 1 , 2.0000000 ] , [ 1 , 2 , 1.0000000 ] , [ 3 , 2 , 7.5000000 ] , [ 4 , 2 , 7.5000000 ] , [ 3 , 3 , 10.000000 ] , [ 4 , 3 , 1.6619038 ] , [ 6 , 3 , 1.1661904 ] , [ 2 , 4 , 10.000000 ] , [ 4 , 4 , 15.000000 ] [ 6 , 4 , 1.5000000 ] , [ 3 , 5 , 1.0000000 ] , [ 9 , 5 , 1.0000000 ] , [ 10 , 5 , 38.338096 ] , [ 11 , 5 , 1.0000000 ] , [ 12 , 5 , 21.661904 ] , [ 4 , 6 , 1.000000 ] , [ 13 , 6 , 1.00000 ] , [ 14 , 6 , 36.388097 ] , [ 15 , 6 , 1.0000000 ] , [ 7 , 7 , 1.0000000 ] , [ 9 , 7 , 1.0000000 ] , [ 7 , 8 , 1.00000 ] , [ 11 , 8 , 1.0000000 ] , [ 8 , 9 , 1.0000000 ] , [ 13 , 9 , 1.0000000 ] , [ 8 , 10 , 1.0000000 ] , [ 15 , 10 , 1.00000 ] , [ 5 , 11 , 2.00000 ] , [ 7 , 11 , 1.00000 ] , [ 6 , 12 , 1.900000 ] , [ 9 , 13 , 1.000000 ] , [ 10 , 13 , 38.338096 ] , [ 11 , 14 , 1.000 ] , [ 12 , 14 , 21.661904 ] , [ 13 , 15 , 1.000 ] , [ 14 , 15 , 36.388096 ] , [ 15 , 16 , 1.000 ] , [ 16 , 16 , 23.611904 ] , .
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 r = 1 K 1 Q ( x j ( r ) , x j ( K r ) ) as follows:
r = 1 K 1 Q ( x j ( r ) , x j ( K r ) ) = 1.0 × 10 20 × [ 0 , 0 , 0 , 0.30175640 , 0 , 0.03017564 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ] T .
Then we can get r = 1 K 1 Q ( x j ( r ) , x j ( K r ) ) = 3.0326143 × 10 21 .
In view of (14), we can get the maximum step size Δ s max :
Δ s max = 5.3199695 .
The detailed step size for each point between 1th and 2th step is listed in Table A1. Figure out every point between 1 t h and 2 t h step by using (9).
Table A1. Step size for each point between 1th and 2th.
Table A1. Step size for each point between 1th and 2th.
Iteration Number q Step Size ( Δ s ) 1 q
10
20.3799978
30.7599956
41.1399935
51.5199913
61.8999891
72.2799869
82.6599847
93.0399826
103.4199804
113.7999782
124.1799760
134.5599739
144.9399717
155.3199695

References

  1. Wang, X.-F.; Song, Y.-H.; Irving, M. Modern Power Systems Analysis; Springer: New York, NY, USA, 2008; ISBN 978-0-387-72852-0. [Google Scholar]
  2. Feng, Z.; Ajjarapu, V.; Long, B. Identification of voltage collapse through direct equilibrium tracing. IEEE Trans. Power Syst. 2000, 15, 342–349. [Google Scholar] [CrossRef]
  3. Fnaiech, N.; Jendoubi, A.; Zoghlami, M.; Bacha, F. Continuation power flow of voltage stability limits and a three dimensional visualization approach. In Proceedings of the 2015 16th International Conference on Sciences and Techniques of Automatic Control and Computer Engineering (STA), Monastir, Tunisia, 21–23 December 2015; pp. 163–168. [Google Scholar]
  4. Kumar, S.; Kumar, A.; Sharma, N.K. Analysis of power flow, continuous power flow and transient stability of IEEE-14 bus integrated wind farm using PSAT. In Proceedings of the 2015 International Conference on Energy Economics and Environment (ICEEE), Greater Noida, India, 27–28 March 2015; pp. 1–6. [Google Scholar]
  5. Garcia, P.A.N.; Vinagre, M.P.; Oliveira, E.J. Fault analysis using continuation power flow and phase coordinates. In Proceedings of the IEEE Power Engineering Society General Meeting, Denver, CO, USA, 6–10 June 2004; Volume 2, pp. 872–874. [Google Scholar]
  6. Wang, M.; Xia, Y.; Chen, Y.; Huang, S. GPU-based power flow analysis with continuous Newton’s method. In Proceedings of the 2017 IEEE Conference on Energy Internet and Energy System Integration (EI2), Beijing, China, 26–28 November 2017; pp. 1–5. [Google Scholar]
  7. Cao, B.; Bo, M.; Pan, J.; Liu, Y. Static Voltage Stability Analysis Based on the Combination of Dynamic Continuous Power Flow and Adaptive Chaotic Particle Swarm Optimization. In Proceedings of the 2018 IEEE 3rd Advanced Information Technology, Electronic and Automation Control Conference (IAEAC), Chongqing, China, 12–14 October 2018; pp. 2217–2221. [Google Scholar]
  8. Canizares, C.A.; Alvarado, F.L. Point of collapse and continuation methods for large AC/DC systems. IEEE Trans. Power Syst. 1993, 8, 1–8. [Google Scholar] [CrossRef]
  9. Zhao, J.; Zhang, B. Reasons and countermeasures for computation failures of continuation power flow. In Proceedings of the 2006 IEEE Power Engineering Society General Meeting, Montreal, QC, Canada, 18–22 June 2006; p. 6. [Google Scholar]
  10. Zhao, J.; Zhou, C.; Chen, G. A novel bus-type extended Continuation Power Flow considering remote voltage control. In Proceedings of the 2013 IEEE Power & Energy Society General Meeting, Vancouver, BC, Canada, 21–25 July 2013; pp. 1–5. [Google Scholar]
  11. Zhu, P.; Taylor, G.; Irving, M. A novel Q-Limit guided Continuation Power Flow method. In Proceedings of the 2008 IEEE Power and Energy Society General Meeting—Conversion and Delivery of Electrical Energy in the 21st Century, Pittsburgh, PA, USA, 20–24 July 2008; pp. 1–7. [Google Scholar]
  12. Yue, X.; Venkatasubramanian, V.M. Complementary limit induced bifurcation theorem and analysis of Q limits in power-flow studies. In Proceedings of the 2007 iREP Symposium—Bulk Power System Dynamics and Control—VII. Revitalizing Operational Reliability, Charleston, SC, USA, 19–24 August 2007; pp. 1–8. [Google Scholar]
  13. Avalos, R.J.; Canizares, C.A.; Milano, F.; Conejo, A.J. Equivalency of Continuation and Optimization Methods to Determine Saddle-Node and Limit-Induced Bifurcations in Power Systems. IEEE Trans. Circuits Syst. I Regul. Pap. 2009, 56, 210–223. [Google Scholar] [CrossRef]
  14. Pirnia, M.; Cañizares, C.A.; Bhattacharya, K. Revisiting the power flow problem based on a mixed complementarity formulation approach. IET Gener. Transm. Distrib. 2013, 7, 1194–1201. [Google Scholar] [CrossRef]
  15. De Rubira, T.T.; Wigington, A. Extending complementarity-based approach for handling voltage band regulation in power flow. In Proceedings of the 2016 Power Systems Computation Conference (PSCC), Genoa, Italy, 20–24 June 2016; pp. 1–6. [Google Scholar]
  16. Lin, Y.; Ju, Y. A Robust Three-Phase Power Flow for Active Distribution Network Embedded Autonomous Voltage Regulating Strategies. In Proceedings of the 2018 2nd IEEE Conference on Energy Internet and Energy System Integration (EI2), Beijing, China, 20–22 October 2018; pp. 1–8. [Google Scholar]
  17. Chiang, H.-D.; Flueck, A.J.; Shah, K.S.; Balu, N. CPFLOW: A practical tool for tracing power system steady-state stationary behavior due to load and generation variations. IEEE Trans. Power Syst. 1995, 10, 623–634. [Google Scholar] [CrossRef]
  18. Li, S.-H.; Chiang, H.-D. Nonlinear predictors and hybrid corrector for fast continuation power flow. IET Gener. Transm. Distrib. 2008, 2, 341. [Google Scholar] [CrossRef]
  19. Fischer, A. A special newton-type optimization method. Optimization 1992, 24, 269–284. [Google Scholar] [CrossRef]
  20. Yang, X.; Zhou, X. Application of asymptotic numerical method with homotopy techniques to power flow problems. Int. J. Electr. Power Energy Syst. 2014, 57, 375–383. [Google Scholar] [CrossRef]
  21. Zimmerman, R.D.; Murillo-Sanchez, C.E.; Thomas, R.J. MATPOWER: Steady-State Operations, Planning, and Analysis Tools for Power Systems Research and Education. IEEE Trans. Power Syst. 2011, 26, 12–19. [Google Scholar] [CrossRef]
  22. Ju, Y. Constraints Exchange Data. Available online: https://gitee.com/cieeyodaju/Constraints-Exchange-Data (accessed on 6 August 2019).
  23. Ruan, C.; Wang, X.; Wang, X.; Gao, F.; Li, Y. Improved Continuation Power Flow Calculation Method Based on Coordinated Combination of Parameterization. In Proceedings of the 2018 IEEE 2nd International Electrical and Energy Conference (CIEEC), Beijing, China, 4–6 November 2018; pp. 207–211. [Google Scholar]
  24. Alam, M.M.; Moreira, C.; Islam, M.R.; Mehedi, I.M. Continuous Power Flow Analysis for Micro-Generation Integration at Low Voltage Grid. In Proceedings of the 2019 International Conference on Electrical, Computer and Communication Engineering (ECCE), Cox’sBazar, Bangladesh, 7–9 February 2019; pp. 1–5. [Google Scholar]
Figure 1. V-lambda curve of bus 9 for both methods.
Figure 1. V-lambda curve of bus 9 for both methods.
Energies 12 03493 g001
Figure 2. Voltage transient stability analysis results of bus 1, 2, and 3.
Figure 2. Voltage transient stability analysis results of bus 1, 2, and 3.
Energies 12 03493 g002
Figure 3. Vertical increment at different points.
Figure 3. Vertical increment at different points.
Energies 12 03493 g003
Figure 4. Maximum step size at different points.
Figure 4. Maximum step size at different points.
Energies 12 03493 g004
Figure 5. V-lambda curve of bus 1.
Figure 5. V-lambda curve of bus 1.
Energies 12 03493 g005
Figure 6. V-lambda curve of bus 9 for both methods.
Figure 6. V-lambda curve of bus 9 for both methods.
Energies 12 03493 g006
Figure 7. Step size s = 0.001 for conventional predictor-corrector method.
Figure 7. Step size s = 0.001 for conventional predictor-corrector method.
Energies 12 03493 g007
Figure 8. V-lambda curve of bus 1 when μ = 10 7 .
Figure 8. V-lambda curve of bus 1 when μ = 10 7 .
Energies 12 03493 g008
Figure 9. V-lambda curve of bus 1 when μ = 10 20 .
Figure 9. V-lambda curve of bus 1 when μ = 10 20 .
Energies 12 03493 g009
Figure 10. V-lambda curve of bus 30 and 58.
Figure 10. V-lambda curve of bus 30 and 58.
Energies 12 03493 g010
Figure 11. Conventional CPF corrector iteration at λ = 0.01 .
Figure 11. Conventional CPF corrector iteration at λ = 0.01 .
Energies 12 03493 g011
Table 1. Calculation time of one point in PV curve for both methods.
Table 1. Calculation time of one point in PV curve for both methods.
MethodANM-CPFConventional CPF
Average CPU time for each point/s 4.01 × 10 4 2.34 × 10 2
Table 2. Comparison on the central processing unit (CPU) time for different μ .
Table 2. Comparison on the central processing unit (CPU) time for different μ .
Slack Variable μ = 10 7 μ = 10 20
Average CPU time for each point/s 3.60 × 10 4 3.52 × 10 4

Share and Cite

MDPI and ACS Style

Huang, Y.; Ju, Y.; Zhu, Z. An Asymptotic Numerical Continuation Power Flow to Cope with Non-Smooth Issue. Energies 2019, 12, 3493. https://doi.org/10.3390/en12183493

AMA Style

Huang Y, Ju Y, Zhu Z. An Asymptotic Numerical Continuation Power Flow to Cope with Non-Smooth Issue. Energies. 2019; 12(18):3493. https://doi.org/10.3390/en12183493

Chicago/Turabian Style

Huang, Yan, Yuntao Ju, and Zeping Zhu. 2019. "An Asymptotic Numerical Continuation Power Flow to Cope with Non-Smooth Issue" Energies 12, no. 18: 3493. https://doi.org/10.3390/en12183493

APA Style

Huang, Y., Ju, Y., & Zhu, Z. (2019). An Asymptotic Numerical Continuation Power Flow to Cope with Non-Smooth Issue. Energies, 12(18), 3493. https://doi.org/10.3390/en12183493

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop