Critical investigation of preconditioned GMRES via incomplete LU factorization applied to power flow simulation

Electrical Power and Energy Systems 33 (2011) 1695–1701

Contents lists available at SciVerse ScienceDirect

Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes

Critical investigation of preconditioned GMRES via incomplete LU factorization applied to power flow simulation José Eduardo O. Pessanha a,⇑, Carlos Portugal b, Ricardo Prada b, Alex R. Paz b a b

Power Quality Laboratory, Institute of Electrical Engineering, UFMA, Campus do Bacanga, São Luís, Ma 65080-040, Brazil PUC-Rio-DEE, Department of Electrical Engineering, Rua Marquês de São Vicente, 225, Gávea, RJ 22453-900, Brazil

a r t i c l e

i n f o

Article history: Received 25 September 2007 Received in revised form 16 March 2011 Accepted 12 August 2011 Available online 10 October 2011 Keywords: GMRES Preconditioning Incomplete factors Ordering Power flow

a b s t r a c t For solving the power flow sublinear problem efficiently by the GMRES preconditioned via incomplete LU factorization (ILU), this paper investigates causes associated to the preconditioner low quality and proposes a method to improve it and the GMRES convergence rate as well. The goal is provide a well-organized ILUGMRES for solving linear systems of difficult solution comprising ill-conditioned coefficient matrices, normally associated to heavy load power systems. The investigations reveal that a dropping rule for nonzero elements (fill-ins) based on a relative tolerance may introduce large errors during the preconditioner construction, lowering its quality and the GMRES performance. Based on that, it is proposed a fill-in dropping rule making use of two criteria; one based on the resulting error and the other based on a relative tolerance, applied to the preconditioner lower (L) and upper (U) triangular matrices, respectively. Ordering schemes are also considered. Numerical experiments taking into account different power system configurations operating under heavy load conditions corroborate the efficiency of such strategies. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction In order to supply and reach more costumers with a reliable and continuous service, power systems around the world are continuously increasing in terms of dimension, complexity and load level. As a result, computer simulations are facing large and complex mathematical systems, such as those found in the solution of large sparse linear systems of type (1), frequently the most timeconsuming part of the computation, demanding not only for high quality computer technologies but also for robust numerical methods. In general, these systems are solved through direct methods and despite of their robustness, they tend to have need of a predictable amount of resources in terms of computational time and storage [1,2], or even fail for ill-conditioned and/or very large coefficient matrices (A) [3]. Preconditioned Krylov subspace iterative methods have proven to be suitable alternatives for these cases; in general, their requirements for memory storage are usually smaller by orders of magnitude [4] and present relative simplicity of implementation in parallel algorithms, conversely for direct methods [5–7].

Ax¼b

ð1Þ

⇑ Corresponding author. Address: Avenida dos Portugueses s/n, Campus do Bacanga, São Luís, Ma 65080-040, Brazil. Tel.: +55 98 3301 9204. E-mail addresses: [email protected] (José Eduardo O. Pessanha), portugal@ele. puc-rio.br (C. Portugal), [email protected] (R. Prada), [email protected] (A.R. Paz). 0142-0615/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijepes.2011.08.010

One of the first applications of preconditioned iterative methods in power system problems was presented by Decker et al. [8]. The authors combined the LU factorization method with the Conjugate Gradient in transient stability analysis using parallel computation. In load flow analysis, Ref. [9] applied the Preconditioned Conjugate Gradient in dc and fast decoupled models. However, as iterative methods and preconditioned techniques improved, new contributions came up, manly related to the preconditioned Generalized Minimum Residual method – GMRES [10] and variants, such as GMRES(m) and GMRES-E [11]. The GMRES performance has been investigated in load flow applications [12,13], state estimation, power system security and transient stability [14]. Ref. [15] proposes an adaptive preconditioner constructed for Jacobian-free Newton-GMRES(m) method for solving coordination equations in distributed simulations of power systems. An extension of this work is presented in [16] where a new algorithm for distributed transient stability simulation of interconnected power systems based on a Jacobian-free Newton-GMRES(m) is proposed. Finally, Ref. [17] proposes a method comprising continued power flow and GMRES for voltage collapse analysis referred to as CPF-GMRES method. The present work is focused on the iterative solution of the power flow sublinear problem (Appendix A) of type (1) using a well-organized preconditioned GMRES method (Appendix B), with the preconditioning matrix constructed by means of incomplete LU factorization based on the power system Jacobian matrix calculated at the very first Newton–Raphson iteration. Normally, if the power system is heavy loaded, the condition number of the Jacobian matrix is very high resulting in a low quality ILU

