- Email: [email protected]
A direct method based on tensor calculation to determine maximum loadability power flow solutions
Electric Power Systems Research 103 (2013) 114–121
Contents lists available at SciVerse ScienceDirect
Electric Power Systems Research journal homepage: www.elsevier.com/locate/epsr
A direct method based on tensor calculation to determine maximum loadability power flow solutions Roberto S. Salgado ∗ , Anesio F. Zeitune Department of Electrical Engineering, Power System Laboratory (Labspot), Federal University of Santa Catarina, Florianópolis, SC, Brazil
a r t i c l e
i n f o
Article history: Received 9 November 2012 Received in revised form 15 May 2013 Accepted 17 May 2013 Available online 15 June 2013 Keywords: Maximum loadability Direct methods Tensor term Saddle-node bifurcation Sensitivity analysis
a b s t r a c t In this paper, the determination of critical loadability points of the power flow equations is formulated as an optimisation problem. A quadratic parameterisation scheme is used to model the load change, which is equivalent to implicitly including a non-negativity constraint in the load variation. The power flow equations are expressed in rectangular coordinates. This is suitable to exploit the second order information (tensor term) of the equations representing the optimality conditions. The use of this term improves the convergence of the iterative process and facilitates the manipulation of the reactive power generation limits. Simulation results obtained with a number of power systems, including real networks, are used to illustrate the main features of the proposed methodology. © 2013 Elsevier B.V. All rights reserved.
1. Introduction Several aspects of voltage instability have been explained by the Bifurcation Theory of nonlinear systems. In particular, two types of local bifurcations have been frequently associated with the steady state aspects of the voltage collapse [1,2]. The first is the saddle-node bifurcation, in which the Jacobian matrix of the differential-algebraic equations representing the power system is singular [3,4]. The second, named limit-induced bifurcation, is related to changes in the stability condition, as a consequence of a particular variable reaching its limit. In the present work, the power flow solution corresponding to the maximum (here also referred to as critical) loadability is characterised as a saddle-node bifurcation point (also named limit point, fold point and turning point). Although very interesting from the practical point of view of power system stability, other types of bifurcation are not focused on here. A number of methods for determining the critical loadability point can be found in the literature. Continuation methods have been widely used to determine a sequence of solutions for the parameterised power flow equations from a base load to the point of maximum loadability [5–7]. A predictor-corrector strategy is applied to calculate the set of equilibrium points that compose the PV-curve of each bus. Additionally, the predictor tangent vector is obtained as a by-product of the computational process, with
∗ Corresponding author. E-mail address: [email protected] (R.S. Salgado). 0378-7796/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.epsr.2013.05.011
reduced computational effort through the use of the techniques presented in [8]. Direct methods determine the saddle-node bifurcation point in one step, providing simultaneously the right or left eigenvectors of the singular power flow Jacobian matrix. Some of the relevant approaches found in the literature are mentioned as follows. In [9], an optimisation problem is solved by Newton’s method, in which (without loss of generality) the active power level is constant, such that only the QV subproblem is solved. In [7,10], the set of equations representing the Transversality Conditions are extended to integrate HVDC links, area interchange power control and other power systems particularities. The resulting set of nonlinear equations is solved to provide the so-called Point of Collapse. Refs. [11,12] propose direct methods to compute the saddle-node bifurcation point closest to the base case in the load parameter space. These methods combine the computation of both the load power margin and the vector normal to the limit surface of the feasible region of the power flow equations. They can provide the saddle-node bifurcation point in a fixed direction of load increase as well as the bifurcation point locally closest to the current operating load level. Ref. [13] proposes the determination of the loadability margin by solving the Moore–Spence augmented system. Its major contribution is the application of a matrix decomposition strategy to speed up the computational process. Ref. [14] combines the Continuation method with optimisation techniques to detect the limit of solvability of the power flow equations. A least squares technique is used to calculate the Lagrange multipliers, which are used as a figure of merit to find the power flow solvability limit. More recently, Yang et al. [15] proposes the use of block elimination techniques to solve the Moore–Spence augmented system. Similarly
R.S. Salgado, A.F. Zeitune / Electric Power Systems Research 103 (2013) 114–121
to [13], the contribution of this approach is related to the implementation aspects. The similarities and differences between these methodologies concern five main points: (1) the analytical formulation, (2) the type of load parameterisation, (3) the solution method, (4) the initial estimates, and (5) the by-products of the computation. First, the sets of equations solved to calculate fold points are basically those which represent alternatively the Transversality Conditions [7,12] the Moore–Spence determining system [13,15], or the optimality conditions of the maximum loadability optimisation problem [9]. In case of the optimisation model, adequate assumptions must be taken into account in order to obtain fold points similar to those produced by Continuation methods. Except for this aspect there are no considerable differences between these models. Second, with respect to the parameterisation, all of these methodologies have used the linear power load parameterisation to model the load increase. Unless a non-negativity constraint on the load parameter is added to the formulation, there is a risk of determining undesirable critical loadability solutions with load decrease (negative load parameter), if the initial estimate is not well chosen. Referring to the solution method, due to the particular structure and features of the maximum loadability problem, all of these approaches apply Newton’s method to solve the set of nonlinear equations that provide the fold point. This method is well defined and quadratically (although only locally) convergent. It is efficient to deal with equality and inequality constraints [9], but requires a good initial estimate. For this purpose, strategies like iterations of the Inverse Power Method have been used in [7,10], information about the largest eigenvalue of the Jacobian matrix available in the current operating point is employed in [11,12] and the Continuation method is applied in [13–15] to provide the initial guesses for power flow variables and eigenvectors. From this point of view, the exception is the approach presented in [9], which claims no need for the solution of any separate system to obtain improved initial estimates. The numerical results obtained through the application of these direct methods include the fold point and the sensitivity relationships represented by the eigenvectors and Lagrange multipliers, which are used to define critical buses and areas. The maximum loadability problem can be also formulated as a general Optimal Power Flow problem, as in [16,17], for example. The analytical models used in these approaches include other types of operational constraints (voltage magnitude, power flows, etc.) and optimisation variables (generated active and reactive powers, transformer taps, etc.). This brings additional difficulties (in terms of size and complexity of the problem) to obtain critical solutions and requires different numerical methods to solve the optimisation problem (as Interior Point Methods and/or Trust Region based-methods, for example). The maximum loadability solutions produced through the application of these methodologies is completely different from those focused on the present work. For this reason this type of approach is out of the scope of this paper. The present paper proposes a methodology for the direct determination of the turning points of the power flow equations. This is formulated as a constrained static optimisation problem, which is solved by an extended version of Newton’s method. For this purpose, the power flow equations are expressed in rectangular coordinates. The main contributions reported here are: (a) the use of an alternative load parameterisation scheme to model the load variation, so as to ensure power flow solutions with load level greater than that of the base case; (b) the inclusion of the secondorder information (tensor component) in the search direction, to improve the robustness of the iterative process and facilitate the treatment of the reactive power generation limits. The proposed approach differs from those previously mentioned with respect to: (1) the type of load parameterisation, (2) the strategy to solve the optimisation problem, which is an extension of Newton’s method
115
to exploit the second-order information. Similarly to [9], there is no need of any special strategy to choose initial guesses for the Lagrange multipliers, although we also show two alternative procedures to compute the initial estimates of the Lagrange multipliers. Numerical results obtained for power systems with sizes ranging from 24 to 1916 buses are used to illustrate the proposed methodology. This paper is organised as follows. Section 2 describes the theoretical basis of the maximum loadability problem, including the solution through static optimisation algorithms. Section 3 presents the analytical model proposed herein, with emphasis on the use of the information conveyed by the tensor term. Section 4 presents the numerical results obtained from the application of the methodology proposed here and Section 5 summarises the main conclusions. 2. Preliminary concepts The determination of the power flow solution at the maximum loadability level can be formulated as an optimisation problem expressed in terms of power system variables as: Maximize
subject to
(Pg0 + 2 Pgi ) − (Pd0 + 2 Pdi ) − Pi (e, f ) = i
i
Qg0 − (Qd0 + 2 Qdi ) − Qi (e, f ) = i
ref
Vi
i
2
0 (PV, PQ buses)
0 (PQ buses)
(1)
= ei2 + fi2 (PV buses)
Qgi = (Qd0 + 2 Qdi ) + Qi (e, f ) ≤ Qglim (PV buses) i
Pg0i
Qg0i
i
i
i
where and are the active and reactive power generation of bus i in the base case and Pgi represents the active power generation change rate, which indicates how the active power generation follows the increase of the loadability level (distribution factors based on economic and/or security criteria can be used for this purpose); Pd0 and Qd0 refer to the active and reactive power demand of the ith bus in the base case, Pdi and Qdi represent the prespecified direction of increase of the active and reactive power load ref of bus i, Vi is the reference voltage magnitude of the ith PV bus, and Pi (e, f) and Qi (e, f) are the nodal active and reactive power injections expressed as functions of the real and imaginary components of complex voltages. The inequality constraints refer to limits of generated reactive power of the PV buses. The optimisation variables are the components of the complex voltage (ei , fi ) and the load parameter . If the conventional linear parameterisation is used to model the load variation and a non-negativity constraint is not imposed on , the optimisation problem of Eq. (1) can provide undesirable critical solutions (with negative loading factor, for example), which are useless from the practical point of view. This is usually due to a poor initial estimate of the optimisation variables, in particular the Lagrange multipliers. In order to overcome this difficulty, the load change is quadratically parameterised, as shown in Eq. (1). This type of parameterisation is equivalent to imposing a non-negativity constraint in the load variation. The constraints of the optimisation problem of Eq. (1) are expressed in compact form as: g(x, )
=
y0 + 2 r − g0 (x) = 0
hqg (x, )
=
yq0 + 2 rqg − h0 (x) ≤
hlim qg
(2)
where the components of vectors (y0 + 2 r) and (yq0 + 2 rqg ) are related to the bus power injections, the squared voltage magnitude of the PV buses in the base case (y0 and yq0 ) and the load increase direction (r and rqg ). Since the PV bus voltage magnitude does not depend on the load variation (see the third equality constraint in Eq. (1)), the components of the vector r corresponding to the squared voltage magnitudes are zero. Vectors g0 (x) and h0 (x)
116
R.S. Salgado, A.F. Zeitune / Electric Power Systems Research 103 (2013) 114–121
represent the power injections and the PV buses squared voltage magnitude expressed in terms of vector x, which stands for the real and imaginary components of the complex voltages, and hlim qg represents the reactive power generation limits. If only the equality constraints (the treatment of inequality constraints will be described later) are taken into account, the optimality conditions for the problem stated by Eq. (1) are:
∇ x £(x, , ) = 0 =
J(x)t
∇ £(x, , ) = 0 =
− 1 + 2rt
∇ £(x, , ) = 0 =
− (y0 + 2 r − g0 (x))
(3)
The most important aspect of Eq. (6) is that there are no terms with order higher than the second one, so that given a direction d the function F(ze + d) can be accurately computed.
2.2. Newton’s method Newton’s method uses a linear approximation of Eq. (6), such that at each iteration the solution of the following linear system is required: W(ze )dn = −F(ze )
where J(x) is the Jacobian matrix of equations g0 (x) with respect to the variables x. Since the rate of load increase (r) is non-zero, from the second optimality condition of Eq. (3), = / 0. Thus, equation J(x)t = 0 indicates that the Jacobian matrix is singular at the point of critical loadability, and that vector is co-linear with the left eigenvector corresponding to the zero eigenvalue of matrix J(x).
(9)
where the subscript n denotes the component of the search direction obtained by Newton’s method. The non-singularity of the coefficient matrix W(ze ) of Eq. (9) (formally shown in [3]) allows the point of critical loadability to be accurately determined. To speed up the iterative process the decomposition strategy of this matrix, shown in [13,20], can be applied.
2.1. Use of the quadratic form If the power flow equations are expressed in rectangular coordinates, the functions g0 (x) of Eq. (2) can be represented by a quadratic function expressed as [18] 1 g0 (x) = xt T0 x 2 where T0 is a three-dimensional array, called the tensor, which depends only on the transmission line parameters; such that its first and second derivatives are given by:
∂g0 (x) = J(x) = T0 x and ∂x
3. Proposed methodology 3.1. Initial condition In general, the initial estimates of the Lagrange multipliers are considerably far from those of the critical loadability solution, which can be a drawback to the application of Newton’s method. To specify initial values for these dual variables, the literature [21] suggests the least squares solution given by:
2
∂ g0 (x) ∂J(x) = = T0 ∂x2 ∂x
(4)
The second order Taylor series expansion of Eq. (3) at the point (xe , e , e ), in the direction x, , is given by: t
= (J0 (x)t J0 (x))
−1
J0 (x)t b
(10)
where
t
t
J(xe )t e + J(x) e + J(xe )t + J(x)
J(xe + x) (e + ) = −[1 + 2(e + )rt (e + )] =
− (1 + 2e rt e ) − 2rt e − 2e rt − 2rt
−g(xe + x, e + ) =
− g(xe , e ) + T0 xe x − 2e r
(5)
1 + xt T0 x − 2 r 2 where, from the linear nature of the Jacobian matrix stated by Eq. (4), t
t e T0 x
t
J(x) e = x T0 e =
t
= J(e ) x = G0 x
J0 (x) =
and G0 = T0 e is a symmetric matrix. Eq. (5) can be represented in a compact form as, F(ze + d) = F(ze ) + W(ze )d + g1 (d) where zte = [xte
⎡
e t
J(xe ) e
t
e ], dt = [xt
⎤
F(ze ) = ⎣ −(1 + 2e rt e ) ⎦ ,
⎡ ⎢ ⎣
g1 (d) = ⎢
t ],
⎡
G0
W(ze ) = ⎣ 0
−g(xe , e )
and
(6)
J(xe )
t
J(x) −2rt
t
0
J(xe )
−2rt e
−2e rt
−2e r
0
⎤ ⎦ (7)
1 t x T0 x − 2 r 2 where g1 (d) is the so-called tensor term [19].
(8)
J(x)t
b=
2rt
0
−1
The computational difficulties inherent to Eq. (10) are well known, particularly for the term involving the product (J0 (x)t J0 (x)). The present study proposes two alternative strategies to estimate the initial values of the Lagrange multipliers. In the first, the following linear system is solved:
J(x)t 2rt
⎤ ⎥ ⎥ ⎦
2r 0
=b
(11)
where b is the same previously defined, J(x) is calculated for the initial estimate of the power flow variables and can be interpreted as a relaxation factor applied to the first optimality condition stated by Eq. (3). Note that if the Jacobian matrix is singular, the solution of this linear system provides = 0.
R.S. Salgado, A.F. Zeitune / Electric Power Systems Research 103 (2013) 114–121
The second strategy is a variation of the technique presented in [14] and consists of solving the following least squares problem: 1 t ss 2
Minimize
(1 + 2rt )
(12)
=0
where the equality constraints are the first two optimality conditions given by Eq. (3) and s is a kind of residual of the first equality constraint. The solution of this problem is obtained by solving the linear system given by:
⎡
0
⎢ J(x)t ⎣
J(x)
−2r
−I
0
−2rt
0
The effective search direction is defined with the scalar ˇ (0 ≤ ˇ ≤ 1) that minimises the function: (ˇ) =
− J(x)t + s = 0
subject to
⎤⎡
⎤
⎡ ⎤ 0
⎥ ⎣ 1 ⎦ ⎣ ⎦ = 0 ⎦ 2
0
(13)
1
where 1 and 2 are the Lagrange multipliers corresponding to the equality constraints, and I is the identity matrix of adequate size. This can be seen as an alternative to obtain the solution stated by Eq. (10), which has the advantage of preventing the difficulties previously mentioned. Note that in both cases only the solution of one linear system is required. 3.2. Use of tensor components With the aim of using the second-order Taylor series information, a component based on the tensor term of Eq. (6) is added to the increment obtained via Newton’s method. The calculation of this extra component consists of solving the linear system: W(ze )ds = −g1 (dn )
d = dn + ˇds
1 F(ze + d)t F(ze + d) 2
(18)
which basically requires the determination of the roots of a cubic equation, as shown in [18]. Observe that the basic direction of search is defined by Newton algorithm. The tensor component is used only if it improves the step in terms of the optimality conditions. Note also, that the additional computational effort required for this purpose corresponds to the calculation of the four last terms of Eq. (16). Previous works in the literature show that, in spite of these additional calculations, the use of the tensor term improves the convergence of the iterative process to solve both sets of nonlinear equations [22] and Sequential Quadratic Programming problems [23]. 3.3. Reactive power generation limits The reactive power generation limits are generically represented by the inequality hqg (x, ) ≤ hlim qg
(19)
where the indices qg and lim denote the reactive power generation and its limit, respectively. From Eq. (2), the generation of reactive power is given by: hqg (x, ) = yq0 + 2 rqg − h0 (x)
(20)
with: h0 (x) =
(14)
where the subscript s denotes the component of the search direction based on the tensor term and, g1 (dn ) is computed with the aid of Eq. (8). The composed search direction is defined as
117
1 t x T1 x 2
expansion where T1 is a three-dimensional array. The Taylor series of hqg (x, ) at the point (xe , e ), in the direction ˛x, ˛ is given by: hqg (xe + ˛x, e + ˛)
=
hqg (xe , e ) + ˛[−Jqg (xe )x + 2e rqg ] (21)
(15)
− ˛2 (h0 (x) + 2 rqg )
where the directions dn and ds are obtained by solving the linear systems of Eqs. (9) and (14) and ˇ is a scalar (0 ≤ ˇ ≤ 1). The function representing the optimality conditions is given by
where Jqg (x) is the Jacobian matrix of h0 (x) and ˛ is a scalar. The combination of Eqs. (19) and (21) provides:
F(ze + d) = F(ze ) + W(ze )(dn + ˇds ) + g1 (dn + ˇds )
aq + ˛bq + ˛2 cq ≤ 0
such that the appropriate arrangement of the terms of the last equation provides,
where,
F(ze + d)
=
F(ze ) + W(ze )dn + g1 (dn )
+ˇ2 g where
⎡
(16)
1 (ds )
t
t
J(xn ) s + J(xs ) n
2g0 (xn , xs ) − 2n s r t
and g0 (xn , xs ) = 12 xn T0 xs . Eq. (16) can be expressed in the compact form, F(ze + d) = a + ˇb + ˇ2 c = 0 where, b = −g1 (dn ) + g2 (dn , ds )
c = g1 (ds )
(x)+2 r
bq = −(Jqg (xe )x + 2e rqg )
(23)
qg )
For each reactive power generation device there is a quadratic equation of the form aqi + ˛bqi + ˛2 cqi ≤ 0
⎤
⎢ ⎥ t t g2 (dn , ds ) = ⎢ −2n r s − 2s r n ⎥ ⎣ ⎦
a = g1 (dn )
aq = hqg (xe , e )−hlim qg cq = −(h0
+ˇ(W(ze )ds + g2 (dn , ds ))
(22)
(17)
which is a component of Eq. (22). At each iteration this set of quadratic equations is solved and the minimum positive value of ˛ (0 ≤ ˛ ≤ 1), denoted ˛* , is selected as a step factor to adjust the increments ei , fi , i , i . From this point onward, this bus is converted into a PQ bus. Since the power flow solution of the base case is supposed to satisfy the reactive power generation limits, this procedure ensures that at each iteration only one bus eventually reaches the reactive power generation limit and that no limit of the reactive power generation will be violated. This way of handling the inequalities is equivalent to building a subset of inequality constraints, selecting one at each iteration, to be the so-called working set. These inequalities are imposed as equalities whereas all other inequalities are ignored, such that an
R.S. Salgado, A.F. Zeitune / Electric Power Systems Research 103 (2013) 114–121
equality constrained optimisation subproblem is solved at each iteration. In case of the power flow equations expressed in rectangular coordinates, the modification in the working set consists only of an interchange between the squared voltage magnitude equations and the reactive power balance equations of the PV buses. Thus the size of the set of the equality constraints does not change during the iterative process. In case of large systems, this strategy is computationally less expensive than solving the optimisation problem considering all constraints simultaneously [21]. Here, we take the advantage of using the information provided by the tensor term to compute more accurate estimates of the reactive power changes, which provides a better working set of inequalities. In the end of the iterative process, the sign of the Lagrange multipliers corresponding to the active inequality constraints is checked with respect to the Karush–Kuhn–Tucker conditions. If these conditions are not satisfied (the Lagrange multiplier of one inequality with the wrong sign), the corresponding constraint is excluded from the working set and the iterative process continues. Due to the characteristics of the maximum loadability optimisation problem, in the vast majority of the cases no constraint exclusion is needed, such that this procedure has been very effective in practice. 3.4. Algorithm
Summary of the iterative process 1 With base case PF,only Newton With base case PF, with Tensor 0.9
0.8
0.7
0.6 max |F(z)|
118
0.5
0.4
0.3
0.2
0.1
0
2
4
6
8
10
12
14
16
Iterations
Fig. 1. 1916-bus power system – convergence characteristics (a).
1. Specify initial values for the optimisation variables: (e0 ; f0 ; 0 ; 0 ); 2. Check convergence of the iterative process: does the component of F(z) with the largest magnitude satisfy a pre-specified tolerance? Do the signs of the Lagrange multipliers corresponding to inequality constraints satisfies the optimality conditions? If so, finish the iterative process. Otherwise, proceed to next step; 3. Solve linear systems (9) and (14) to obtain the search direction components dn and ds , and determine the scalar ˇ, minimizing the function represented in Eq. (18); 4. Calculate ˛ to limit the reactive power generation, according to Eq. (23), and make the PV–PQ bus conversion, if necessary; 5. Update the optimisation variables and return to step 2;
and qi = −1.0. The critical power flow solution resulting from the linear parameterisation corresponds to * = −132.62, with the bus voltage magnitudes ranging from 0.954 pu and 1.259 pu, which is clearly a useless power flow solution. This result is attributed to both the choice of the initial Lagrange multiplier estimates and the type of parameterisation, since a satisfactory solution is obtained with positive Lagrange multipliers. Subsequently, the maximum loadability power flow solution was obtained under the conditions stated previously, but with the power load variation modelled through quadratic parameterisation. In this case, the critical load parameter is = 8.220 (corresponding to an increase of 8 . 222 = 67.568% in the power load) and the bus voltage magnitude in the range 0.905 pu to 1.100 pu. Similar results are obtained even if all initial Lagrange multipliers estimates are equal to −1.0, indicating that the quadratic parameterisation effectively prevents negative load variation regardless the choice of the initial estimates of the Lagrange multipliers.
4. Simulation results
4.2. Effect of the tensor term
The proposed methodology was implemented in MatLab, version R2009a. The tests were focused on two main aspects. The first was the accuracy, robustness and computational effort of the iterative process. The second concerns the use of the sensitivity relationships both to identify critical regions and to define control actions at the maximum loadability point. The main features of the optimisation process were observed with respect to those of Newton’s method. The reason for this choice is that most of the direct methods use basically Newton strategy to define the search direction, as previously mentioned. In all cases the tolerance for convergence was adopted 1.0 × 10−8 and the rate of load increase was 1.0% of the base load, with a constant power factor. Although the active power generation could be adjusted according to Eq. (1), the slack bus was assumed to compensate the load increase.
In order to assess the effectiveness of the tensor term, the operational constraints were ignored and the initial estimates of the Lagrange multipliers were assumed equal to the unity. The maximum component (in magnitude) of the vector function representing the optimality conditions (Eq. (3)) was used to check the convergence of the iterative process. The results obtained with the 1916-bus system, which is an equivalent network of the Southern Region of Brazil, are summarised as follows. In all four cases analysed, the critical value of the load parameter at the convergence was * = 4.157. Figs. 1 and 2 show the trajectories from the initial estimates to the critical loadability solution. A base case power flow solution (Fig. 1) and the flat start profile (Fig. 2) were used as initial estimates of the complex bus voltages. The trajectories correspond to the use of Newton’s method (only) and the inclusion of the tensor component. The analysis of Fig. 1 indicates that the use of the second order term provides a decrease from 17 (Newton) to 10 (tensor) in the number of iterations to reach the optimal solution, with a corresponding reduction in the CPU time from 49.36 (Newton) s to 36.93 (tensor) s. Note that a power flow solution is available in the first iteration and thus the figure of merit used to check the convergence is relatively small in the first iteration. The peak in the second iteration is a consequence of the adjustment of the Lagrange multipliers (all equal to unity
The procedure for applying the proposed methodology can be summarised in the following steps:
4.1. Effect of the parameterisation In order to illustrate the effect of the parameterisation on the search for the critical loadability solution, we summarise the numerical results obtained with the IEEE 24-bus test system. For this purpose, the optimisation problem was solved through the conventional Newton’s method. The initial estimates were obtained from a base case power flow and multipliers pi = 1.0
R.S. Salgado, A.F. Zeitune / Electric Power Systems Research 103 (2013) 114–121
119
Table 1 Number of iterations required for convergence: quadratic parameterisation. System
0
IEEE24-bus IEEE57-bus IEEE118-bus IEEE300-bus Sys749-bus Sys1916-bus
1 T
N
T
N
T
NQg
16 14 37 30 20 28
11(4) 12(5) 35(26) 27(12) 18(7) 23(8)
13 13 38 32 17 32
11(7) 12(7) 35(8) 28(9) 16(4) 32(4)
16 15 46 33 20 35
13(4) 13(6) 44(5) 30(5) 18(10) 28(8)
7 6 33 22 9 19
Summary of the iterative process without base case PF,only Newton without base case PF, with Tensor 45
40
35
max |F(z)|
30
25
20
15
10
5
2
4
6
i
N
50
0
Qgi = Qglim
2
8
10
12 Iterations
14
16
18
20
22
Fig. 2. 1916-bus power system – convergence characteristics (b).
in the first iteration) to satisfy the first optimality condition. After iteration 2, the figure of merit decreases until the convergence is reached, in spite of the other small magnitude peaks at iterations 8 (tensor) and 12 (Newton). In case of the flat start (Fig. 2), the figure of merit is considerable larger in the beginning of the iterative process (the flat start and the unit Lagrange multipliers are more distant from the maximum loadability power flow solution), with a similar peak in the second iteration. The number of iterations required for determining these solutions were 23 (Newton) and 10 (tensor) with respective CPU times of 70.40 s (Newton) and 36.71 s (Tensor). Although the considerable decrease in the in the number of iterations, which can be attributed to the use of the tensor term, some additional computation is needed (the computations of the tensor term and the solution of the second linear system), reason by which the CPU time has not decreased proportionally. 4.3. Effect of the initial estimates The objective of this test was to verify how different initial estimates affect the robustness of the iterative process. The initial
estimate of the complex bus voltages was obtained from a power flow solution for the base case and three types of Lagrange multipliers were tested: 0 (all multipliers equal to unity), 1 (multipliers calculated through Eq. (11), and 2 (multipliers calculated through Eq. (13)). The critical loadability solution was determined for power networks with sizes ranging from 24 to 1916 buses. The convergence of the iterative process was observed with respect to the use of both Newton’s direction only (denoted N) and its combination with the tensor term (denoted T). In all cases, the reactive power generation limits were taken into account. Table 1 shows the main features of the iterative process in terms of the number of iterations required for convergence. In case of the tensor strategy, the number in parenthesis indicates how many times the tensor term was included in the search direction defined by Eq. (15). The last column of this table shows the number of buses whose reactive power generation limit was reached. The analysis of Table 1 reveals that the use of the second-order information on the tensor term generally reduces the number of iterations required for convergence (columns 3, 5 and 7) relative to the use of the Newton direction only (columns 2, 4 and 6). It must be pointed out that the region of feasible power flow solutions tends to become smaller with the inclusion of reactive power generation limits. For this reason, the use of the tensor term tends to be less effective during the iterations in which some inequality constraint becomes active. The results obtained for the 1916-bus system corresponding to the case without (Fig. 1) and with (Table 1) inequality constraints illustrate this remark. It is also important to observe that the fold points determined from distinct initial guesses have no significative differences. Although the strategies to calculate 1 and 2 are relatively simple, this indicates the robustness of the iterative process with respect to the initial solution, and no need of additional strategies to chose the initial Lagrange multipliers. Aiming at observing the accuracy and the computational effort of the results obtained via the proposed method, the fold points determined through Newton’s method, Continuation method and Sequential Quadratic Programming (SQP) of the software Knitro (version 7.0, developed by Ziena Optimisation, Inc.) were taken as references. Table 2 shows the parameter at the critical loadability and the computational effort required for convergence in terms of processing time. Direct methods provide the power flow solution at the critical loadability in one step. For this reason they cannot be compared to the Continuation method in terms of cpu time. Note also that Knitro is a general constrained optimisation software.
Table 2 Load parameter and CPU time (in s) required for convergence. System
IEEE24-bus IEEE57-bus IEEE118-bus IEEE300-bus Sys749-bus Sys1916-bus
Proposed methodology
Newton’s method
Continuation method
Knitro
*
CPU (s)
*
CPU (s)
*
CPU (s)
*
CPU (s)
67.58 44.10 61.92 2.47 107.43 4.12
0.091 0.191 0.557 0.968 9.186 66.272
67.58 44.10 61.92 2.49 107.41 4.15
0.121 0.132 0.467 2.321 12.38 75.44
65.19 44.10 61.38 2.48 106.32 4.10
0.772 1.459 11.742 3.872 122.674 331.360
65.25 44.10 63.92 2.47 107.41 4.13
0.873 1.304 6.523 10.760 58.192 249.004
120
R.S. Salgado, A.F. Zeitune / Electric Power Systems Research 103 (2013) 114–121
Fig. 3. IEEE 118 bus test system.
Thus the results shown in Table 1 should be used only as a reference with respect to the accuracy of the critical loadability parameter. The cpu times presented here are aimed only at providing a general idea about the computational effort required by the proposed technique. They must not be taken for strict comparison purposes. The differences between the critical load parameters obtained through these four methods are attributed mainly to the way the inequality constraints are handled. In the proposed approach as well as in
Newton’s method, the treatment of the reactive power generation limits is based on the strategy described in Section 3.3. On the other hand, the continuation power flow uses the conventional strategy of conversion PV–PQ bus, for each power flow solution, as the load increases. The SQP algorithm used in Knitro uses a scheme based on the selection of the Active Set to manipulate the inequality constraints. These strategies are prone to provide different trajectories to distinct (although very close) critical loadability points.
Table 3 Voltage magnitude of the PV buses with Qgi = Qglim at the bifurcation point. i
Bus
1
6
12
15
18
19
32
34
36
46
49
Vi (pu)
(0)
0.955
0.990
0.990
0.970
0.973
0.963
0.964
0.986
0.980
1.005
1.025
Vi (pu)
(∗)
0.919
0.976
0.962
0.914
0.923
0.903
0.961
0.868
0.862
0.877
0.910
Bus
54
55
56
59
61
62
65
66
70
72
74
0.955
0.952
0.954
0.985
0.995
0.998
1.005
1.050
0.984
0.980
0.958
0.966
0.953
0.954
0.984
1.034
1.019
0.942
1.010
0.753
0.939
0.700
(0) Vi (pu) (∗) Vi (pu)
Bus
76
77
80
85
92
99
100
103
104
105
110
(0) Vi (pu) (∗) Vi (pu)
0.943
1.006
1.040
0.985
0.993
1.010
1.017
1.001
0.971
0.965
0.973
0.694
0.854
0.944
0.930
0.963
0.998
0.959
0.940
0.917
0.20
0.954
R.S. Salgado, A.F. Zeitune / Electric Power Systems Research 103 (2013) 114–121 Table 4 IEEE118-bus test system: Lagrange multipliers and tangent vector. Bus
p
Bus
q
Bus
tgnv
1 117 2 13 3 14 15 19 12 16 33 11 7 18 6
1.0000 0.9922 0.9876 0.9853 0.9778 0.9761 0.9676 0.9675 0.9661 0.9634 0.9604 0.9587 0.9572 0.9549 0.9492
44 45 47 43 48 51 50 52 58 53 57 38 63 67 64
1.0000 0.9836 0.8011 0.7471 0.7343 0.7120 0.7063 0.7056 0.6637 0.6524 0.6512 0.6194 0.5690 0.5550 0.5246
45 44 52 47 51 53 58 48 50 57 71 63 43 67 60
1.0000 0.9652 0.8898 0.8875 0.8818 0.8701 0.8642 0.8623 0.8482 0.8476 0.7014 0.6662 0.6623 0.6536 0.6433
121
number of inequality constraints become active during the optimisation process, and (3) there is no need of special procedures to calculate the initial guesses of the Lagrange multipliers. Current work in this subject consists of developing strategies: (1) to improve the selection of the working inequality constraints set, with the aim of increasing the influence of the tensor term in the search direction even if an inequality constraint becomes active and (2) to include other types of inequality constraints in the optimisation problem and identify limit induced bifurcation points. Acknowledgements The authors acknowledge the financial support from the Brazilian Government Funding Agency Conselho Nacional de Desenvolvimento Científico e Tecnológico, CNPq/Brazil. References
4.4. The use of Lagrange multipliers The system used to obtain the numerical results described in this subsection is the IEEE118-bus, whose line diagram is shown in Fig. 3. The total load of the base case is 3668.0 MW and 1438.0 Mvar. The loadability margin determined by the proposed approach was 61.922% (Pd∗ = 5939.29 MW; Qd∗ = 2328.44 Mvar). The minimum and maximum bus voltage magnitudes at the critical loadability point were found at buses 118 (Vmin = 0.691 pu) and 25 (Vmax = 1.050 pu). At this point, 33 of 54 generation buses encountered their maximum reactive power generation limit. Table 3 presents these buses with their voltage magnitude in the base case and at the maximum loadability point. The purpose of this test is to verify how the sensitivity relationships between the power system variables can be exploited to identify critical areas and to determine control actions. Several studies have pointed the tangent vector as a well-accepted sensitivity relationship, suitable to check the critical areas in voltage stability studies (see Refs. [7–9]), reason by which it was chosen for this purpose. Table 4 shows the buses with the largest Lagrange multipliers and components of the tangent vector. These values were normalised with respect to the corresponding maximum of each type. From columns 3 and 5 we can infer that qi and the tangent vector indicate similar sets of buses with the largest sensitivities (columns 3–6). The buses with the largest qi (or tgnv) compose the critical area in terms of voltage stability, which is indicated by the shaded part in Fig. 3. Observe also that six of the thirty-three voltage control buses which reached the reactive power limit (buses 34, 42, 46, 49, 54 and 56 in Table 3) are in the neighbourhood of the critical voltage stability area. 5. Conclusions This paper proposes the determination of critical loadability power flow solutions through direct methods. The main features of the approach presented here are: (a) the use of a scheme of quadratic parameterisation to model the load variation; (b) the use of rectangular coordinates to represent the power flow equations; (c) the exploration of the tensor term of the power flow equations; (d) the use of alternative strategies to obtain initial estimates of the Lagrange multipliers; (e) the use of a scheme based on the step factor to treat the reactive power generation limits. The numerical results obtained from the application of the proposed direct method show that: (1) undesirable fold points with negative loading factor can be avoided by using the quadratic parameterisation of the load increase, (2) the use of the tensor term improves the convergence of the iterative process, in spite of the additional calculations required, and although it is less effective when a large
[1] S. Ayasun, C. Nwankpa, H.G. Kwatny, Computation of singular and singularity induced bifurcation points of differential-algebraic power system model, IEEE Transactions on Circuits and Systems – Regular Papers 51 (2004) 1525–1538. [2] R.J. Avalos, C.A. Ca nizares, F. Milano, A.J. Cpnejo, Equivalency of continuation and optimization methods to determine saddle node and limit-induced bifurcations in power systems, IEEE Transactions on Circuits and Systems – Regular Papers 56 (2009) 210–232. [3] P. Deufhard, Newton Methods for Nonlinear Problems, Springer-Verlag, Berlin, Heildelberg, 2005. [4] T. Van Cutsem, C. Vournas, Voltage Stability of Electric Power Systems, Kluwer Academic Publishers, Boston, USA, 1998. [5] V. Ajjarapu, C. Christy, The continuation power flow: a tool for steady state voltage stability analysis, IEEE Transactions on Power Systems 10 (1991) 304–310. [6] K. Iba, H. Suzuki, M. Egawa, T. Watanabe, Calculation of critical loading condition with nose curve using homotopy Continuation method, IEEE Transactions on Power Systems 6 (1991) 584–593. [7] C.A. Ca nizares, F.L. Alvarado, Point of collapse and Continuation method for large AC/DC systems, IEEE Transactions on Power Systems 8 (1993) 1–8. [8] A.C.Z. de Souza, C.A. Ca nizares, V.H. Quintana, New techniques to speed up voltage collapse computations using tangent vectors, IEEE Transactions on Power Systems 12 (1997) 1380–1387. [9] T. Van Cutsem, A method to compute reactive power margins with respect to voltage collapse, IEEE Transactions on Power Systems 6 (1991) 145–156. [10] C. Ca nizares, F. Alvarado, C. DeMarco, I. Dobson, W. Long, Point of collapse methods applied to AC/DC power systems, IEEE Transactions on Power Systems 7 (1992) 673–683. [11] I. Dobson, L. Lu, Y. Hu, A direct method for computing a closest saddle node bifurcation in the load power parameter space of an electric power system, in: Proc. of the IEEE International Symposium on Circuits and Systems, vol. 5, 1991, pp. 3019–3022. [12] I. Dobson, L. Lu, New methods for computing a closest saddle node bifurcation and worst case load power margin for voltage collapse, IEEE Transactions on Power Systems 8 (1993) 905–913. [13] Z. Yan, F. Wu, Y. Ni, Method for direct calculation of quadratic turning points, in: Proc. of the IET Generation, Transmission and Distribution, vol. 151, 2004, pp. 83–89. [14] F.M. Echavarren, E. Lobato, L. Rouco, A power flow solvability identification and calculation algorithm., Electric Power System Research 76 (2006) 242–250. [15] X. Yang, X. Zhou, Z. Du, Efficient solution algorithms for computing fold points of power flow equations., Electric Power System Research 33 (2011) 229–235. [16] G.D. Irisarri, X. Wang, J. Tong, S. Mokhtari, Maximum loadability of power systems using non linear Interior Point Method, IEEE Transactions on Power Systems 12 (1997) 162–172. [17] Y. Dai, J.D. McCalley, V. Vittal, Simplification, expansion and enhancement of direct interior point algorithm for power system maximum loadability, IEEE Transactions on Power Systems 15 (2000) 1014–1021. [18] S. Iwamoto, Y. Tamura, A load flow calculation method for ill-conditioned power systems, IEEE Transactions on Power Apparatus and Systems 100 (1981) 1736–1743. [19] R.B. Schnabel, P.D. Frank, Tensor methods for nonlinear equations, SIAM Journal on Numerical Analysis 21 (1984) 815–843. [20] R.S. Salgado, A.F. Zeitune, A framework to study critical loadability solutions, in: Proc. of the IEEE PowerTech Conference, Trondheim-Norway, 2011, pp. 1–6. [21] J. Nocedal, S.J. Wright, Numerical Optimization, Springer, USA, 1999. [22] R.S. Salgado, A.F. Zeitune, Power flow solutions through tensor methods, in: Proc. of the IET Generation, Transmission and Distribution, vol. 3, 2009, pp. 413–424. [23] D. Feng, R.B. Schnabel, Tensor methods for equality constrained optimization, SIAM Journal on Numerical Analysis 6 (1996) 653–673.
Contents lists available at SciVerse ScienceDirect
Electric Power Systems Research journal homepage: www.elsevier.com/locate/epsr
A direct method based on tensor calculation to determine maximum loadability power flow solutions Roberto S. Salgado ∗ , Anesio F. Zeitune Department of Electrical Engineering, Power System Laboratory (Labspot), Federal University of Santa Catarina, Florianópolis, SC, Brazil
a r t i c l e
i n f o
Article history: Received 9 November 2012 Received in revised form 15 May 2013 Accepted 17 May 2013 Available online 15 June 2013 Keywords: Maximum loadability Direct methods Tensor term Saddle-node bifurcation Sensitivity analysis
a b s t r a c t In this paper, the determination of critical loadability points of the power flow equations is formulated as an optimisation problem. A quadratic parameterisation scheme is used to model the load change, which is equivalent to implicitly including a non-negativity constraint in the load variation. The power flow equations are expressed in rectangular coordinates. This is suitable to exploit the second order information (tensor term) of the equations representing the optimality conditions. The use of this term improves the convergence of the iterative process and facilitates the manipulation of the reactive power generation limits. Simulation results obtained with a number of power systems, including real networks, are used to illustrate the main features of the proposed methodology. © 2013 Elsevier B.V. All rights reserved.
1. Introduction Several aspects of voltage instability have been explained by the Bifurcation Theory of nonlinear systems. In particular, two types of local bifurcations have been frequently associated with the steady state aspects of the voltage collapse [1,2]. The first is the saddle-node bifurcation, in which the Jacobian matrix of the differential-algebraic equations representing the power system is singular [3,4]. The second, named limit-induced bifurcation, is related to changes in the stability condition, as a consequence of a particular variable reaching its limit. In the present work, the power flow solution corresponding to the maximum (here also referred to as critical) loadability is characterised as a saddle-node bifurcation point (also named limit point, fold point and turning point). Although very interesting from the practical point of view of power system stability, other types of bifurcation are not focused on here. A number of methods for determining the critical loadability point can be found in the literature. Continuation methods have been widely used to determine a sequence of solutions for the parameterised power flow equations from a base load to the point of maximum loadability [5–7]. A predictor-corrector strategy is applied to calculate the set of equilibrium points that compose the PV-curve of each bus. Additionally, the predictor tangent vector is obtained as a by-product of the computational process, with
∗ Corresponding author. E-mail address: [email protected] (R.S. Salgado). 0378-7796/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.epsr.2013.05.011
reduced computational effort through the use of the techniques presented in [8]. Direct methods determine the saddle-node bifurcation point in one step, providing simultaneously the right or left eigenvectors of the singular power flow Jacobian matrix. Some of the relevant approaches found in the literature are mentioned as follows. In [9], an optimisation problem is solved by Newton’s method, in which (without loss of generality) the active power level is constant, such that only the QV subproblem is solved. In [7,10], the set of equations representing the Transversality Conditions are extended to integrate HVDC links, area interchange power control and other power systems particularities. The resulting set of nonlinear equations is solved to provide the so-called Point of Collapse. Refs. [11,12] propose direct methods to compute the saddle-node bifurcation point closest to the base case in the load parameter space. These methods combine the computation of both the load power margin and the vector normal to the limit surface of the feasible region of the power flow equations. They can provide the saddle-node bifurcation point in a fixed direction of load increase as well as the bifurcation point locally closest to the current operating load level. Ref. [13] proposes the determination of the loadability margin by solving the Moore–Spence augmented system. Its major contribution is the application of a matrix decomposition strategy to speed up the computational process. Ref. [14] combines the Continuation method with optimisation techniques to detect the limit of solvability of the power flow equations. A least squares technique is used to calculate the Lagrange multipliers, which are used as a figure of merit to find the power flow solvability limit. More recently, Yang et al. [15] proposes the use of block elimination techniques to solve the Moore–Spence augmented system. Similarly
R.S. Salgado, A.F. Zeitune / Electric Power Systems Research 103 (2013) 114–121
to [13], the contribution of this approach is related to the implementation aspects. The similarities and differences between these methodologies concern five main points: (1) the analytical formulation, (2) the type of load parameterisation, (3) the solution method, (4) the initial estimates, and (5) the by-products of the computation. First, the sets of equations solved to calculate fold points are basically those which represent alternatively the Transversality Conditions [7,12] the Moore–Spence determining system [13,15], or the optimality conditions of the maximum loadability optimisation problem [9]. In case of the optimisation model, adequate assumptions must be taken into account in order to obtain fold points similar to those produced by Continuation methods. Except for this aspect there are no considerable differences between these models. Second, with respect to the parameterisation, all of these methodologies have used the linear power load parameterisation to model the load increase. Unless a non-negativity constraint on the load parameter is added to the formulation, there is a risk of determining undesirable critical loadability solutions with load decrease (negative load parameter), if the initial estimate is not well chosen. Referring to the solution method, due to the particular structure and features of the maximum loadability problem, all of these approaches apply Newton’s method to solve the set of nonlinear equations that provide the fold point. This method is well defined and quadratically (although only locally) convergent. It is efficient to deal with equality and inequality constraints [9], but requires a good initial estimate. For this purpose, strategies like iterations of the Inverse Power Method have been used in [7,10], information about the largest eigenvalue of the Jacobian matrix available in the current operating point is employed in [11,12] and the Continuation method is applied in [13–15] to provide the initial guesses for power flow variables and eigenvectors. From this point of view, the exception is the approach presented in [9], which claims no need for the solution of any separate system to obtain improved initial estimates. The numerical results obtained through the application of these direct methods include the fold point and the sensitivity relationships represented by the eigenvectors and Lagrange multipliers, which are used to define critical buses and areas. The maximum loadability problem can be also formulated as a general Optimal Power Flow problem, as in [16,17], for example. The analytical models used in these approaches include other types of operational constraints (voltage magnitude, power flows, etc.) and optimisation variables (generated active and reactive powers, transformer taps, etc.). This brings additional difficulties (in terms of size and complexity of the problem) to obtain critical solutions and requires different numerical methods to solve the optimisation problem (as Interior Point Methods and/or Trust Region based-methods, for example). The maximum loadability solutions produced through the application of these methodologies is completely different from those focused on the present work. For this reason this type of approach is out of the scope of this paper. The present paper proposes a methodology for the direct determination of the turning points of the power flow equations. This is formulated as a constrained static optimisation problem, which is solved by an extended version of Newton’s method. For this purpose, the power flow equations are expressed in rectangular coordinates. The main contributions reported here are: (a) the use of an alternative load parameterisation scheme to model the load variation, so as to ensure power flow solutions with load level greater than that of the base case; (b) the inclusion of the secondorder information (tensor component) in the search direction, to improve the robustness of the iterative process and facilitate the treatment of the reactive power generation limits. The proposed approach differs from those previously mentioned with respect to: (1) the type of load parameterisation, (2) the strategy to solve the optimisation problem, which is an extension of Newton’s method
115
to exploit the second-order information. Similarly to [9], there is no need of any special strategy to choose initial guesses for the Lagrange multipliers, although we also show two alternative procedures to compute the initial estimates of the Lagrange multipliers. Numerical results obtained for power systems with sizes ranging from 24 to 1916 buses are used to illustrate the proposed methodology. This paper is organised as follows. Section 2 describes the theoretical basis of the maximum loadability problem, including the solution through static optimisation algorithms. Section 3 presents the analytical model proposed herein, with emphasis on the use of the information conveyed by the tensor term. Section 4 presents the numerical results obtained from the application of the methodology proposed here and Section 5 summarises the main conclusions. 2. Preliminary concepts The determination of the power flow solution at the maximum loadability level can be formulated as an optimisation problem expressed in terms of power system variables as: Maximize
subject to
(Pg0 + 2 Pgi ) − (Pd0 + 2 Pdi ) − Pi (e, f ) = i
i
Qg0 − (Qd0 + 2 Qdi ) − Qi (e, f ) = i
ref
Vi
i
2
0 (PV, PQ buses)
0 (PQ buses)
(1)
= ei2 + fi2 (PV buses)
Qgi = (Qd0 + 2 Qdi ) + Qi (e, f ) ≤ Qglim (PV buses) i
Pg0i
Qg0i
i
i
i
where and are the active and reactive power generation of bus i in the base case and Pgi represents the active power generation change rate, which indicates how the active power generation follows the increase of the loadability level (distribution factors based on economic and/or security criteria can be used for this purpose); Pd0 and Qd0 refer to the active and reactive power demand of the ith bus in the base case, Pdi and Qdi represent the prespecified direction of increase of the active and reactive power load ref of bus i, Vi is the reference voltage magnitude of the ith PV bus, and Pi (e, f) and Qi (e, f) are the nodal active and reactive power injections expressed as functions of the real and imaginary components of complex voltages. The inequality constraints refer to limits of generated reactive power of the PV buses. The optimisation variables are the components of the complex voltage (ei , fi ) and the load parameter . If the conventional linear parameterisation is used to model the load variation and a non-negativity constraint is not imposed on , the optimisation problem of Eq. (1) can provide undesirable critical solutions (with negative loading factor, for example), which are useless from the practical point of view. This is usually due to a poor initial estimate of the optimisation variables, in particular the Lagrange multipliers. In order to overcome this difficulty, the load change is quadratically parameterised, as shown in Eq. (1). This type of parameterisation is equivalent to imposing a non-negativity constraint in the load variation. The constraints of the optimisation problem of Eq. (1) are expressed in compact form as: g(x, )
=
y0 + 2 r − g0 (x) = 0
hqg (x, )
=
yq0 + 2 rqg − h0 (x) ≤
hlim qg
(2)
where the components of vectors (y0 + 2 r) and (yq0 + 2 rqg ) are related to the bus power injections, the squared voltage magnitude of the PV buses in the base case (y0 and yq0 ) and the load increase direction (r and rqg ). Since the PV bus voltage magnitude does not depend on the load variation (see the third equality constraint in Eq. (1)), the components of the vector r corresponding to the squared voltage magnitudes are zero. Vectors g0 (x) and h0 (x)
116
R.S. Salgado, A.F. Zeitune / Electric Power Systems Research 103 (2013) 114–121
represent the power injections and the PV buses squared voltage magnitude expressed in terms of vector x, which stands for the real and imaginary components of the complex voltages, and hlim qg represents the reactive power generation limits. If only the equality constraints (the treatment of inequality constraints will be described later) are taken into account, the optimality conditions for the problem stated by Eq. (1) are:
∇ x £(x, , ) = 0 =
J(x)t
∇ £(x, , ) = 0 =
− 1 + 2rt
∇ £(x, , ) = 0 =
− (y0 + 2 r − g0 (x))
(3)
The most important aspect of Eq. (6) is that there are no terms with order higher than the second one, so that given a direction d the function F(ze + d) can be accurately computed.
2.2. Newton’s method Newton’s method uses a linear approximation of Eq. (6), such that at each iteration the solution of the following linear system is required: W(ze )dn = −F(ze )
where J(x) is the Jacobian matrix of equations g0 (x) with respect to the variables x. Since the rate of load increase (r) is non-zero, from the second optimality condition of Eq. (3), = / 0. Thus, equation J(x)t = 0 indicates that the Jacobian matrix is singular at the point of critical loadability, and that vector is co-linear with the left eigenvector corresponding to the zero eigenvalue of matrix J(x).
(9)
where the subscript n denotes the component of the search direction obtained by Newton’s method. The non-singularity of the coefficient matrix W(ze ) of Eq. (9) (formally shown in [3]) allows the point of critical loadability to be accurately determined. To speed up the iterative process the decomposition strategy of this matrix, shown in [13,20], can be applied.
2.1. Use of the quadratic form If the power flow equations are expressed in rectangular coordinates, the functions g0 (x) of Eq. (2) can be represented by a quadratic function expressed as [18] 1 g0 (x) = xt T0 x 2 where T0 is a three-dimensional array, called the tensor, which depends only on the transmission line parameters; such that its first and second derivatives are given by:
∂g0 (x) = J(x) = T0 x and ∂x
3. Proposed methodology 3.1. Initial condition In general, the initial estimates of the Lagrange multipliers are considerably far from those of the critical loadability solution, which can be a drawback to the application of Newton’s method. To specify initial values for these dual variables, the literature [21] suggests the least squares solution given by:
2
∂ g0 (x) ∂J(x) = = T0 ∂x2 ∂x
(4)
The second order Taylor series expansion of Eq. (3) at the point (xe , e , e ), in the direction x, , is given by: t
= (J0 (x)t J0 (x))
−1
J0 (x)t b
(10)
where
t
t
J(xe )t e + J(x) e + J(xe )t + J(x)
J(xe + x) (e + ) = −[1 + 2(e + )rt (e + )] =
− (1 + 2e rt e ) − 2rt e − 2e rt − 2rt
−g(xe + x, e + ) =
− g(xe , e ) + T0 xe x − 2e r
(5)
1 + xt T0 x − 2 r 2 where, from the linear nature of the Jacobian matrix stated by Eq. (4), t
t e T0 x
t
J(x) e = x T0 e =
t
= J(e ) x = G0 x
J0 (x) =
and G0 = T0 e is a symmetric matrix. Eq. (5) can be represented in a compact form as, F(ze + d) = F(ze ) + W(ze )d + g1 (d) where zte = [xte
⎡
e t
J(xe ) e
t
e ], dt = [xt
⎤
F(ze ) = ⎣ −(1 + 2e rt e ) ⎦ ,
⎡ ⎢ ⎣
g1 (d) = ⎢
t ],
⎡
G0
W(ze ) = ⎣ 0
−g(xe , e )
and
(6)
J(xe )
t
J(x) −2rt
t
0
J(xe )
−2rt e
−2e rt
−2e r
0
⎤ ⎦ (7)
1 t x T0 x − 2 r 2 where g1 (d) is the so-called tensor term [19].
(8)
J(x)t
b=
2rt
0
−1
The computational difficulties inherent to Eq. (10) are well known, particularly for the term involving the product (J0 (x)t J0 (x)). The present study proposes two alternative strategies to estimate the initial values of the Lagrange multipliers. In the first, the following linear system is solved:
J(x)t 2rt
⎤ ⎥ ⎥ ⎦
2r 0
=b
(11)
where b is the same previously defined, J(x) is calculated for the initial estimate of the power flow variables and can be interpreted as a relaxation factor applied to the first optimality condition stated by Eq. (3). Note that if the Jacobian matrix is singular, the solution of this linear system provides = 0.
R.S. Salgado, A.F. Zeitune / Electric Power Systems Research 103 (2013) 114–121
The second strategy is a variation of the technique presented in [14] and consists of solving the following least squares problem: 1 t ss 2
Minimize
(1 + 2rt )
(12)
=0
where the equality constraints are the first two optimality conditions given by Eq. (3) and s is a kind of residual of the first equality constraint. The solution of this problem is obtained by solving the linear system given by:
⎡
0
⎢ J(x)t ⎣
J(x)
−2r
−I
0
−2rt
0
The effective search direction is defined with the scalar ˇ (0 ≤ ˇ ≤ 1) that minimises the function: (ˇ) =
− J(x)t + s = 0
subject to
⎤⎡
⎤
⎡ ⎤ 0
⎥ ⎣ 1 ⎦ ⎣ ⎦ = 0 ⎦ 2
0
(13)
1
where 1 and 2 are the Lagrange multipliers corresponding to the equality constraints, and I is the identity matrix of adequate size. This can be seen as an alternative to obtain the solution stated by Eq. (10), which has the advantage of preventing the difficulties previously mentioned. Note that in both cases only the solution of one linear system is required. 3.2. Use of tensor components With the aim of using the second-order Taylor series information, a component based on the tensor term of Eq. (6) is added to the increment obtained via Newton’s method. The calculation of this extra component consists of solving the linear system: W(ze )ds = −g1 (dn )
d = dn + ˇds
1 F(ze + d)t F(ze + d) 2
(18)
which basically requires the determination of the roots of a cubic equation, as shown in [18]. Observe that the basic direction of search is defined by Newton algorithm. The tensor component is used only if it improves the step in terms of the optimality conditions. Note also, that the additional computational effort required for this purpose corresponds to the calculation of the four last terms of Eq. (16). Previous works in the literature show that, in spite of these additional calculations, the use of the tensor term improves the convergence of the iterative process to solve both sets of nonlinear equations [22] and Sequential Quadratic Programming problems [23]. 3.3. Reactive power generation limits The reactive power generation limits are generically represented by the inequality hqg (x, ) ≤ hlim qg
(19)
where the indices qg and lim denote the reactive power generation and its limit, respectively. From Eq. (2), the generation of reactive power is given by: hqg (x, ) = yq0 + 2 rqg − h0 (x)
(20)
with: h0 (x) =
(14)
where the subscript s denotes the component of the search direction based on the tensor term and, g1 (dn ) is computed with the aid of Eq. (8). The composed search direction is defined as
117
1 t x T1 x 2
expansion where T1 is a three-dimensional array. The Taylor series of hqg (x, ) at the point (xe , e ), in the direction ˛x, ˛ is given by: hqg (xe + ˛x, e + ˛)
=
hqg (xe , e ) + ˛[−Jqg (xe )x + 2e rqg ] (21)
(15)
− ˛2 (h0 (x) + 2 rqg )
where the directions dn and ds are obtained by solving the linear systems of Eqs. (9) and (14) and ˇ is a scalar (0 ≤ ˇ ≤ 1). The function representing the optimality conditions is given by
where Jqg (x) is the Jacobian matrix of h0 (x) and ˛ is a scalar. The combination of Eqs. (19) and (21) provides:
F(ze + d) = F(ze ) + W(ze )(dn + ˇds ) + g1 (dn + ˇds )
aq + ˛bq + ˛2 cq ≤ 0
such that the appropriate arrangement of the terms of the last equation provides,
where,
F(ze + d)
=
F(ze ) + W(ze )dn + g1 (dn )
+ˇ2 g where
⎡
(16)
1 (ds )
t
t
J(xn ) s + J(xs ) n
2g0 (xn , xs ) − 2n s r t
and g0 (xn , xs ) = 12 xn T0 xs . Eq. (16) can be expressed in the compact form, F(ze + d) = a + ˇb + ˇ2 c = 0 where, b = −g1 (dn ) + g2 (dn , ds )
c = g1 (ds )
(x)+2 r
bq = −(Jqg (xe )x + 2e rqg )
(23)
qg )
For each reactive power generation device there is a quadratic equation of the form aqi + ˛bqi + ˛2 cqi ≤ 0
⎤
⎢ ⎥ t t g2 (dn , ds ) = ⎢ −2n r s − 2s r n ⎥ ⎣ ⎦
a = g1 (dn )
aq = hqg (xe , e )−hlim qg cq = −(h0
+ˇ(W(ze )ds + g2 (dn , ds ))
(22)
(17)
which is a component of Eq. (22). At each iteration this set of quadratic equations is solved and the minimum positive value of ˛ (0 ≤ ˛ ≤ 1), denoted ˛* , is selected as a step factor to adjust the increments ei , fi , i , i . From this point onward, this bus is converted into a PQ bus. Since the power flow solution of the base case is supposed to satisfy the reactive power generation limits, this procedure ensures that at each iteration only one bus eventually reaches the reactive power generation limit and that no limit of the reactive power generation will be violated. This way of handling the inequalities is equivalent to building a subset of inequality constraints, selecting one at each iteration, to be the so-called working set. These inequalities are imposed as equalities whereas all other inequalities are ignored, such that an
R.S. Salgado, A.F. Zeitune / Electric Power Systems Research 103 (2013) 114–121
equality constrained optimisation subproblem is solved at each iteration. In case of the power flow equations expressed in rectangular coordinates, the modification in the working set consists only of an interchange between the squared voltage magnitude equations and the reactive power balance equations of the PV buses. Thus the size of the set of the equality constraints does not change during the iterative process. In case of large systems, this strategy is computationally less expensive than solving the optimisation problem considering all constraints simultaneously [21]. Here, we take the advantage of using the information provided by the tensor term to compute more accurate estimates of the reactive power changes, which provides a better working set of inequalities. In the end of the iterative process, the sign of the Lagrange multipliers corresponding to the active inequality constraints is checked with respect to the Karush–Kuhn–Tucker conditions. If these conditions are not satisfied (the Lagrange multiplier of one inequality with the wrong sign), the corresponding constraint is excluded from the working set and the iterative process continues. Due to the characteristics of the maximum loadability optimisation problem, in the vast majority of the cases no constraint exclusion is needed, such that this procedure has been very effective in practice. 3.4. Algorithm
Summary of the iterative process 1 With base case PF,only Newton With base case PF, with Tensor 0.9
0.8
0.7
0.6 max |F(z)|
118
0.5
0.4
0.3
0.2
0.1
0
2
4
6
8
10
12
14
16
Iterations
Fig. 1. 1916-bus power system – convergence characteristics (a).
1. Specify initial values for the optimisation variables: (e0 ; f0 ; 0 ; 0 ); 2. Check convergence of the iterative process: does the component of F(z) with the largest magnitude satisfy a pre-specified tolerance? Do the signs of the Lagrange multipliers corresponding to inequality constraints satisfies the optimality conditions? If so, finish the iterative process. Otherwise, proceed to next step; 3. Solve linear systems (9) and (14) to obtain the search direction components dn and ds , and determine the scalar ˇ, minimizing the function represented in Eq. (18); 4. Calculate ˛ to limit the reactive power generation, according to Eq. (23), and make the PV–PQ bus conversion, if necessary; 5. Update the optimisation variables and return to step 2;
and qi = −1.0. The critical power flow solution resulting from the linear parameterisation corresponds to * = −132.62, with the bus voltage magnitudes ranging from 0.954 pu and 1.259 pu, which is clearly a useless power flow solution. This result is attributed to both the choice of the initial Lagrange multiplier estimates and the type of parameterisation, since a satisfactory solution is obtained with positive Lagrange multipliers. Subsequently, the maximum loadability power flow solution was obtained under the conditions stated previously, but with the power load variation modelled through quadratic parameterisation. In this case, the critical load parameter is = 8.220 (corresponding to an increase of 8 . 222 = 67.568% in the power load) and the bus voltage magnitude in the range 0.905 pu to 1.100 pu. Similar results are obtained even if all initial Lagrange multipliers estimates are equal to −1.0, indicating that the quadratic parameterisation effectively prevents negative load variation regardless the choice of the initial estimates of the Lagrange multipliers.
4. Simulation results
4.2. Effect of the tensor term
The proposed methodology was implemented in MatLab, version R2009a. The tests were focused on two main aspects. The first was the accuracy, robustness and computational effort of the iterative process. The second concerns the use of the sensitivity relationships both to identify critical regions and to define control actions at the maximum loadability point. The main features of the optimisation process were observed with respect to those of Newton’s method. The reason for this choice is that most of the direct methods use basically Newton strategy to define the search direction, as previously mentioned. In all cases the tolerance for convergence was adopted 1.0 × 10−8 and the rate of load increase was 1.0% of the base load, with a constant power factor. Although the active power generation could be adjusted according to Eq. (1), the slack bus was assumed to compensate the load increase.
In order to assess the effectiveness of the tensor term, the operational constraints were ignored and the initial estimates of the Lagrange multipliers were assumed equal to the unity. The maximum component (in magnitude) of the vector function representing the optimality conditions (Eq. (3)) was used to check the convergence of the iterative process. The results obtained with the 1916-bus system, which is an equivalent network of the Southern Region of Brazil, are summarised as follows. In all four cases analysed, the critical value of the load parameter at the convergence was * = 4.157. Figs. 1 and 2 show the trajectories from the initial estimates to the critical loadability solution. A base case power flow solution (Fig. 1) and the flat start profile (Fig. 2) were used as initial estimates of the complex bus voltages. The trajectories correspond to the use of Newton’s method (only) and the inclusion of the tensor component. The analysis of Fig. 1 indicates that the use of the second order term provides a decrease from 17 (Newton) to 10 (tensor) in the number of iterations to reach the optimal solution, with a corresponding reduction in the CPU time from 49.36 (Newton) s to 36.93 (tensor) s. Note that a power flow solution is available in the first iteration and thus the figure of merit used to check the convergence is relatively small in the first iteration. The peak in the second iteration is a consequence of the adjustment of the Lagrange multipliers (all equal to unity
The procedure for applying the proposed methodology can be summarised in the following steps:
4.1. Effect of the parameterisation In order to illustrate the effect of the parameterisation on the search for the critical loadability solution, we summarise the numerical results obtained with the IEEE 24-bus test system. For this purpose, the optimisation problem was solved through the conventional Newton’s method. The initial estimates were obtained from a base case power flow and multipliers pi = 1.0
R.S. Salgado, A.F. Zeitune / Electric Power Systems Research 103 (2013) 114–121
119
Table 1 Number of iterations required for convergence: quadratic parameterisation. System
0
IEEE24-bus IEEE57-bus IEEE118-bus IEEE300-bus Sys749-bus Sys1916-bus
1 T
N
T
N
T
NQg
16 14 37 30 20 28
11(4) 12(5) 35(26) 27(12) 18(7) 23(8)
13 13 38 32 17 32
11(7) 12(7) 35(8) 28(9) 16(4) 32(4)
16 15 46 33 20 35
13(4) 13(6) 44(5) 30(5) 18(10) 28(8)
7 6 33 22 9 19
Summary of the iterative process without base case PF,only Newton without base case PF, with Tensor 45
40
35
max |F(z)|
30
25
20
15
10
5
2
4
6
i
N
50
0
Qgi = Qglim
2
8
10
12 Iterations
14
16
18
20
22
Fig. 2. 1916-bus power system – convergence characteristics (b).
in the first iteration) to satisfy the first optimality condition. After iteration 2, the figure of merit decreases until the convergence is reached, in spite of the other small magnitude peaks at iterations 8 (tensor) and 12 (Newton). In case of the flat start (Fig. 2), the figure of merit is considerable larger in the beginning of the iterative process (the flat start and the unit Lagrange multipliers are more distant from the maximum loadability power flow solution), with a similar peak in the second iteration. The number of iterations required for determining these solutions were 23 (Newton) and 10 (tensor) with respective CPU times of 70.40 s (Newton) and 36.71 s (Tensor). Although the considerable decrease in the in the number of iterations, which can be attributed to the use of the tensor term, some additional computation is needed (the computations of the tensor term and the solution of the second linear system), reason by which the CPU time has not decreased proportionally. 4.3. Effect of the initial estimates The objective of this test was to verify how different initial estimates affect the robustness of the iterative process. The initial
estimate of the complex bus voltages was obtained from a power flow solution for the base case and three types of Lagrange multipliers were tested: 0 (all multipliers equal to unity), 1 (multipliers calculated through Eq. (11), and 2 (multipliers calculated through Eq. (13)). The critical loadability solution was determined for power networks with sizes ranging from 24 to 1916 buses. The convergence of the iterative process was observed with respect to the use of both Newton’s direction only (denoted N) and its combination with the tensor term (denoted T). In all cases, the reactive power generation limits were taken into account. Table 1 shows the main features of the iterative process in terms of the number of iterations required for convergence. In case of the tensor strategy, the number in parenthesis indicates how many times the tensor term was included in the search direction defined by Eq. (15). The last column of this table shows the number of buses whose reactive power generation limit was reached. The analysis of Table 1 reveals that the use of the second-order information on the tensor term generally reduces the number of iterations required for convergence (columns 3, 5 and 7) relative to the use of the Newton direction only (columns 2, 4 and 6). It must be pointed out that the region of feasible power flow solutions tends to become smaller with the inclusion of reactive power generation limits. For this reason, the use of the tensor term tends to be less effective during the iterations in which some inequality constraint becomes active. The results obtained for the 1916-bus system corresponding to the case without (Fig. 1) and with (Table 1) inequality constraints illustrate this remark. It is also important to observe that the fold points determined from distinct initial guesses have no significative differences. Although the strategies to calculate 1 and 2 are relatively simple, this indicates the robustness of the iterative process with respect to the initial solution, and no need of additional strategies to chose the initial Lagrange multipliers. Aiming at observing the accuracy and the computational effort of the results obtained via the proposed method, the fold points determined through Newton’s method, Continuation method and Sequential Quadratic Programming (SQP) of the software Knitro (version 7.0, developed by Ziena Optimisation, Inc.) were taken as references. Table 2 shows the parameter at the critical loadability and the computational effort required for convergence in terms of processing time. Direct methods provide the power flow solution at the critical loadability in one step. For this reason they cannot be compared to the Continuation method in terms of cpu time. Note also that Knitro is a general constrained optimisation software.
Table 2 Load parameter and CPU time (in s) required for convergence. System
IEEE24-bus IEEE57-bus IEEE118-bus IEEE300-bus Sys749-bus Sys1916-bus
Proposed methodology
Newton’s method
Continuation method
Knitro
*
CPU (s)
*
CPU (s)
*
CPU (s)
*
CPU (s)
67.58 44.10 61.92 2.47 107.43 4.12
0.091 0.191 0.557 0.968 9.186 66.272
67.58 44.10 61.92 2.49 107.41 4.15
0.121 0.132 0.467 2.321 12.38 75.44
65.19 44.10 61.38 2.48 106.32 4.10
0.772 1.459 11.742 3.872 122.674 331.360
65.25 44.10 63.92 2.47 107.41 4.13
0.873 1.304 6.523 10.760 58.192 249.004
120
R.S. Salgado, A.F. Zeitune / Electric Power Systems Research 103 (2013) 114–121
Fig. 3. IEEE 118 bus test system.
Thus the results shown in Table 1 should be used only as a reference with respect to the accuracy of the critical loadability parameter. The cpu times presented here are aimed only at providing a general idea about the computational effort required by the proposed technique. They must not be taken for strict comparison purposes. The differences between the critical load parameters obtained through these four methods are attributed mainly to the way the inequality constraints are handled. In the proposed approach as well as in
Newton’s method, the treatment of the reactive power generation limits is based on the strategy described in Section 3.3. On the other hand, the continuation power flow uses the conventional strategy of conversion PV–PQ bus, for each power flow solution, as the load increases. The SQP algorithm used in Knitro uses a scheme based on the selection of the Active Set to manipulate the inequality constraints. These strategies are prone to provide different trajectories to distinct (although very close) critical loadability points.
Table 3 Voltage magnitude of the PV buses with Qgi = Qglim at the bifurcation point. i
Bus
1
6
12
15
18
19
32
34
36
46
49
Vi (pu)
(0)
0.955
0.990
0.990
0.970
0.973
0.963
0.964
0.986
0.980
1.005
1.025
Vi (pu)
(∗)
0.919
0.976
0.962
0.914
0.923
0.903
0.961
0.868
0.862
0.877
0.910
Bus
54
55
56
59
61
62
65
66
70
72
74
0.955
0.952
0.954
0.985
0.995
0.998
1.005
1.050
0.984
0.980
0.958
0.966
0.953
0.954
0.984
1.034
1.019
0.942
1.010
0.753
0.939
0.700
(0) Vi (pu) (∗) Vi (pu)
Bus
76
77
80
85
92
99
100
103
104
105
110
(0) Vi (pu) (∗) Vi (pu)
0.943
1.006
1.040
0.985
0.993
1.010
1.017
1.001
0.971
0.965
0.973
0.694
0.854
0.944
0.930
0.963
0.998
0.959
0.940
0.917
0.20
0.954
R.S. Salgado, A.F. Zeitune / Electric Power Systems Research 103 (2013) 114–121 Table 4 IEEE118-bus test system: Lagrange multipliers and tangent vector. Bus
p
Bus
q
Bus
tgnv
1 117 2 13 3 14 15 19 12 16 33 11 7 18 6
1.0000 0.9922 0.9876 0.9853 0.9778 0.9761 0.9676 0.9675 0.9661 0.9634 0.9604 0.9587 0.9572 0.9549 0.9492
44 45 47 43 48 51 50 52 58 53 57 38 63 67 64
1.0000 0.9836 0.8011 0.7471 0.7343 0.7120 0.7063 0.7056 0.6637 0.6524 0.6512 0.6194 0.5690 0.5550 0.5246
45 44 52 47 51 53 58 48 50 57 71 63 43 67 60
1.0000 0.9652 0.8898 0.8875 0.8818 0.8701 0.8642 0.8623 0.8482 0.8476 0.7014 0.6662 0.6623 0.6536 0.6433
121
number of inequality constraints become active during the optimisation process, and (3) there is no need of special procedures to calculate the initial guesses of the Lagrange multipliers. Current work in this subject consists of developing strategies: (1) to improve the selection of the working inequality constraints set, with the aim of increasing the influence of the tensor term in the search direction even if an inequality constraint becomes active and (2) to include other types of inequality constraints in the optimisation problem and identify limit induced bifurcation points. Acknowledgements The authors acknowledge the financial support from the Brazilian Government Funding Agency Conselho Nacional de Desenvolvimento Científico e Tecnológico, CNPq/Brazil. References
4.4. The use of Lagrange multipliers The system used to obtain the numerical results described in this subsection is the IEEE118-bus, whose line diagram is shown in Fig. 3. The total load of the base case is 3668.0 MW and 1438.0 Mvar. The loadability margin determined by the proposed approach was 61.922% (Pd∗ = 5939.29 MW; Qd∗ = 2328.44 Mvar). The minimum and maximum bus voltage magnitudes at the critical loadability point were found at buses 118 (Vmin = 0.691 pu) and 25 (Vmax = 1.050 pu). At this point, 33 of 54 generation buses encountered their maximum reactive power generation limit. Table 3 presents these buses with their voltage magnitude in the base case and at the maximum loadability point. The purpose of this test is to verify how the sensitivity relationships between the power system variables can be exploited to identify critical areas and to determine control actions. Several studies have pointed the tangent vector as a well-accepted sensitivity relationship, suitable to check the critical areas in voltage stability studies (see Refs. [7–9]), reason by which it was chosen for this purpose. Table 4 shows the buses with the largest Lagrange multipliers and components of the tangent vector. These values were normalised with respect to the corresponding maximum of each type. From columns 3 and 5 we can infer that qi and the tangent vector indicate similar sets of buses with the largest sensitivities (columns 3–6). The buses with the largest qi (or tgnv) compose the critical area in terms of voltage stability, which is indicated by the shaded part in Fig. 3. Observe also that six of the thirty-three voltage control buses which reached the reactive power limit (buses 34, 42, 46, 49, 54 and 56 in Table 3) are in the neighbourhood of the critical voltage stability area. 5. Conclusions This paper proposes the determination of critical loadability power flow solutions through direct methods. The main features of the approach presented here are: (a) the use of a scheme of quadratic parameterisation to model the load variation; (b) the use of rectangular coordinates to represent the power flow equations; (c) the exploration of the tensor term of the power flow equations; (d) the use of alternative strategies to obtain initial estimates of the Lagrange multipliers; (e) the use of a scheme based on the step factor to treat the reactive power generation limits. The numerical results obtained from the application of the proposed direct method show that: (1) undesirable fold points with negative loading factor can be avoided by using the quadratic parameterisation of the load increase, (2) the use of the tensor term improves the convergence of the iterative process, in spite of the additional calculations required, and although it is less effective when a large
[1] S. Ayasun, C. Nwankpa, H.G. Kwatny, Computation of singular and singularity induced bifurcation points of differential-algebraic power system model, IEEE Transactions on Circuits and Systems – Regular Papers 51 (2004) 1525–1538. [2] R.J. Avalos, C.A. Ca nizares, F. Milano, A.J. Cpnejo, Equivalency of continuation and optimization methods to determine saddle node and limit-induced bifurcations in power systems, IEEE Transactions on Circuits and Systems – Regular Papers 56 (2009) 210–232. [3] P. Deufhard, Newton Methods for Nonlinear Problems, Springer-Verlag, Berlin, Heildelberg, 2005. [4] T. Van Cutsem, C. Vournas, Voltage Stability of Electric Power Systems, Kluwer Academic Publishers, Boston, USA, 1998. [5] V. Ajjarapu, C. Christy, The continuation power flow: a tool for steady state voltage stability analysis, IEEE Transactions on Power Systems 10 (1991) 304–310. [6] K. Iba, H. Suzuki, M. Egawa, T. Watanabe, Calculation of critical loading condition with nose curve using homotopy Continuation method, IEEE Transactions on Power Systems 6 (1991) 584–593. [7] C.A. Ca nizares, F.L. Alvarado, Point of collapse and Continuation method for large AC/DC systems, IEEE Transactions on Power Systems 8 (1993) 1–8. [8] A.C.Z. de Souza, C.A. Ca nizares, V.H. Quintana, New techniques to speed up voltage collapse computations using tangent vectors, IEEE Transactions on Power Systems 12 (1997) 1380–1387. [9] T. Van Cutsem, A method to compute reactive power margins with respect to voltage collapse, IEEE Transactions on Power Systems 6 (1991) 145–156. [10] C. Ca nizares, F. Alvarado, C. DeMarco, I. Dobson, W. Long, Point of collapse methods applied to AC/DC power systems, IEEE Transactions on Power Systems 7 (1992) 673–683. [11] I. Dobson, L. Lu, Y. Hu, A direct method for computing a closest saddle node bifurcation in the load power parameter space of an electric power system, in: Proc. of the IEEE International Symposium on Circuits and Systems, vol. 5, 1991, pp. 3019–3022. [12] I. Dobson, L. Lu, New methods for computing a closest saddle node bifurcation and worst case load power margin for voltage collapse, IEEE Transactions on Power Systems 8 (1993) 905–913. [13] Z. Yan, F. Wu, Y. Ni, Method for direct calculation of quadratic turning points, in: Proc. of the IET Generation, Transmission and Distribution, vol. 151, 2004, pp. 83–89. [14] F.M. Echavarren, E. Lobato, L. Rouco, A power flow solvability identification and calculation algorithm., Electric Power System Research 76 (2006) 242–250. [15] X. Yang, X. Zhou, Z. Du, Efficient solution algorithms for computing fold points of power flow equations., Electric Power System Research 33 (2011) 229–235. [16] G.D. Irisarri, X. Wang, J. Tong, S. Mokhtari, Maximum loadability of power systems using non linear Interior Point Method, IEEE Transactions on Power Systems 12 (1997) 162–172. [17] Y. Dai, J.D. McCalley, V. Vittal, Simplification, expansion and enhancement of direct interior point algorithm for power system maximum loadability, IEEE Transactions on Power Systems 15 (2000) 1014–1021. [18] S. Iwamoto, Y. Tamura, A load flow calculation method for ill-conditioned power systems, IEEE Transactions on Power Apparatus and Systems 100 (1981) 1736–1743. [19] R.B. Schnabel, P.D. Frank, Tensor methods for nonlinear equations, SIAM Journal on Numerical Analysis 21 (1984) 815–843. [20] R.S. Salgado, A.F. Zeitune, A framework to study critical loadability solutions, in: Proc. of the IEEE PowerTech Conference, Trondheim-Norway, 2011, pp. 1–6. [21] J. Nocedal, S.J. Wright, Numerical Optimization, Springer, USA, 1999. [22] R.S. Salgado, A.F. Zeitune, Power flow solutions through tensor methods, in: Proc. of the IET Generation, Transmission and Distribution, vol. 3, 2009, pp. 413–424. [23] D. Feng, R.B. Schnabel, Tensor methods for equality constrained optimization, SIAM Journal on Numerical Analysis 6 (1996) 653–673.