An Efficient and More Reliable Second Order Load Flow Solution Method with Cubic Interpolation Technique

An Efficient and More Reliable Second Order Load Flow Solution Method with Cubic Interpolation Technique

Asst. Prof. Dr. Hassan Kubba*

Electrical Engineering Department, Engineering College, Baghdad University

Al-Jadreyia Complex, Baghdad, Iraq,

Abstract:- This paper presents a fast, reliable, and new method for solving the load flow problem in electrical power systems. The proposed method is a second order load flow technique based on the "Taylor series expansion" of a multivariable function. This approach takes the first three terms of the Taylor series.

The method has advantages over Newton's method in terms of speed of solution (no. of iterations), and reliability of convergence. By inserting a minimization technique in this proposed method, the algorithm exhibits a control of the convergence. By means of this control, the method converges for cases when conventional Newton's method and some other popular methods diverge.

Also this paper presents a comparison between the proposed method and the Newton-Raphson method according to the major criteria, namely reliability of convergence, and speed of solution. Two test systems (five busbar typical test system and forty busbar typical test system) are used to examine the performance of each method.

Keywords:- Load flow problem, Taylor series expansion, Second order load flow model, Cubic interpolation techniques.

1. INTRODUCTION

   The load flow studies are the backbone in the planning of a power system. They are the means by which the future operation of the system is known ahead of time. A load flow study involves the determination of voltage, current, power, and power factor or reactive power at various points in an electrical network under existing or contemplated conditions of normal operation, subject to the regulating capability of generators, condensers, and tap changing under load transformers as well as specified net interchange between individual operating systems. This information is essential for the continuous evaluation of the current performance of a power system and for analyzing the effectiveness of alternative plans for system expansion to meet increased load demand. They are being increasingly used to solve very large systems, to solve multiple cases for such purposes as outage security assessment, and within more complicated calculation as optimization and stability.

   The studies for the load flow calculations started with the Ward and Hale method[1] in 1956, and currently the Newton

   Raphson method with the sparsity technique and sub-optimal ordering are being widely used[2]. The conventional (N-R) method has their relative advantages and disadvantages.

   A load flow problem consists of solving a set of nonlinear equations. The conventional (N-R) technique uses the first two terms of the Taylor series. This approach transforms the nonlinear load flow equations to a linear form before a solution is attempted. The proposed method formulate the load flow problem by using the first three terms of the Taylor series. In other words, second order terms, which are not insignificant, can be included in the algorithm and can be used during digital computation. The significant second order terms are found to be minor variations of the terms of the Jacobian matrix. It is shown in this paper that the coefficients of the second order terms are not required to be separately stored. In the proposed technique, the state vector is first calculated by an iteration of the conventional Newton-Raphson technique. Using the calculated state vector and elements of the Jacobian matrix, second order terms are estimated and subtracted from the residual vector. The modified residual vector obtained in this manner is then used to compute a new state vector. This procedure is repeated till the elements of the latest state vector are within permissible tolerance of those previously calculated. The magnitudes and phase angle of bus voltages are then updated. The total procedure is then repeated until a converged solution is obtained.

   The next problem sought to be solved in the proposed method is how to solve ill-conditioned systems or determine the existence of load flow solutions the proposed method is a combination of second order load flow algorithm and the cubic interpolation technique to determine an optimal multiplier for improving the load flow calculations. The algorithm exhibits a control of the convergence process and a non-divergent characteristic for any problem. Also the existence of the solution from the initial estimate can be easily judged.