1696

J.E.O. Pessanha et al. / Electrical Power and Energy Systems 33 (2011) 1695–1701

preconditioner, making the iterative process time consuming and the convergence may not be reached exceeding maximum iterations. In order to achieve a high quality ILU preconditioner, a critical investigation over its drawbacks was careful performed. For the sublinear power flow problem, it was found that a well know ILU dropping rule based on a relative tolerance might introduce large errors, resulting in a low quality preconditioner. Based on that, it is proposed a fill-in dropping rule during the ILU preconditioner construction making use of two criteria; one based on the resulting error and the other on a relative tolerance. The first criterion is applied only over the lower triangular matrix (L) and the second only over the upper (U). Another proposal, which in fact improves the preconditioner quality, is ordering the Jacobian matrix before the preconditioner construction. Therefore, an appropriate ordering scheme and the composite dropping rule form the proposed method to improve the ILU quality and the GMRES performance. Preconditioning is considered by many as the most important element in the development of efficient iterative algorithms for challenging problems in scientific computation, and the importance of preconditioning is destined to increase even further [18]. This statement gives a special motivation to the development of the present work. 2. Incomplete LU preconditioner The accuracy of the ILU preconditioners with no fill-in may be insufficient to offer an adequate rate of convergence, such as the ILU(0) preconditioner. However, allowing some fills they become more accurate, efficient and reliable. To allow a controlled increase of fill-ins, one can make use of the ILU(k) preconditioner, where k specifies how many levels of fill-ins are permitted, but these methods do not offer robustness in terms of reliability and fast convergence. A few alternative methods are available based on dropping elements in the Gaussian elimination process according to their magnitude rather than their locations. In a well known incomplete factorization preconditioner with threshold [19], known as ILUT(q, s), fill-ins are dropped according to a relative tolerance (strategy proper for badly scaled matrices [18]) obtained by multiplying an absolute tolerance (s) by the original 2-norm of the ith row keeping only the q largest elements in L and U parts of the row (Appendix C). Despite of being considered by many [4,18,20–22], this dropping rule may introduce errors during the preconditioner construction, decreasing its quality, as shown bellow. 3. ILUT drawback

Table 1 Dropped fill-ins (conventional rule). lij

|lij|

Tolerances

l10,2 l6,5 l8,6 l9,5 l8,5 l9,6 l7,2 l10,4 l9,8 l8,7 l10,7 l7,3 l4,3 l2,1

1.31E18 7.62E16 7.66E16 7.66E16 6.62E15 6.62E15 6.68E15 0.0095 0.0156 0.0241 0.0247 0.0328 0.0332 0.0612

Error

s1

s2

v

X X X X X X X – – – – – – –

X X X X X X X – – – – X – X

0.0010 0.0018 0.0018 0.0018 0.0076 0.0076 0.0077 0.0328 3.2377 4.8525 4.9734 8.0816 8.1728 4.0129

the Floating Point Operations per Second (FLOPS), shown in Table 2. This table also includes information regarding to unpreconditioned GMRES (UN) and to a complete factor preconditioner (LU). The performance of the UN-GMRES is very low presenting a high number of operations (16,164). With a LU preconditioner, the number of operations was reduced to 413, with 199 associated to the preconditioner construction and 214 to GMRES. Making use of ILUT(s1), the number of operations associated to the preconditioner construction was reduced even further (155), but did not change to GMRES (214), totalizing 369 operations. For ILUT(s2), the number of operations associated to its construction was also reduced (143), but increased to GMRES (428), resulting in 571 operations. This happened because the next elements to be dropped should be l10,4, l9,8, . . . , l2,1, but instead elements l7,3 and l2,1 were dropped resulting in large errors, lowering the preconditioner quality, justifying the increasing in GMRES operations. 4. The composite dropping rule In order to avoid the previous problems, it is proposed a fill-in dropping rule for the elements of L based on the resulting error. Assuming L as the complete lower triangular factor and L as the incomplete after dropping a single element (lik), and rewriting (2) as in (3), one can verify that ðL  LÞ has only the fill-in lik (dropped in L), as shown in (4) [18,23].

Errorðlik Þ ¼ kðL  LÞ  Uk1 ; i > k 0 6 .. . . 6. . 6 6 ðL  LÞ ¼ 6 0    lik 6 .. 6 .. 4. .

Errorðlik Þ ¼ kL  U  L  Uk1 ; i > k

Table 2 FLOPS.

ð2Þ

Elements above the solid black line are smaller in magnitude than the absolute tolerances, and those bellow are larger. For s1 were dropped the seven smaller elements and for s2, two more elements were also dropped out of sequence (l7,3 and l2,1). The solution was obtained but the consequences of these droppings above the absolute tolerance can be better visualized computing

3

2

To illustrate the ILUT dropping rule weakness when applied to the power flow sublinear problem, two numerical experiments are performed using a small-port (5-bus and two generators) hypothetical power system model for different absolute tolerances, s1 = 2.0E04 and s2 = 2.5E05, and since it results in a small-port linear system the parameter q is not taking into account. Table 1 shows the absolute values of each element of L (columns 1 and 2), the dropped ones (X – marked, column 3) and the resulting error given by (2) associated to each dropping (column 4), where L represents the incomplete lower triangular matrix [18,23].

0 

ð3Þ

..

0

.

7 7 7 7 7; i > k 7 7 5

ð4Þ

 0

Replacing (4) in (3):

Preconditioner GMRES Total

UN

LU

– 16,164 16,164

199 214 413

ILUT

s1

s2

155 214 369

143 428 571

J.E.O. Pessanha et al. / Electrical Power and Energy Systems 33 (2011) 1695–1701

2

0  0 .. 6 .. 6. . 6 6 ðL  LÞ  U ¼ 6 0    lik  ukk 6 .. 6 .. 4. . 0 

0



0 .. .

   lik  uki .. . ...

0



0 .. .