2. TAYLOR SERIES EXPANSION

   A function can be evaluated by using the Taylor series expansion. The procedure is simple and well known for functions consisting of a single variable. Taylor series of multivariable functions can also be defined.

   For an (n) variable function, the series is expressed as follows:

   f x1 x1, x2 x2 ,……, xn x2 f x1, x2 ,……, xn

   m 1

   i

   i! x

   x1 x

   x2 …… x

   xn

   f x1 , x2 ,….., xn Rm

   (1)

   i1 1 2

   Where

   1

   n

   m1

   Rm m 1! x

   x1 x

   x2 ….. x xn

   f x1 a1x1, x2 a2x2 ,….., xn anxn

   1 2 n

   And 0 < ai < 1 ; i = 1, 2, ,, n

   Because the summation in equation (1) consists of (m) terms, a residue, Rm , is introduced to take care of the summation from (m+1)th term to infinity. This residue is not known exactly because definite values cannot be assigned to (ai's). Neglecting the third and higher terms in equation (1), equation (2) is obtained.

   f (x1 x1 , x2 x2 ……. xn xn ) f (x1 , x2 ,….., xn ) x x1 x x2 …… x xn f x1 , x2 ,….., xn

   1

   1 2 n

   2

   2 x

   x1 x

   x2 …. x

   xn

   f x1, x2 ,…., xn

   (2)

   1 2 n

   Equation (2) can be rearranged to give equation (3) which can be expanded as equation (4).

   f

   f x1 x1 , x2 x2 ,……, xn xn f x1 , x2 ,…., xn

   1

   2

   f x x x x … x x

   f x, x,…x

   2 x

   x1 x

   x2 …… x

   xn

   f x1, x2 ,…., xn

   1 2

   n

   (3)

   f n

   f x1, x2 ,…., xn x 1 n

   2 f x , x ,…., x

   2

   2

   i

   i

   1 2 n x

   2

   2

   i 1

   xi

   i

   i 1

   x2 i

   n

   n

   n 1

   2 f x , x ,….., x

   1 2 n x x

   i j

   i j

   (4)

   i 1 j i 1

   xixj

3. THE PROPOSED METHOD

   1. Second order load flow model

      In a power system, real and reactive power injected into a bus, say K, is equal to the net flow in the elements connected to this bus. Power flowing in an element connecting two buses, say K and m, can be defined in terms of the magnitudes and phase angles of the voltages at these buses and the parameters of the element. A load flow is, therefore, a problem of solving a set of nonlinear equations consisting of the magnitudes and phase angles of the system bus voltages as variables and the parameters of the system elements as constant coefficients. Let the power mismatch at a bus, K, of an (n bus system be defined as the difference between the scheduled power injection into this bus and the sum of the calculated power flows in all the elements connected to this bus. Real power mismatch can be defined in terms of the above mentioned variables by equation (5). This second order equation is similar to the Taylor series expansion excluding third and higher order terms in equation (4).

      Pk

      n

      n

      m1

      Pk

      m

      m

      m

      n

      n

      m1

      Pk V

      m

      m

      Vm

      1

      n

      n

      2 m1

      2

      P

      P

      k ( )2

      2 m m

      n

      n

      1

      k V

      n1

      n

      n

      2

      2

      P

      P

      k

      m r

      m r

      P

      P

      2

      2

      2 m1

      V 2 m m

      m1 r m1

      m r

      n n 2 P n n 2 P

      m r

      m r

      k V k V V

      (5)

      m1

      r 1

      m

      Vr

      m r

      m1

      r m1

      Vm

      Vr

      Where Pk is the real power injected into bus K,

      m is the phase angle of the voltage at bus m, Vm is the magnitude of the voltage at bus m,

      defines small changes in the variables.

      n

      n

      m

      m

      It is interesting to note that the terms of first two series in equation (5) are similar to the terms of the Jacobian matrix used in the Newton-Raphson method. The remaining five series constitute the second order terms. In a similar manner, reactive power mismatch can be defined as follows:-

      n

      n

      Q

      Qk

      Qk V

      1

      2Q 2

      k

      m1

      m

      m

      n

      n

      m1

      Vm

      2 m1

      k

      k

      2

      2

      m

      m

      m

      n

      n

      1

      2Q

      V

      n1

      n

      n

      2

      2Q

      m r

      k

      k

      m r

      k

      k

      2 m1

      V m 2

      m

      m1 r m1

      m r

      n

      n

      n

      n

      2

      Q

      Q

      k

      n1

      V

      2

      n

      n

      Q

      Q

      m r

      m r

      k V V

      (6)

      m1

      r 1

      m

      Vr

      m r

      m1

      r m1

      Vm

      Vr

      Expressing equation (5) for all system buses except the slack bus and equation (6) for all load buses, a set of equations is obtained in the following form.

      P

      J1

      J 2

      i

      Vh

      S1

      S2 S2

      i J

      V V

      (7)

      Q J

      J S

      S S h l

      3

      4 Vh 4 5

      6

      Vt

      s

      s

      The submatrices S1 through S6 include series of second order terms of all buses similar to those given in equations (5) and (6). Suitable values are assigned to subscripts i, j, h, l, s, and t. Many elements of the second order coefficient matrix, [S], are very small and can be neglected. Also, it is not essential to compute the elements of matrix [S] and store it in computer memory. These two aspects will now be discussed.

      Real and reactive power injection into a bus K of an (n) bus system can be mathematically expressed by equation (8) and (9) respectively. Second order coefficients can be derived from these equations and can be grouped into twenty categories. These coefficients are given in equations (10.1) through (10.20). a comparison with the elements of the Jacobian matrix indicates that the second order terms given in equations (10.1) through (10.10) are minor modifications of the Jacobian elements. These aspects are also indicated in these equations. The terms defined by equations (10.11) through (10.20) can be neglected as well be discussed later in this section.

      n

      Pk

      m1

      VkVmYkm Cos k

      m

      • km

        (8)

        n

        n

        Qk

        m1

        VkVmYkm Sink m km

        (9)

        Where Ykm is the admittance of the element connecting buses K and m and is given by

        YKm

        Km

        k

        k

        2 P

        V Y sin

        J1k ,k

        (10.1)

        n

        n

        k Vk

        m1

        m km

        k m km

        V

        V

        k

        V

        V

        2 P

        k

        k

        J1k ,m

        k

        k

        k

        Vm

        Vk Ykm sin k m km

        m

        (10.2)

        2 P

        k

        k

        J1k ,m

        V

        V

        m

        k

        k

        Vk

        VmYkm sin k m km

        k

        (10.3)

        V

        V

        2 P

        k

        k

        J1k ,m

        m

        k

        k

        Vm

        Vk Ykm sin k m km

        m

        (10.4)

        k

        k

        2Q

        k

        k

        2

        n

        n

        m1

        mk

        VkVm

        Ykm

        sin k

        m

      • km

        J1

        k, k

        (10.5)

        k

        k

        2Q

        m

        m

        2

        VkVmYkm sin k

        m

        km

        J1 k, m

        (10.6)

        k

        k

        2Q

        k

        k

        V 2

        2Ykk sin kk

        2

        V

        V

        2 J 4 (k, k)

        k

      • Qk

        (10.7)

        2Q

        k

        k

        k

        k

        k m

        VkVmYkm sin k m

        km J1 k, m

        (10.8)

        k

        k

        2Q

        V V

        Ykm sin k m km

        J1 k, m

        V V

        (10.9)

        k

        P

        P

        2

        k

        k

        k

        V 2

        2 P

        m

        2Ykk

        Cos kk

        2 J

        V

        V

        2

        2

        2

        2

        k

        k, k

        k m

        Pk

        (10.10)

        k 0

        m

        m

        V 2

        2Q

        k 0

        m

        m

        V 2

        (10.11)

        (10.12)

        k

        k

        2 P

        k

        k

        2

        n

        n

        m1

        mk

        VkVm

        Ykm

        cos k

        m

      • km

        (10.13)

        P

        P

        2

        k

        m

        m

        2

        VkVm

        Ykm

        cos k

        m

      • km

        (10.14)

        P

        P

        2

        k V V Y

        cos

        (10.15)

        cos

        cos

        (10.16)

        (10.16)

        k

        m k

        m km

        k m km

        k Y

        k Y

        2 P

        Vk Vm

        k

        k

        2Q

        km

        n

        n

        V Y

        k m

        cos

        km

        (10.17)

        k

        Vk

        m1 mk

        m km

        k m km

        k

        k

        n

        n

        2Q

        V Y cos

        (10.18)

        k

        Vm

        m1 mk

        m km

        k m km

        2Q

        k

        k

        m

        k

        k

        Vk

        VmYkm cos k

        m km

        (10.19)

        2Q

        k

        k

        m

        k

        k

        Vm

        Vk Ykm cos k

        m km

        (10.20)

        The elements defined by equations (10.11) and (10.12) are equal to zero. Equations (10.13) through (10.20) define the second order elements which include cos(k m km) as a multiplier. Since (km) is close to (90o) and (k m) is usually small, cos(k m km) is quite small and therefore, these coefficients can be neglected. This assumption and the simple relationship of the second order coefficients with the elements of the Jacobians matrix makes the application of second order load flow mode straight-forward and with minimal additional computing effort. The double summation is equation (5) and (6) have |Ykm| as multipliers. For any bus K, |Ykm| has non-zero values only for the buses (m) connected to it; for all other values of m, |Ykm| is zero. Therefore, the double summations for a bus can be reduced to single summation which include only a few terms depending on the number of busses connected to that particular bus. Deleting the significant terms from equations (5) and (6) and converting the double summations to single summation, equation (11) and (12) are obtained.

        n

        n

        P

        P

        2 P 2

        Pk

        m1

        k

        n

        n

        m

        m

        k

        k

        m1

        k V

        Vm

        k v

        k

        k

        n

        n

        m

        m

        V 2 K

        n

        n

        p

        p

        2

        k V

        2

        n

        n

        P

        P

        k

        V

        2

        P

        P

        k V

        (11)

        k V

        k V

        m1 m K

        k

        Vm

        K m

        m1

        mk

        m

        Vm

        m m m1

        m

        Vk m

        n

        n

        Q

        Qk

        Qk V

        1 2Q 2

        k

        k

        k

        m1

        m

        m

        n

        n

        m1

        Vm

        m 2 V 2 k

        n

        n

        1

        2Q

        2

        k

        k

        2

        2

        2Q

        n

        n

        k m

        k m

        2Q

        V V

        (12)

        2 m1

        m

        m

        n

        n

        m1

        mk

        k

        m

        k m

        m1

        mk

        k

        k

        Vk Vm

   2. Cubic Interpolation Technique [3, 4, 5]

   It is well-Known that the load flow calculation can be regarded as a nonlinear programming problem [6], which determines the direction and magnitude of the solution so that a certain function F(x) may be minimized. The F(x) is the squares of the active and reactive mismatch power. By employing this formulation, the valuable property can be obtained that the solution never diverges. The value of the function F(x) becomes eventually zero if there is a solution from the initial estimate, and stays at a positive value if no solution exists. In nonlinear programming approach (Fletcher-Powell method), (x) is modified by a correction factor () which can be considered as an acceleration factor. The computation of () is made such that minimize F(x). The function to be minimized is:

   FV ,

   kPQ,PV

   P2

   k

   k

   kPQ

   Q2

   (13)

   k

   k

   Where the first summation in (13), k includes load and generator buses, the second summation, k includes load buses only.

   The minimization of F(x) with respect to () in the direction of (x) is a one-dimensional problem. The object is to determine the correction factor () given (x) and the point (x). The problem can be stated as that of finding a value of () that will minimize:

   Fxo x

   And therefore the derivative of (Z) with respect to (x) is :

   ' 2 f yf yx

   y

   (14)

   (15)

   Where

   Y xo x

   A cubic interpolation technique is used to find (x) as follows:

   1. A step "a" is chosen as:

      a = min (1, -2 (Fx Fo ) / Z'x) (16)

      Where (Fo) is an estimate of (F) at the problem optimum, (Fx) is the value of function (F) at point (x), and (Z'x) is the derivative of (Z) with respect to (x) evaluated at point (x).

   2. A step of size (a) is taken to arrive at point (y), Y = x +x, and (Zy') is evaluated to determine whether a change of sign has occurred with respect to (Zx'). Such a change of sign, from negative (Zx') to positive (Zy') would indicate that the minimum is enclosed within these two points. If there is no change of sign, successive step of size (a) are taken until two adjacent points that enclose the minimum are found. Let these two adjacent points be called (w) and (y) which are located at

      distances ( w ) and ( y

      ) from the original point (x).

   3. The distance () from (x) to the minimum point is:

      y

      y

      • w

      y'r s

      (17)

      Where

      y y'Zw'2r

      1

      1

      r s 2 Zw' Zy'2

      (18)

      and

      S 3 Zw Zy Zw'Zy'

      y w

      (19)

   4. The point (xo + x) is accepted as the minimum point if the function F(xo +x) = Z(x) is smaller than both (Zw) and (Zy). If this is not true, the interpolation is repeated using the point (xo +x) and either (w) or (y), chosen so that the minimum is enclosed. This decision is based on the sign of Z'(). The cubic interpolation process is shown in fig (1).

      Fig (1) Cubic interpolation to minimize F(x) with respect to ( ) in the direction of ( x )

      The interpolation process can be simply implemented with the second order load flow approach.

      Equation (15) must be evaluated for several values of ()

      When 0

      f(y) = f(x) (20)

      And

      f y

      y

      f x

      x

      (21)

      Therefore

      2 f yf yx 2 f xf xx

      (22)

      y x

      From the equation f(x) = -J x Equation (21) becomes:

      Z ' 2 f x2 2Fx

      (23)

      The minus sign of equation (23) shows that the direction at minimization (x) always points in a direction which reduces F(x). The terms f(x), which represent the mismatch power, are calculated from the second order algorithm. In the

      proposed method, a simple alteration of the second order algorithm's program will allow a re-evaluation of the Jacobian at

      f y

      point (y) in order to obtain the term

      , also f(y) represents the mismatch power at point (y).

      y

      Then the interpolation process calculates the optimum () and the new point(x1 = x0 +x) is used for the start of the next iteration. The additional requirement of this method is the re-evaluation of the Jacobian matrix which is small, using sparsity techniques. There is no need for extra Gauss-Elimination and back substitutions, just additional evaluation of the Jacobian.

4. ITERATIVE ALGORITHM OF THE PROPOSED METHOD

The basic iterative algorithm for solving the load flow problem by the proposed method is as follows:

1. Set the initial voltage magnitudes of the busbars equal to that of the slack bus and usually (VK = V1 = 1 pu.) . All voltage angles are set equal to the slack bus ( K= 1 = 0), where bus (1) represents the slack bus.

2. Compute the real and reactive power mismatches and elements of the Jacobian matrix using the specified loads and

   p p

   p p

   generations, system parameters and estimated magnitudes and phase angles of bus voltages.

   p n

   p n

   P p jQ p V pe j k

   Y V pe j k V pe j k

   (24)

   k k k

   m1 mk

   km k

   m

   k,m=1,2..,n

   Where (*) complex conjugate

   And (n) no. of busses, and (p) represent the iteration index.

   Pp Psp Pp

   k k k

   Qp Qsp Qp

   (25)

   k k k

   (sp) indicates the specified value.

3. Neglect the second order terms, now the problem is the same as in Newton Raphson load flows. Using Gauss- elimination technique, evaluate the state vector [ V/V ]T.

4. For each P and Q of the residual vector, compute the second order terms (defined in equations 11 and 12) using

   's and V's evaluated in step 3.

   At this stage the residual vector is known and the sum of second order series for each residual term has been estimated. Transferring the second order terms to the left hand sde, equation (7) becomes

   P S S

   S i j J

   J i

   1 2

   3 V V 1

   2 Vk

   (26)

   Q S

   S S k

   L J

   J

   4

   5 6

   V 3

   4 Vk

   s t

   As a second order series is evaluated, it is subtracted from the corresponding residual term which procedure provides a modified residual vector defined by the left hand side of equation (26). It is important to realize that the elements of the matrix [S] are not computed and therefore, no additional storage is used for this matrix.

5. Using the modified residual vector and equation (26), recalculate the state vector [ V/V]T The triangularized Jacobian used in step (3) is reused at this stage.

6. (p+1) and (Vp+1) are modified by a correction factor () which can be considered as an acceleration factor. The computation of () is made such that it minimize F(x), which is the sum of the squares of the active and reactive mismatch power. An optimal correction is determined by finding an optimal () following the cubic interpolation technique.

7. Calculate the new values of the voltage magnitudes and angles for all busbars except the slack as follows:

   p1 p p1

   k k k

   And,

   Vk Vk V

   Vk Vk V

   p1 p p1 k

8. An iterative process is required until the expressions

   (27)

   k

   k

   P p1

   k

   k

   Q p1

   (28)

   are satisfied for all busses except the slack, where ( ) is a small power tolerance value. When the expression (28) are satisfied, the problem is solved. If not, the procedure is repeated with the next iteration, computing the elements of the Jacobian matrix and proceed to step (3).

   5. ANALYSIS BASED ON NUMERICAL EXAMPLES

   A 5-bus system taken from reference [7], and 40 bus system based on Iraqi National Grid were chosen to test the load flow solution methods. Tables (1) and (2) shows that the load flow problem was solved by Newton's method in 4 iterations to an accuracy of 10-3 for each individual mismatch power. In the proposed method 3 iterations were required with the same accuracy. The values of optimum () obtained for each iteration in the proposed method were close to one, and in the final iteration of a convergent load flow case, the value of () would be very close to one or exactly equal to 1.0. My experience has shown that optimum () is either close to 1.0 or very close to 0.0. The cubic interpolation formula will produce an appropriate value of () even in the case where the optimum is near, but outside the interpolating limits. So optimum () may sometimes be slightly greater than one. If the interpolation is performed between zero and one, the correction () value would be determined for all cases without any extra Jacobian calculation per iteration, thus saving computation time.

   The principle value of the proposed method lies in the control of the convergence process for both ill-conditioned and data error cases, whereas using the conventional Newton-Raphson method alone during the iterations of a load flow problem may result in poor solution or divergence.

   A study of the maximum power mismatches for both (N-R) method and the proposed method indicates what the maximum mismatches are smaller when second order algorithm (proposed method) are used. This indicates that the solution descends into the final state at a faster rate in the second order load flows than in the conventional (N-R) cases. Also, the contributions of the second order terms were noticed to be negligible after the first iteration.

   The performance of the proposed method with an ill-conditioned system was tested on the 5-busbar typical test system with insertion of series capacitors in four lines connecting the busbars for two different cases. The value of the capacitor is about 0.01 pu. on the same base of the other admittances. The solution by conventional Newton-Raphson method in the case, in which the diagonal elements of the [Y] matrix are less effective by the insertion of series capacitors, converged, while in the other case, the solution diverged. The proposed method converged for all cases.

   Also, a 40-busbar typical test system with insertion of series capacitors in ten lines connecting the busbars, including the slack, was chosen for this test. The solution by conventional (N-R) method diverged, while the proposed method converged with the value of optimum () as in table (3). The convergence characteristics of the 40-bus ill-conditioned system was shown in fig (2). From the figure, the following observation is made. In the case of optimum () is not applied, the solution continues to oscillate. However, if the optimum () is used the system converged easily.

   Table (1) Comparison of the proposed method and the conventional (N-R) method

   with the values of optimum ( ) for 5 bus bar typical test system, 4 load busbar and 1 slack .

   Method	No. of iterations	Optimum / iteration Proposed method
   Newton-Raphson Proposed method	4 3	1	2	3
   		1.02	1.02	1.003

   Notes: 1- Power tolerance = 10-3 for each individual mismatch power .

   1. This system is the same as that used in reference [7].

   2. Polar coordinates are used for both the proposed method and N-R method.

Table (2) Comparison of the proposed method and the conventional (N-R) method

with the values of optimum () for 40 busbar typical test system, 30 load busbar, 9 generator busbar, and 1

slack basbar[8].

Notes: 1- Power tolerance = 10-3 for each individual mismatch power .

1. This system is based on the Iraqi National Grid system.

2. Polar coordinates are used for both the proposed method and N-R method.

Table (3): Values of optimum () "Proposed method"

Fig (2) Convergence characteristics of 40 bus ill-conditioned system

102

101

10-0

10-1

10-2

10-3

Proposed method N – R

0 1 2 3 4

Fig (3). Convergence patterns of load flows of the 40 bus test system by (a) the Newton – Raphson approach

and, (b) proposed method.

1. CONCLUSION

   A very simple mathematical model and method for solving load flow problems has been developed. It would be concluded that the combined use of second order algorithm and minimization method by cubic interpolation technique has many attractive characteristics:

   1. The second order load flow model has some advantages over the conventional (N-R) approach. It is obvious that the proposed technique has better convergence characteristic, the magnitudes and phase angles of bus voltages descend into the final solution at a faster rate than that observed in (N-R) method. Therefore, in many cases, the second order approach requires lesser iterations than (N-R) technique. Moreover, theadditional computing effort is only slightly greater compared to a (NR) load flow. It has also been shown that the elements of the second order coefficient matrix need not be stored separately

   2. A more rapid convergence and a non-divergent characteristic for any problem. A control of the convergence process for both well-conditioned and ill-conditioned systems by the optimum

      correction factor ().

   3. The method can be applied efficiently to large systems.

   APPENDIX

   Nodal admittance matrix elements for 40 busbars "Typical Test System", ykm=G-jB Per unit quantity = 132 MVA, 132 KV all data in P.U.

   Bus to Bus	G(P.U), conductance	B(P.U), susceptance
   1-1	73.20302	-193.05308
   1-7	-43.763	109.66543
   1-36	-21.26757	53.28579
   1-38	-8.17245	30.10186
   2-2	15.35563	-32.0213
   2-3	-1.05799	4.58466
   2-20	-5.31241	11.95284
   3-3	9.62425	-26.04736
   3-15	-3.07535	7.70527
   3-28	-5.4909	13.75742
   4-4	32.25366	-95.89834
   4-14	-8.001	29.47035
   4-16	-4.80817	17.70843
   4-19	-12.8506	32.19863
   4-29	-6.59387	16.52092
   5-5	16.97918	-59.86251
   5-7	-8.30093	29.2611
   5-24	-8.67825	30.59639
   6-6	13.90746	-51.44682
   6-15	-6.53796	24.08149
   6-17	-7.4295	27.36532
   7-7	13.81057	-46.88945
   7-22	-1.76607	4.4249
   7-35	-3.74355	13.19844
   8-8	5.66517	-19.64528
   8-20	-2.48434	8.63067
   8-21	-3.18083	11.0146
   9-9	6.94882	-12.10935
   9-11	-5.14539	8.9666
   9-30	-1.80343	3.14275
   10-10	18.01201	-45.70547
   10-37	-9.62425	26.04736
   10-39	-5.31241	11.95284
   10-40	-3.07535	7.70527
   11-11	11.37208	-24.5679
   11-29	-6.22669	15.60129
   12-12	22.89817	-24.3483
   12-20	-3.20594	8.03246
   12-23	-19.69224	34.31584
   13-13	14.12767	-52.03695
   13-15	-4.97453	18.32287
   13-17	-9.15314	33.71408

   Bus to Bus	G(P.U.), conductance	B(P.U.), susceptance
   14-14	14.65299	-53.97187
   14-15	-6.65199	24.50151
   15-16	32.89171	-110.77638
   15-15	-8.17245	30.10186
   15-34	-3.47942	6.06339
   16-16	12.98062	-47.81029
   17-17	93.63839	-244.98614
   17-27	-5.05963	8.81716
   17-33	-6.96555	12.1384
   18-18	6.60826	-12.71954
   18-25	-0.90372	2.03338
   18-27	-4.23597	7.38181
   18-31	-1.46855	3.30425
   19-19	34.11817	-85.48443
   20-20	11.00269	-28.61598
   21-21	15.63407	-32.7162
   21-32	-3.56803	6.21779
   22-22	5.26614	-12.30004
   22-30	-1.59093	3.57961
   22-31	-1.90912	4.29553
   23-23	27.24222	53.23229
   23-26	-7.54998	18.91645
   24-24	9.75381	-33.01641
   24-25	-1.07556	2.42001
   25-25	11.97929	-4.4534
   26-26	9.43165	-32.55575
   26-28	-1.88166	4.6393
   27-27	9.29561	-16.19897
   28-28	7.37256	-18.39672
   29-29	12.82056	-32.12222
   30-30	3.39437	-6.72236
   31-31	3.37768	-7.59978
   32-32	3.56803	-6.21779
   33-33	6.96555	-12.1384
   34-34	3.47942	-6.06339
   35-35	3.74355	-13.19844
   36-36	21.26757	-53.28579
   37-37	9.62425	-26.04736
   38-38	21.02305	-62.30049
   38-40	-12.8506	32.19863
   39-39	5.31241	-11.95284
   40-40	15.92595	-39.9039

   Bus to Bus	G(P.U), conductance	B(P.U), susceptance
   1-1	73.20302	-193.05308
   1-7	-43.763	109.66543
   1-36	-21.26757	53.28579
   1-38	-8.17245	30.10186
   2-2	15.35563	-32.0213
   2-3	-1.05799	4.58466
   2-20	-5.31241	11.95284
   3-3		9.62425	-26.04736
   3-15	-3.07535	7.70527
   3-28	-5.4909	13.75742
   4-4	32.25366	-95.89834
   4-14	-8.001	29.47035
   4-16	-4.80817	17.70843
   4-19	-12.8506	32.19863
   4-29	-6.59387	16.52092
   5-5	16.97918	-59.86251
   5-7	-8.30093	29.2611
   5-24	-8.67825	30.59639
   6-6	13.90746	-51.44682
   6-15	-6.53796	24.08149
   6-17	-7.4295	27.36532
   7-7	13.81057	-46.88945
   7-22	-1.76607	4.4249
   7-35	-3.74355	13.19844
   8-8	5.66517	-19.64528
   8-20	-2.48434	8.63067
   8-21	-3.18083	11.0146
   9-9	6.94882	-12.10935
   9-11	-5.14539	8.9666
   9-30	-1.80343	3.14275
   10-10	18.01201	-45.70547
   10-37	-9.62425	26.04736
   10-39	-5.31241	11.95284
   10-40	-3.07535	7.70527
   11-11	11.37208	-24.5679
   11-29	-6.22669	15.60129
   12-12	22.89817	-24.3483
   12-20	-3.20594	8.03246
   12-23	-19.69224	34.31584
   13-13	14.12767	-52.03695
   13-15	-4.97453	18.32287
   13-17	-9.15314	33.71408

   Bus to Bus	G(P.U.), conductance	B(P.U.), susceptance
   14-14	14.65299	-53.97187
   14-15	-6.65199	24.50151
   15-16	32.89171	-110.77638
   15-15	-8.17245	30.10186
   15-34	-3.47942	6.06339
   16-16	12.98062	-47.81029
   17-17	93.63839	-244.98614
   17-27	-5.05963	8.81716
   17-33	-6.96555	12.1384
   18-18	6.60826	-12.71954
   18-25	-0.90372	2.03338
   18-27	-4.23597	7.38181
   18-31	-1.46855	3.30425
   19-19	34.11817	-85.48443
   20-20	11.00269	-28.61598
   21-21	15.63407	-32.7162
   21-32	-3.56803	6.21779
   22-22	5.26614	-12.30004
   22-30	-1.59093	3.57961
   22-31	-1.90912	4.29553
   23-23	27.24222	53.23229
   23-26	-7.54998	18.91645
   24-24	9.75381	-33.01641
   24-25	-1.07556	2.42001
   25-25	11.97929	-4.4534
   26-26	9.43165	-32.55575
   26-28	-1.88166	4.6393
   27-27	9.29561	-16.19897
   28-28	7.37256	-18.39672
   29-29	12.82056	-32.12222
   30-30	3.39437	-6.72236
   31-31	3.37768	-7.59978
   32-32	3.56803	-6.21779
   33-33	6.96555	-12.1384
   34-34	3.47942	-6.06339
   35-35	3.74355	-13.19844
   36-36	21.26757	-53.28579
   37-37	9.62425	-26.04736
   38-38	21.02305	-62.30049
   38-40	-12.8506	32.19863
   39-39	5.31241	-11.95284
   40-40	15.92595	-39.9039

2. REFERENCES

1. Stott B, "Review of load flow calculation method", proc. IEEE, Vol-62, pp. 916-929, July 1974.

2. Sasson A. M., "An optimal ordering algorithm for sparse matrix application", IEEE Trans. power App., Vol. PAS-97, No. 6, pp. 860-867, Nov. 1978.

3. Kubba H.A.,"A rapid and more reliable load flow solution method for ill-conditioned power system", Engineering and technology Journal Vol-17 , No.5, pp. 550-568. 1998.

4. Kubba H. A.,"A improved and more reliable decoupled load flow method", Engineering, Scientific Journal of Engineering college / Baghdad University, No.3 , Vol.7 , pp. 35-37. Sept. 2001.

5. Kubba H. A.,"An improved Newton method for radial distribution system load flow analysis, Engineering, Scientific Journal of Engineering College/ Baghdad University, No. 4, Vol. 12, pp. 1122-1135, September 2006.

6. A. K. Laha, K. E. Bollinger, R. Billinton and S. B. Dhar, "Modified form of Newton's method for faster load flow solutions", proc. IEE, Vol-121, pp. 849-853, Aug. 1974.

7. Stagg G. W., El-Abiad A. H., "Computer methods in power system analisi", Mc Graw Hill book Company. 1968

8. Kubba H. A.,"comparative study of different load flow solution methods", Al-Muhandis, Journal, Vol-107, pp. 25-46, 1991.

9. Aoki K. Nishikori A., "An algorithm for constrained load flow" IEEE transaction on pow. App. And systems. Vol. PAS-103, No.5, pp. 963-973, May 1984.

10. Kubba H. A., "An improved method of automatic adjustment of transformer and phase shifter Taps for constrained load flow", Al-Muhandis Journal, serial -148, No.4, pp.34-50. December 2001.

11. Charles A. Gross,"power system analysis", second edition, John Wiley publ. comp.1986.