(

3

7 7 7    lik  ukn 7 7; i > k 7 .. 7 . 5 

0 ð5Þ

Resulting for the error:

Errorðlik Þ ¼ jlik j  maxfukm g; i > k; m ¼ ½k; . . . ; n

ð6Þ

Note that the 1-norm is used in (3) in order to reduce the number of numerical operations associated to the dropping rule used during the preconditioner construction. This norm takes advantage of the matrix structure in (5) avoiding additional operations, needing only comparisons among the elements ukk, uk,k+1, . . . , uki, . . . , ukn. The ILU factorization used here is based on Gaussian elimination (Appendix C), known as IKJ [20] and it is suitable for sparse matrices, like those normally found in power flow and time domain analysis. In order to compute the error from dropping lik, all ukm elements already calculated in previous steps are needed. On the other hand, for the error associated to the dropping of uik, one must know all lkm, but this is not possible due to the IKJ scheme. Therefore, the dropping rule based on error is applied only over fill-ins of the lower triangular matrix L and a variant of the conventional rule over fill-ins of the upper triangular matrix U, resulting on the following composite dropping rule:  A lik element is dropped if the resulting error (6) is smaller than a v tolerance used for all rows of L, as given by (7). Otherwise the element is calculated using IKJ. An element uik is dropped if its absolute value is less than a relative tolerance obtained by multiplying an absolute tolerance (s) by the arithmetic average of the ith row elements of the original Jacobian matrix, as presented by (8).

lik

¼

k¼ð1;...;i1Þ

uik ¼



0; a0ik ukk

Errorðlik Þ < v ; Errorðlik Þ P v

xk < s xk ; xk P s k ¼ ði þ 1; . . . ; nÞ

0;

1697

ð7Þ

ð8Þ

The IKJ factorization allows the application of the dropping rule based on the error only over L fill-ins. However, if L presents more fill-ins, the use of the composite dropping rule, which the resulting preconditioner is identified in this paper as ILU(v, s), may be still attractive. Using three different electric power systems it is possible to set up a rate (%) between L and U fill-ins for values bellow an absolute tolerance (s1 and s2 previously chosen) that should be dropped. The results are illustrated in Fig. 1 and corroborate the fact that the greater part of fill-ins to be dropped is found in L.

5. Ordering and preconditioning the Jacobian matrix Ordering a sparse matrix is not a new practice for improving the performance of a direct method for solving linear systems [1,2]. In this paper, ordering the power system Jacobian matrix before the preconditioner construction results in substantial benefits, once this strategy reduces fill-ins, improve the stability of the incomplete factorization and GMRES convergence rate [24]. Investigations performed by the authors compared different ordering schemes, such as nested dissection [25], Sloan’s [26], Gibbs et al. [27], reverse Cuthill-McKee [28] and minimum degree [29]. The RCM was already considered by Portugal [3] and Flueck and Chiang [30] in the construction of an ILU preconditioner for the power flow problem and minimum degree has been considered to improve direct methods performance [31]. The Jacobian matrix J is ordered through symmetrical permutations using reverse Cuthill McKee or minimum degree, where P is the permutation matrix, resulting in the modified linear system

(a) IEEE 118 - bus test system (τ 1 and τ 2)

(b) Brazilian North-Northeast Power System 257 - bus test system ( τ1 and τ2 )

(c) Brazilian Power System - 3315 - bus test system ( τ 1 and τ 2 ) Fig. 1. Rate of dropped fill-ins in L and U triangular matrices.

1698

J.E.O. Pessanha et al. / Electrical Power and Energy Systems 33 (2011) 1695–1701

given by (9) to be solved at the kth Newton iteration, which solution is the vector DZ.

½PT  J  P  ½DZ ¼ PT 



½DP=Vn

k

½DQ =Vn

ð9Þ

2n

After ordering the original linear system, it is left-preconditioned through the inverse of the transformation matrix M (preconditioner), as given by (10). As closer the transformation matrix is to the inverse of the ordered Jacobian matrix (M1 ffi (PT  J  P)-1 ffi [I]), as easier the solution is obtained (11).

M

1

T

 ½P  J  P  ½DZ ¼ M

½DZ ffi M

1

P

T



½DP=Vn ½DQ =Vn

1

T



P 

½DP=V

k

½DQ =Vn

Table 3 Preconditioner parameters – Experiment I. Parameter

ILU(k) MD

ILU(k) RCM

ILUT(q, s) MD

ILUT(q, s) RCM

s q

– – 20

– – 10

4.1E12 95 –

1.0E07 18 –

k

Table 4 Preconditioners and ordering schemes performance. Preconditioner

ð10Þ

2n

k

ILU(k) ILUT(q, s)

CPU time (s)

MFLOPS

None

RCM

MD

RCM

MD

27.980 25.437

0.390 0.375

0.203 0.187

11.621 10.129

2.902 2.672

ð11Þ

2n

Because the calculation of the complete triangular factors (LU) demands a substantial computational effort, incomplete factors ðLUÞ are considered to construct the preconditioner (or transformation) matrix, as given by (12). Once the system is ordered and preconditioned, the coefficient matrix A and the right-side vector b (13) are ready for the GMRES method, which residual r for each linear iteration i is given by (14).

M ¼ L  U ffi PT  J  P A ¼ M1  ½PT  J  P b ¼ M1  PT 

ð12Þ 

½DP=Vn ½DQ =Vn

k

6.1. Experiment I: ordering strategies and dropping rules

ð13Þ

2n

  ½DP=Vn k  M1  ½PT  J  P  ½DZi  ri ¼ M1  PT  ½DQ =Vn 2n |fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl} |ffl{zffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} xi A

ð14Þ

b

Since the transformation matrix is not explicitly known, the residuals are calculated considering incomplete triangular factors, as presented in (15), applying forward-substitution to calculate L1 and backward substitution to calculate U1 .

A  xi ¼ U1  ½L1  ½ðPT  J  PÞ  xi  b   k !! ½DP=Vn ¼ U1  L1  PT  ½DQ =Vn 2n

by Eqs. (9)–(15) for solving the power flow sublinear problem via preconditioned GMRES. All experiments were run on a 1.73 GHz Intel Centrino, 1 GB ram computer and the code was written using Compaq Visual FORTRAN Professional Edition 6.1.0. The tolerances for the nonlinear (Newton–Raphson) and linear (GMRES) systems have been set, respectively, as 103 (MW and MVAr power mismatches) and 105. The preconditioners are kept fixed in all experiments.

ð15Þ

6. Numerical experiments In this section several numerical experiments are performed in order to verify the efficiency of the proposed strategies formulated

The RCM and minimum degree (MD) ordering schemes are evaluated taking into account ILU(k) and ILUT(q, s) preconditioners for solving the sublinear power flow problem using the 3515-bus and 301 generators (fixed MVAr control) Brazilian power system, operating under heavy and real load conditions according to the Brazilian National Power System Operator for the year of 2006 (http://www.ons.org.br). The preconditioners parameters given in Table 3 were established after several simulations, with the CPU time used as criterion (the lower the better), but a simple loop made the search inexpensive. 6.1.1. Results Fig. 2 illustrates the GMRES convergence rate and Table 4 presents the associated CPU time and MFLOPS spent for solving the sublinear problem. The highest GMRES convergence rate is achieved when considering minimum degree (MD) ordering scheme and the ILUT preconditioner. Therefore, for power flow problems, avoiding fill-ins (minimum degree) rather than confining nonzero elements in a narrow bandwidth (RCM) results in a better quality ILU preconditioner.

Fig. 2. GMRES convergence rate using different ordering schemes. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

1699

J.E.O. Pessanha et al. / Electrical Power and Energy Systems 33 (2011) 1695–1701 Table 5 Performance of both dropping rules. Preconditioner

CPU time (s)

Fill-ins

MFLOPS

ILUT(q, s), q = 95 s = 4.1E12 ILU(v, s), v = 1.0E02 s = 1.0E04

0.187 0.130

15,314 14,720

2.6723 2.4591

Table 6 IEEE-118 bus – CPU time (seconds). MD + ILU(v, s) + GMRES, v = 1.0E02 s = 1.0E04 ILU(v, s) + GMRES, v = 1.0E02 s = 1.0E04 MD + ILUT(q, s) + GMRES, q = 20 s = 10E07 RCM + ILUT(q, s) + GMRES, q = 15 s = 10E06 SDPF UN-GMRES MD + LDU

Fig. 3. Fill-ins percent rate. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

The convergence rate is improved because the GMRES residual (15) is modified after ordering and preconditioning the original Jacobian matrix and since this rate changes as the residual changes, and once the residual is used to construct the Krylov subspace, one can conclude that the Krylov subspace is modified after ordering and preconditioning the original matrix, shorting the distance to the solution. Fig. 3 illustrates the fill-ins percentage droppings where one can see that the smaller occurs for minimum degree ordering scheme. There are two reasons for that. First, a large amount of fill-ins have already been dropped (see Fig. 4). Second, the ILUT(q, s) dropping rule seems to be inefficient in this case. Using the ILU(v, s) instead, more fill-ins are dropped during the preconditioner construction, reducing the number of operations and the CPU time even further improving the overall efficiency (see Table 5). 6.2. Experiment II: voltage instability The present set of experiments considers the IEEE-118 bus power (http://www.ee.washington.edu/research/pstca/pf118/pg_tca118-

0.008 s 0.031 s 0.062 s 0.093 s 0.563 s Fail Fail

bus.htm) system with 54 generators (fixed MVAr control) operating extremely loaded and despite of using a hypothetical model, the scenario is associated to voltage instability, which may lead to voltage collapse as experienced by many power systems around the world [32,33]. This load level was obtained through a continued power flow computer program [34] corresponding to the last converged load flow. Besides the different iterative combinations, two other sources are used for comparative purposes; one based on LDU factorization with a minimum degree ordering scheme (MD + LDU) and the other using a new solution concept, known as Synthetic Dynamics Power Flow – SDPF [35]. The authors do not intend to disqualify these sources, but since both are used by many Brazilian power system utilities and universities as well, they seem to be reliable tools to corroborate the proposed method. 6.2.1. Results Table 6 presents the resulting CPU time in seconds where one can see that the proposed method reached the solution presenting the lowest CPU time, followed by ILU(v, s) preconditioner without any ordering scheme. The following iterative combinations differ based on the ordering scheme and once again the results corroborate the minimum degree performance over RCM. The SDPF reaches the solution spending 0.563 seconds, while UN-GMRES and MD + LDU fail. Fig. 5a shows the eigenvalues distribution of the original Jacobian matrix and the number in parenthesis is the corresponding condition number, both indicating that the matrix is originally very ill-conditioned justifying the two failures and the high CPU time for SDPF. After ordering and preconditioning via MD + ILU(v, s), the eigenvalues of the modified Jacobian matrix become tightly clustered (see Fig. 5b) improving the GMRES convergence rate. The problem is solved by all iterative combinations with MD + ILU(v, s) presenting an outstanding performance. The accuracy of ILU(v, s) and ILUT(q, s) preconditioners is given in Table 7 and it is obtained using the Frobenius norm (N1) [24] which measures the distance between the Jacobian matrix J (or coefficient matrix A) and the product of the incomplete factors. It is also shown the total number of linear/nonlinear iterations and MFLOPS. The GMRES reached convergence presenting the same number of Table 7 Preconditioner accuracy.

Fig. 4. Preconditioner fill-ins. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Preconditioner

N 1 ¼ kJ  LUkF

Linear iterations (GMRES)

Nonlinear iterations (Newton– Raphson)

MFLOPS

ILU(v, s) ILUT(q, s)

30 600

14 14

14 14

0.1978 0.2360

1700

J.E.O. Pessanha et al. / Electrical Power and Energy Systems 33 (2011) 1695–1701

(a) Original (2,792.687)

(b) Modified (1.1046)

Fig. 5. Jacobian matrix eigenvalues distribution.

linear iterations (as well as for nonlinear iterations) for both preconditioners. However, the Frobenius norm indicates that the quality of ILU(v, s) is much better than the quality of ILUT(q, s) justifying the GMRES fast convergence with a smaller MFLOPS. 7. Conclusions This paper proposes a method for constructing an efficient ILU preconditioner for solving the power flow sublinear system by the GMRES method using serial computation, focusing on very ill-conditioned Jacobian matrices associated to extreme load conditions, closer to voltage collapse. It was shown that a fill-in dropping rule based on relative tolerance may introduce errors during the preconditioner construction lowering its quality and the GMRES performance as well. An efficient dropping rule based on the error is proposed to overcome the drawbacks and to reduce the amount of fill-ins before the preconditioner construction the method makes use of a minimum degree scheme for ordering the original Jacobian matrix (preconditioner image). Numerical experiments showed that as more loaded the power system is, more difficult the solution of the resulting linear system will be because the original Jacobian matrix presents a very high condition number. The results corroborate the best performance of the proposed method over selected ones, converging rapidly within a low number of linear iterations, smaller CPU time and MFLOPS as well, avoiding the necessity for restarting the GMRES. Despite of the method be focused for planning and operation of power systems, it may be applicable to other scientific problems.

The above nonlinear equations can be solved by Newton–Raphson method (A3). The N  N Jacobian matrix (J) is normally highly sparse and formed by four submatrices H, N, M and L, as shown in (A4), including the active and reactive power mismatches (DP and DQ), as well bus voltage magnitudes and phase angles. These submatrices are structurally identical and symmetric; therefore the Jacobian matrix is also structurally symmetric, despite of numerically asymmetric. The power flow sublinear problem is then given by (A5).

xði þ 1Þ ¼ xðiÞ þ J1 ðiÞfy  f ½xðiÞg " J¼  H  M

@ DP @h @ DQ @h

@ DP @V @ DQ @V

#

N



L

     DhðiÞ DPðiÞ ¼  L DVðiÞ DQ ðiÞ

N

ðA4Þ

ðA5Þ

Appendix B. GMRES The Generalized Minimum Residual method, or simply GMRES, was introduced by Decker et al. [8] for solving non-symmetric linear equations. This method searches the solution over a minimum residual norm, identifying the solution xk such that its Euclidian norm is minimal over the Krylov subspace, which bases are orthonomalized by Arnoldi process. An approximated solution xk is estimated using (B1), and V (B2) forms the base vectors of the Krylov subspace after orthonomalized.

xk ¼ x0 þ V  yk h ½Vnm ¼ v f1 v f2

Acknowledgments

 H ¼ M

ðA3Þ

...

v

f k

i

ðB1Þ ðB2Þ

The authors are grateful to the following for providing financial support to this work: Brazilian Federal Research Agency, CNPq, under the Grant PADCT 620097/2004-3, and Brazilian North Utility – ELETRONORTE, under the Grant P&D 45000044031.

Eq. (A2) is used by Arnoldi’s algorithm to orthonomalize the  vector v ikþ1 (B3) referred to other k vectors v f1 ; v f2 ; . . . ; v fk , previously orthonomalized.

Appendix A. The power flow sublinear problem

vfkþ1 ¼

Under balanced three-phase steady-state conditions, a basic power flow simulation computes voltage magnitude (V) and angle (h) at each bus k in a power system, real (A1) and reactive (A2) power flows through admittance Ykn interconnecting the bus-k to bus-n [36].

Pk ¼ V k

N X

Y kn V n cosðdk  dn  hkn Þ

ðA1Þ

n¼1

Qk ¼ Vk

N X n¼1

Y kn V n sinðdk  dn  hkn Þ k ¼ 1; 2; . . . ; N

ðA2Þ

P

vikþ1  kr¼1 ½v fr  vikþ1   v fr P kv ikþ1  kr¼1 ½v fr  v ikþ1   v fr k

ðB3Þ

Once V and its base vectors are estimated, one calculates yk based on the minimum residual norm, minimizing the residual Euclidian norm to obtain a good approximation to the solution x  (B4). Using (B1) and (B3), the residual Euclidian norm may be represented by (B5) [13]. As the norm must be minimized at each GMRES k iteration, then vector yk is estimated from solution of (B6) (minimum mean square).

krk k2 ¼ kb  A  xk k2

ðB4Þ

krk k ¼ kv kþ1  ðkr0 k2  ^e1  Hk  y2 Þk2

ðB5Þ

J.E.O. Pessanha et al. / Electrical Power and Energy Systems 33 (2011) 1695–1701 Table A1 Experiment I - Ordering Strategies. Parameter

ILU(k) MD

ILU(k) RCM

ILUT(q, s) MD

ILUT(q, s) RCM

sA q

– – 20

– – 10

4.1E12 95 –

1.0E07 18 –

k

Hk  yk ¼ kr0 k2  ^e1

ðB6Þ

From (B4)–(B6), e1 is a canonical vector of dimension-k, Hk is a Hessenberg matrix of dimension (k + 1)  k, which elements are calculated during Arnoldi process (B7). This equation set results in the GMRES algorithm (no preconditioned) illustrated bellow.

hi;k ¼ v fi  ½A  v fk      k X  f f f f hkþ1;k ¼ A  v k  ½v i  ½A  v k   v i    i¼1

ðB7Þ

GMRES algorithm Step Step Step Step Step Step Step Step

1 2 3 4 5 6 7 8

Step Step Step Step

9 10 11 12

Compute: r0 = b  Ax0, b = ||r0||2 and v1 = r0/b For j = 1, 2, . . . , k, Do: Compute: wj = Avj For i = 1, . . . , j, Do: hij = (wj, vi) wj = wj  hijvi End-Do Compute: hj+1,j = ||wj||2. If hj+1,j = 0 adjust m = j and go to step 10 Compute: vj+1 = wj/hj+1,j End-Do Set up Hessenberg matrix b k  yk and Find yk to minimize: kb  e1  H 2

xk = x0 + Vk  yk

Appendix C. ILUT preconditioner algorithm

Step Step Step Step Step Step Step Step Step Step Step Step Step Step

1 2 3 4 5 6 7 8 9 10 11 12 13 14

For i = 1, . . . , n Do x = ai For k = 1, . . . , i  1 and when xk – 0, Do xk = xk/akk Apply a dropping rule to xk If xk – 0 then x = x  xk  uk End-If End-Do Apply a dropping rule to row x lij = xj for j = 1 ,. . . , i  1 uij = xj for j = 1, . . . , n x=0 End-Do

Appendix D. Numerical experiment data Table A1. Experiment I – ordering strategies. References [1] Duff IS, Erisman AM, Reid JK. Direct methods for sparse matrices. Oxford: Claredon; 1989. [2] George A, Liu JW. Computer solution of large sparse positive definite systems. Englewood Cliffs, NJ: Prentice-Hall; 1981.

1701

[3] Portugal CE. An iterative solver for linear systems: Application in the power flow problem. D.Sc. thesis. Catholic University of Rio de Janeiro, Department of Electrical Engineering; 2010 [in Portuguese]. [4] Van der Vorst H. Iterative Krylov Methods for Large Linear Systems. Cambridge monographs on applied and computational mathematic. 1st ed.; 2003. [5] Alvarado FL. Parallel solution of transient problems by trapezoidal integration. IEEE Trans Power Apparat Syst 1979;PAS-98(3):1080–90. [6] La Scala M, Bose A, Tylavsky DJ, Chai JS. A highly parallel method for transient stability analysis. In: Power industry computer application conference, PICA ‘89, conference papers; 1–5 May 1989. p. 380–6. [7] Chai JS, Bose A. Bottlenecks in parallel algorithms for power system stability analysis. IEEE Trans Power Syst 1993;8(1):9–15. [8] Decker IC, Falcão DM, Kaszkurewics E. Parallel Implementation of a power system dynamic simulation methodology using the conjugate gradient method. IEEE Trans Power Syst 1992;7:458–65. [9] Galiana FD, Javidi H, McFee S. On the application of a preconditioned conjugate gradient algorithm to power network analysis. In: Proceedings of PICA conference; 1993. p. 404–10. [10] Saad Y, Schultz MH. GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J Sci Stat Comput 1986;7:856–69. [11] Chianotis D, Pai MA. A new preconditioning technique for the GMRES algorithm in power flow and P–V curve calculation. Electr Power Energy Syst 2002(25):239–45. [12] Semlyen A. Fundamental concepts of Krylov subspace power flow methodology. IEEE Trans Power Syst 1996;11:1528–37. [13] Flueck AJ, Chinag Hsiao-Dong. Solving the nonlinear power flow equations with a Newton process and GMRES. In: IEEE international symposium on circuits and systems, vol. 1; 1996. p. 657–60. [14] Pai MA, Dag H. Iterative solver techniques in large scale power system computation. In: Proceedings of the 36th IEEE conference on decision and control, vol. 4; 1997. p. 3861–6. [15] Chen Y, Shen C. A Jacobian-Free Newton-GMRES(m) method with adaptive preconditioner and its application for power flow calculations. IEEE Trans Power Syst 2006;21(3):1096–103. [16] Chen Y, Shen C, Wang Jian. Distributed transient stability simulation of power systems based on a Jacobian-Free Newton-GMRES method. IEEE Trans Power Syst 2009;24(1):146–56. [17] Jasni J, Hizam H, Kadir MZA, Mariun N, Noor SBM. determination of proximity to static voltage collapse using CPF-GMRES method. In: 2nd IEEE international conference on power and energy (PECon 08), Johor Baharu, Malaysia; December 1–3, 2008. p. 520–5. [18] Benzi M. Preconditioning techniques for large linear systems: a survey. J Comput Phys 2002;182:418–77. [19] Saad Y. ILUT: a dual threshold incomplete LU factorization. Numer Linear Algebra Appl 1994;1(4):387–402. [20] Saad Y. Iterative methods for sparse linear system. 2nd ed. Philadelphia: SIAM, Society for Industrial and Applied Mathematics; 2003. [21] Chen K. Matrix preconditioning techniques and applications. Cambridge monographs on applied and computational mathematics; 2005. [22] Zhang J. A multilevel dual reordering strategy for robust incomplete LU factorization of indefinite matrices. SIAM J Matrix Anal Appl 2000;22(3):925–47. [23] Benzi M, Szyld DB, Van Duin A. Orderings for incomplete factorization preconditioning of nonsymmetric problems. SIAM J Sci Comput 1999;20(5):1652–70. [24] Benzi M, Tuma M. Orderings for factorized approximate inverse preconditioners. SIAM J Sci Comput 2000;21(5):1851–68. [25] George A. Nested dissection of a regular finite element mesh. SIAM J Numer Anal 1973;10:345–63. [26] Sloan SW. An algorithm for profile and wavefront reduction of sparse matrices. Int J Numer Methods Eng 1986;23:1315–24. [27] Gibbs NE, Poole Jr WG, Stockmeyer PK. An algorithm for reducing the bandwidth and profile of a sparse matrix. SIAM J Numer Anal 1976;13(2):236–50. [28] Cuthill E, McKee J. Reducing the bandwidth of sparse symmetric matrices. In: Proceedings of the 24th National Conference ACM, 1969. p. 157-72. [29] Tinney W. Direct solutions of sparse network equations by optimally ordered triangular factorization. Proc IEEE 1967;55(11):1801–9. [30] Flueck AJ, Chiang HD. Solving the nonlinear power flow equations with an inexact Newton method using GMRES. IEEE Trans Power Syst 1998;13(2):267–73. [31] Peschon J, Bree Jr DW, Hajdu LP. Optimal power-flow solutions for power system planning optimal power-flow solutions for power system planning. Proc IEEE 1972;60(1):64–70. [32] NERC. Survey of the voltage collapse phenomenon. North American Reliability Council; 1991. [33] Gustafsson M, Krantz N. Voltage collapse in power systems – analysis of component related phenomena using a power system model. Technical report no. 215L. Chalmers University of Technology, Department of Electrical Power Engineering, Sweden; 1995. [34] CEPEL. Anarede User´s Manual – version V07-08/99; 1999. [35] Jardim J, Stott B. Synthetic dynamics power flow. In: IEEE power engineering society general meeting, vol. 1; 2005. p. 479–84. [36] Glover JD, Sarma M. Power system analysis and design. Boston: PWS Publishers; 1987.