Chapter 5: Direct Methods

CHAPTER 5: DIRECT METHODS IN Ybus FORMATION.

5.1 Positive Sequence Ybus Formation
5.1 The Method of Singular Transformation (MST)
5.2 Example Calculation
5.3 Computer Technique for MST
5.4 Example Calculation, Direct Method, K.Nagappan Style
5.5 Another Computer Technique for MST, Direct Method
5.6 Matrix Multiplication, MST Program (SINGULAR.C)
5.7 Computer Program for Direct Method (DIRECT.C)
5.8 Positive Sequence Ybus Formation
5.9 Example Calculation
5.10 Treatment of Parallel lines with Coupling
5.11 Processing Coupling on a per Group Basis

The Ybus formation for short circuit in the positive sequence is similar as that used in loadflow. Line charging, off-nominal transformer tap, reactive shunts are included unlike in step-by-step Zbus where all these are neglected to simplify the computing procedure. Accuracy can be further increased by treating loads as constant impedance/admittance, just like in transient stability and Zbus loadflow. Converged bus volatge can be used in fault calculation instead of just 1.0 pu.

1) The constant admittance load at bus i is

where is the converged loadflow voltage
= g - jb for inductive load

2) The equations for an off-nominal tap transformer in equivalent pi is

where a= tap on i-side, =transformer impedance, =resulting shunts

3) Transmission line charging is

MVAR= total line charging mvaBASE= 100 mva base
The 2 in the denominator of Eq. (3) simply means that the Mvar charging is equally divided as in pi-model. When Mvar is divided by Base, it results to pu admittance.
4) Reactive shunt is

(-) for reactor, (+) for capacitor, there is no division by 2 here. The positive and zero sequence Ybus matrix has

elements called diagonal Y(i,i) and off-diagonal Y(i,j). For uncoupled zero sequence and positive sequence,

the second part of Eq. (6) is due to load, transformer shunt or line charging. It is also a summation if there are more of them. Eqs. 5 and 6 constitute to what is referred to as nodal admittance matrix. The matrix formed is called as Ybus.

5.2 Example Calculation.

The data is from section 4.2.
--- Eq. (7)
The diagonal elements are;
Y(1,1)= = 1/0.05 + 1/0.10 + 1/0.15 = 36.666
Y(2,2)= = 1/0.2 + 1/0.15 = 11.667
Y(3,3)== 1/0.2 + 1/0.10 = 15.0
Y(4,4)== 1/0.10 + 1/0.10 = 20.0

The off-diagonal elements are;
Y(1,2)== - 1/0.15 = -6.667
Y(3,4)== - 1/0.10 = -10.0,
Y(1,4)== - 1/0.10 = -10.0

Ybus is symmetric. Plugging the computed elements in matrix,

5.3 The Method of Singular Transformation (MST)

Ybus with coupling can be solved using Singular Transformation, see Ref 2. First we define some variables and terms. Branch and link are terminologies in graph theory.

1. n, nodes= or busses, junction of one or more elements. Reference bus is designated as node=0.
2. e, elements=lines, transformer or generator impedance or its equivalent circuits. Synonymous to branches in other topics, except here.
3. branch= element when added to the partial network does not close a loop
4. link= element when added to the partial network closes a loop
5. Partial Network (PN)=a connected graph of a network. All elements are not yet included.
6. Ā= element-node incidence matrix, shows the incidence of elements in a connected network.
Ā has a(i,j)= 1 if ith element is incident to and oriented away from jth node.
= -1 if ith element is incident to and oriented from jth node.
= 0 if not incident to jth node.
The definition above requires a direction of elements. It is similar to the definition in powerflow that a MW flow on a line away from a node is + and if towards the node -. For simplicity a(i,j)=1 for a "From Bus", and a(i,j)= -1 for a "To Bus". Here, the direction is dictated by p-q entry in the input data.

Example Formation of Ā and A
Ā = Element-node Incidence Matrix, and A= Bus Incidence Matrix
See Table 1 and Fig_1 on section 4.2 for the data used in this section.
When the reference bus 0 is deleted from the Ā (element-node incidence matrix whose size is row=6, columns=5), it becomes A (bus incidence matrix whose size is row=6, columns=4). Note that element 4-3 is in conformity with the coupling data. The from bus here is 4 (positive) and the to bus is negative. The shaded portion is the actual matrix. The A (bus incidence matrix) is shown below.

From Ohm's Law, the drawing on Fig_3 v=E-z*I

Reversing the directions, in impedance form,

In admittance form, dividing Eq. (9) by zPQ

and also from Eq. (9)
A set of unconnected elements is defined as primitive network. The performance equation of a primitive network can be derived from Eqs. (12) and (13) by expressing

[i] + [i'] = [y][v]    admittance form --- Eq. (14)
[v] + [e] = [z][i]    impedance form --- Eq. (15)
where [i] = all branch currents, including sources
   [z], [y] = are primitive impedance and admittance.
   z[i,i]= diagonal, is the self impedance
   z[i][j]= off-diagonal, is the mutual impedance of element i to element j.

Pre-multiply Eq. (14) by

where [A] is the bus incidence matrix.
But
, because of Kirchauff's Current Law, "The algebraic sum of currents in a bus is zero."

Eq. (17) refers to the algebraic sum of source currents at a bus, which is the injected current. From Eq. (16)

By power-invariant rule, the sum of powers in the primitive network is equal to the sum of powers in the interconnected network.
From Eq. (19) we need to take the conjugate (*) of [Ibus] then get its transposed. Using Eq. (17);

Eq. (20) is so because

Substituting Eq. (20) in Eq. (19)

Substitute Eq. (21) in Eq. (18) for [v];

--- Eq. (23)
5.7 Another Computer Technique to MST, Direct Method

We can visualize better how each element of the coupling matrix [z] affects Y(i,j) by considering the Ybus of Eq. (7) when line 1-2 is taken out of the matrix. The said z(1,2) here is the self impedance and the matrix does not contain mutual coupling. The figure below is not a real matrix addition as it must be compatible in matrix size. The figure is just a visual representation of how to include the effect of z(1,2) inside the Ybus matrix. By the figure, it says that z(1,2) is added to 4 elements -- Y(1,1), Y(1,2), Y(2,1), and Y(2,2) and taking care of whether it is positive or negative as shown by 1 or -1 inside the 2x2 matrix.
--- Eq. (25)
Shown below is the [y] and by observation, each element in [y] affects 4 elements in Ybus. This is similar to the figure in Eq. (25). This observation can improve steps 5,6, and 7 of section 5.5. As an example, y(1-2,4-3) will affect the following
1) Y(1,4) += y(1-2, 4-3)
2) Y(1,3) -= y(1-2, 4-3)
3) Y(2,4) -= y(1-2, 4-3)
4) Y(2,3) += y(1-2, 4-3)

When all the elements of [y] are taken outside the Ybus, the figure below shows how each element can be added to the matrix. This has a good analogy to Eq. (25).

The programming logic then becomes;

1. Build partial Ybus, those branches not included in a coupling group. The first matrix in Eq. (27). 2. Loop for all coupled groups. 3. Form the primitive impedance matrix [z] 4. Invert [z] to get [y] using Shipley Inversion 5. Each y(p-q,r-s) of [y] affects 4 Y(i,j)'s of Ybus. Perform the add/subtract. 6. Next coupled group.

Search the entire [y] and update 4 elements. We can not just search the upper triangle of [y] to perform the update because by inspection: Y(1,1) has the elements y(1-2,1-4) and y(1-4,1-2), also Y(4,4) has y(1-4,4-3) and y(4-3,1-4). So search the entire [y]. To implement this procedure for the upper triangle portion of Y(), see to it that the 4 elemnets Y(i,j) solved has P<=Q, P<=S, Q<=R and Q<=S and use all elements of [y]. This method is also true for sparse matrix application. Note that for sparse matrix, only the upper triangle is stored and P,Q OR Q,P can be accessed using start[P] or start[Q]. When the program uses node ordering, upper triangle is based on this node ordering. See the topic when a node is left-side or right-side of a diagonal.

Code Snippet for Full Matrix with upper triangle update only for (i=1; i< z_size; i++) {lp = zP[i]; lq = zQ[i]; for (j=1; j<= z_size; j++) {lr = zP[i]; ls = zQ[j]; //there is no zR[ ], zS[ ] here if (lp<=lr) YZ[lp][lq] += cmat[i][j]; if (lp<=ls) YZ[lp][ls] -= cmat[i][j]; if (lq<=lr) YZ[lq][lr] -= cmat[i][j]; if (lq<=ls) YZ[lq][ls] += cmat[i][j]; } } Code Snippet for Sparse Matrix without Node Ordering Y(p,r)=Y(p,r) + y(pq,rs) Y(p,s)=Y(p,s) - y(pq,rs) ^ ^ ^ ^ ^ ^ ^ ^ Y(q,r)=Y(q,r) - y(pq,rs) Y(q,s)=Y(q,s) + y(pq,rs) ^ ^ ^ ^ ^ ^ ^ ^ for (i=1; i< z_size; i++) {lp = zP[i]; lq = zQ[i]; for (j=1; j<= z_size; j++) {lr = zP[i]; ls = zQ[j]; if (lp==lr) YZd[lp] += cmat[i][j]; //diagonal else if (lp< lr) YZ[BrnSrNsrt(lp,lr)] += cmat[i][j]; if (lp==ls) YZd[lp] += cmat[i][j]; //diagonal else if (lp< ls) YZ[BrnSrNsrt(lp,ls)] += cmat[i][j]; if (lq==lr) YZd[lq] += cmat[i][j]; //diagonal else if (lq< lr) YZ[BrnSrNsrt(lq,lr)] += cmat[i][j]; if (lq==ls) YZd[lq] += cmat[i][j]; //diagonal else if (lq< ls) YZ[BrnSrNsrt(lq,ls)] += cmat[i][j]; } }

Program code for Singular Transformation

1: //Program: Method of Singular Transformation file=singular.c
   2: //Copyright 1996. All Rights Reserved   June 9,1996
   3: //Joesel Pante
   4: //Program to perform At*[y]*A, used to verify other programs
   5: //    on direct coupling compensation
   6: #define LMAX 25
   7: #define BMAX 20
   8: #include <stdio.h>
   9: int nb=0,nl=0;
  10: int kp[LMAX],kq[LMAX],ckt[LMAX],At[LMAX][LMAX],A[LMAX][LMAX];
  11: float Y[BMAX][LMAX],Ytemp[BMAX][LMAX],yprim[LMAX][LMAX],x1[LMAX],x0[LMAX];
  12: FILE *infile;
  13: FILE *outfile;  
  14: char source_name[30];
  15: 
  16: void print_A(int A[][LMAX],int row, int col)
  17: {int i,j;
  18:  printf("     ");
  19:  for (i=1;i<=col;i++) printf(" %3d",i); 
  20:  printf("\n");
  21:  for (i=1;i<=row;i++)
  22:   {printf("%2d %2d",kp[i],kq[i]);
  23:    for (j=1;j<=col;j++) printf(" %3d",A[i][j]);
  24:       printf("\n");
  25:   }
  26: }
  27: 
  28: void print_At(int A[][LMAX],int row, int col)
  29: {int i,j;
  30:  printf(" ");
  31:  for (i=1;i<=col;i++) printf(" %3d",kp[i]); 
  32:  printf("\n");
  33:  printf(" ");
  34:  for (i=1;i<=col;i++) printf(" %3d",kq[i]); 
  35:  printf("\n");
  36:  printf(" ");
  37:  for (i=1;i<=col;i++) printf(" %3d",ckt[i]); 
  38:  printf("\n");
  39:   
  40:  for (i=1;i<=row;i++)
  41:   {printf("%2d",i);
  42:    for (j=1;j<=col;j++) printf(" %3d",A[i][j]);
  43:       printf("\n");
  44:   }
  45: }                 
  46: 
  47: void print_zy(float zy[][LMAX],int row,int col)
  48: {int i,j;
  49:  for (i=1;i<=row;i++) printf(" %8d",kp[i]); 
  50:  printf("\n");
  51:  for (i=1;i<=col;i++) printf(" %8d",kq[i]); 
  52:  printf("\n");
  53:  for (i=1;i<=col;i++) printf(" %8d",ckt[i]); 
  54:  printf("\n");
  55:   
  56:  for (i=1;i<=row;i++)
  57:   {for (j=1;j<=col;j++) printf(" %8.3f",zy[i][j]);
  58:    printf("\n");
  59:   }
  60: } 
  61: 
  62: void print(float a[][LMAX],int row,int col)
  63: {int i,j;
  64:  for (i=1;i<=row;i++)
  65:   {for (j=1;j<=col;j++) printf(" %8.3f",a[i][j]);
  66:    printf("\n");
  67:   }
  68: } 
  69: 
  70: void Enter(void)
  71: {char *ch;printf("Enter to continue\n");gets(ch);}
  72: 
  73: int abso(int value)
  74: {if (value<0) return (-value);
  75:    else return(value);
  76: }
  77:     
  78: void shipley(void)
  79: //Shipley Inversion Routine
  80: {int i,j,k;
  81:  for (k=1;k<=nl;k++)
  82:   {yprim[k][k]= -1/yprim[k][k];
  83:    for (i=1;i<=nl;i++)
  84:      if (i!=k) 
  85:         yprim[i][k]*=yprim[k][k];
  86:    for (i=1;i<=nl;i++)
  87:      if (i!=k) 
  88:        for (j=1;j<=nl;j++)
  89:          if (j!=k) 
  90:             yprim[i][j]+=yprim[i][k]*yprim[k][j];
  91:    for (i=1;i<=nl;i++)
  92:      if (i!=k) 
  93:        yprim[k][i]*=yprim[k][k];
  94:   }
  95:  for (i=1;i<=nl;i++)    //reverse sign
  96:    for (j=1;j<=nl;j++)
  97:       yprim[i][j]=-yprim[i][j];
  98: }
  99: 
 100: int Find_pq(int lp,int lq,int cktn)
 101: //Searches lp->lq in [z]. Negative if lq->lp
 102: {int i;
 103:  for (i=1;i<=nl;i++)
 104:    {if (kp[i]==lp && kq[i]==lq && ckt[i]==cktn) return(i);
 105:     if (kp[i]==lq && kq[i]==lp && ckt[i]==cktn) return(-i);
 106:    }
 107:  printf("ERROR on [z].search.stop\n");
 108:  Enter(); return(0);
 109: }
 110: 
 111: void input(void)
 112: {int p,q,ckt1,r,s,ckt2,lrr,lr,lss,ls;
 113:  float xm;
 114:  printf("i:   P   Q   ckt   x1   x0\n");
 115:  while(1)
 116:   {nl++;
 117:    fscanf(infile,"%d %d %d %f %f",
 118:      kp+nl,kq+nl,ckt+nl,x1+nl,x0+nl);
 119:    if (kp[nl]==99) break;
 120:    if (kp[nl]>nb) nb=kp[nl];   //get largest bus
 121:    if (kq[nl]>nb) nb=kq[nl]; 
 122:    printf("%2d %2d %2d %2d %5.2f %5.2f\n",
 123:            nl,kp[nl],kq[nl],ckt[nl],x1[nl],x0[nl]);
 124: //from=+1, to=-1, yprim[diag]=self_z 
 125:    A[nl][kp[nl]]=1; A[nl][kq[nl]]=-1; yprim[nl][nl]=x0[nl];
 126:   }
 127:  nl--; //exclude terminator 99
 128:  Enter(); 
 129:  printf("P  Q  ckt   R  S  ckt   xm\n");
 130:  while(1)
 131:   {fscanf(infile,"%d %d %d %d %d %d %f",
 132:       &p,&q,&ckt1,&r,&s,&ckt2,&xm);
 133:    if (p==99) break;
 134:    printf("%2d %2d %2d   %2d %2d %2d %5.2f\n",
 135:            p,q,ckt1,r,s,ckt2,xm); 
 136: //lrr*lss<0 is the same as (lrr<0 ^ lss<0)-exclusive or
 137:    lrr=Find_pq(p,q,ckt1); lr=abso(lrr);
 138:    lss=Find_pq(r,s,ckt2); ls=abso(lss);
 139:    if (lrr*lss<0) 
 140:      {yprim[lr][ls]=-xm; yprim[ls][lr]=-xm;}
 141:    else 
 142:      {yprim[lr][ls]=xm; yprim[ls][lr]=xm;}
 143:   }
 144:   Enter();
 145: }
 146: 
 147: void main(void)
 148: {int i,j,k;
 149: // printf("Enter input file:");gets(source_name);
 150: //infile =fopen(source_name,"r"); 
 151:  infile=fopen("scm.dat","r");
 152:  outfile=fopen("output.txt","w");
 153:  input(); 
 154:  printf("nl=%d nb=%d\n",nl,nb);
 155:  printf(" Bus Incidence Matrix A:\n");
 156:  print_A(A,nl,nb);Enter();
 157: //get A(transposed)
 158:  for (i=1;i<=nb;i++)
 159:    for (j=1;j<=nl;j++)
 160:      At[i][j]=A[j][i];
 161:  printf("Transposed of A:\n");
 162:  print_At(At,nb,nl);Enter();
 163:  printf("Primitive Impedance Matrix [z]:\n"); 
 164:  print_zy(yprim,nl,nl);
 165:  shipley();Enter();         
 166:  printf("Primitive Admittance Matrix [y]:\n"); 
 167:  print_zy(yprim,nl,nl);
 168:  printf("A(transposed)*[y]:\n");
 169:  for (i=1;i<=nb;i++)
 170:   for (j=1;j<=nl;j++)
 171:    for (k=1;k<=nl;k++)
 172:      Ytemp[i][j]+=At[i][k]*yprim[k][j];
 173:  print(Ytemp,nb,nl);
 174:  printf("Resulting A(t)*[y]*A:\n");
 175:  for (i=1;i<=nb;i++)
 176:   for (j=1;j<=nb;j++)
 177:    for (k=1;k<=nl;k++)
 178:      Y[i][j]+=Ytemp[i][k]*A[k][j];
 179:  print(Y,nb,nb);
 180:  i=fclose(infile);i=fclose(outfile);
 181: }

INPUT/OUTPUT
0  1  1   0.05    0.1
0  2  1   0.2     0.5
0  3  1   0.2     0.5
1  2  1   0.15    0.4
1  4  1   0.1     0.2
4  3  1   0.1     0.3
99 0  0   0       0
1  2  1    1  4  1     0.10
1  2  1    4  3  1     0.20
99 0  0    0  0  0     0 
98 0  0    0  0  0     0
    
OUTPUT PART

i:   P   Q   ckt   x1   x0                               
 1  0  1  1  0.05  0.10                                  
 2  0  2  1  0.20  0.50                                  
 3  0  3  1  0.20  0.50                                  
 4  1  2  1  0.15  0.40                                  
 5  1  4  1  0.10  0.20                                  
 6  4  3  1  0.10  0.30                                  
Enter to continue                                        
                                                         
P  Q  ckt   R  S  ckt   xm                               
 1  2  1    1  4  1  0.10                                
 1  2  1    4  3  1  0.20                                
Enter to continue                                        
                                                         
nl=6 nb=4                                                
Enter to continue

Bus Incidence Matrix A:                                
        1   2   3   4                                   
 0  1  -1   0   0   0                                   
 0  2   0  -1   0   0                                   
 0  3   0   0  -1   0                                   
 1  2   1  -1   0   0                                   
 1  4   1   0   0  -1                                   
 4  3   0   0  -1   1                                   
Enter to continue                                       
                                                        
Transposed of A:                                        
    0   0   0   1   1   4                               
    1   2   3   2   4   3                               
    1   1   1   1   1   1                               
 1  -1   0   0   1   1   0                              
 2   0  -1   0  -1   0   0                              
 3   0   0  -1   0   0  -1                              
 4   0   0   0   0  -1   1                              
Enter to continue                                       
                                                        
Primitive Impedance Matrix [z]:                         
        0        0        0        1        1        4  
        1        2        3        2        4        3  
        1        1        1        1        1        1  
    0.100    0.000    0.000    0.000    0.000    0.000  
    0.000    0.500    0.000    0.000    0.000    0.000  
    0.000    0.000    0.500    0.000    0.000    0.000  
    0.000    0.000    0.000    0.400    0.100    0.200  
    0.000    0.000    0.000    0.100    0.200    0.000  
    0.000    0.000    0.000    0.200    0.000    0.300  
Enter to continue                                       
                                                        
Primitive Admittance Matrix [y]:                        
        0        0        0        1        1        4  
        1        2        3        2        4        3  
        1        1        1        1        1        1  
   10.000    0.000    0.000    0.000    0.000    0.000  
    0.000    2.000    0.000    0.000    0.000    0.000  
    0.000    0.000    2.000    0.000    0.000    0.000  
    0.000    0.000    0.000    4.615   -2.308   -3.077  
    0.000    0.000    0.000   -2.308    6.154    1.538  
    0.000    0.000    0.000   -3.077    1.538    5.385  
A(transposed)*[y]:                                      
  -10.000    0.000    0.000    2.308    3.846   -1.538  
    0.000   -2.000    0.000   -4.615    2.308    3.077  
    0.000    0.000   -2.000    3.077   -1.538   -5.385  
    0.000    0.000    0.000   -0.769   -4.615    3.846  
Resulting A(t)*[y]*A:                                   
   16.154   -2.308    1.538   -5.385                    
   -2.308    6.615   -3.077    0.769                    
    1.538   -3.077    7.385   -3.846                    
   -5.385    0.769   -3.846    8.462

Some Code Explanation

This program could not be used for large data because of matrices involved. This is created only to check the validity of the next computer technique, the Direct Method. The Bus incidence matrix A is formed rightaway while reading the line data in line_123. The diagonal of the primitive impedance matrix [z] is also formed. While reading the coupling data on line_128, the off-diagonal of [z] is also filled up in lines_133 to 137. Proper sign or direction of coupling is taken cared of depending its entry in pq-rs vectors. The matrix At[][] is used to store the transpose of A in line_154 of the main program. In line_159, Shipley is called to invert [z]. Only one matrix is used as this is an in-place inversion. The product At*[y] is stored in Ytemp[][] in line_166. The final product is done in line_172 and Ybus is stored in Y[][].

Program Code for MST, Full Matrix, Upper Triangle only

1: /* Program Direct Method: MST best implemented 
   2:    Copyright 1996. All Rights Reserved July 1996
   3:    Joesel Pante.  BC, Canada.
   4: 
   5:    The constraints are:
   6:      1. contiguous numbering
   7:      2. X only ,not complex
   8:      3. upper triangle only for Ybus matrix
   9:    Compared to Step-by-step Zbus Formation, MST is
  10:      1. a direct method, less memory required for data 
  11:      2. thousandfold faster when implemented on large systems
  12:      3. amenable to sparsity programming
  13:      4. coupling on parallel lines allowed
  14: */
  15: #define LMAX 35
  16: #define BMAX 20
  17: #define GMAX 5
  18: #define CMAX 10
  19: #include <stdio.h>
  20: 
  21: void input(void);                                   
  22: void primitive(int k);
  23: void IncludeDiag(int lp,int lq,int cktn,float xm1); 
  24: void IncludeOffDiag(void);
  25: int Find_pq(int lp,int lq,int cktn);                
  26: 
  27: int nb=0,nl=0,debug=1, kp[LMAX],kq[LMAX],ckt[LMAX];
  28: float YP[BMAX][BMAX],YZ[BMAX][BMAX],x1[LMAX],x0[LMAX];
  29: 
  30: int mg=1,GrStart[GMAX],GrEnd[GMAX]; //coupled group index
  31: int mc=0,p[CMAX],q[CMAX],ckt1[CMAX],r[CMAX],s[CMAX],ckt2[CMAX],
  32:     pq[CMAX],rs[CMAX];
  33: int z_size;         /*cmat(z[]/y[]=primitive coupling matrix) size*/
  34: int mut_yn[LMAX],   zP[10],zQ[10],zC[10];  
  35: /* 10 coupled lines in a group, this is different from number of groups*/
  36: float xm[CMAX],     cmat[10][10];
  37: FILE *infile;
  38: FILE *outfile;
  39: char source_name[30]; //input data filename
  40: /* YP/YZ Ybus positive & zero sequence
  41:    x1/x0 positive & zero sequence impedance
  42: 
  43:    nl=number of lines, nb=number of busses,
  44:    mg=coupled groups,  mc=coupling pairs
  45:    GrStart/GrEnd-pointers to pq-rs for group start/end
  46:    pq[]/rs[] -pointers to self impedance
  47:    mut_yn[]  -line data pointer to coupling group
  48:    zP[]/zQ[] -pq,rs of [z]
  49:    cmat[][] -[z],[y], the primitive impedance/admittance matrix
  50: */
  51: void Enter(void)
  52: {char *ch; printf("Enter to continue\n"); gets(ch);}
  53: 
  54: void print(float A[][BMAX])
  55: {int i,j;
  56:  for (i=1;i<=nb;i++) 
  57:   printf("%10d",i);
  58:  printf("\n");
  59:  for (i=1;i<=nb;i++)
  60:   {for (j=1;j<=nb;j++)
  61:      printf("%10.6f",A[i][j]);
  62:    printf("\n");
  63:   }   
  64:  printf("\n");
  65:  Enter();
  66: }
  67: 
  68: void print_zy()
  69: {int i,j;
  70:  for (i=1;i<=z_size;i++) printf("%10d",zP[i]); printf("\n");
  71:  for (i=1;i<=z_size;i++) printf("%10d",zQ[i]); printf("\n");
  72:  for (i=1;i<=z_size;i++) printf("%10d",zC[i]); printf("\n"); 
  73:  
  74:  for (i=1;i<=z_size;i++)
  75:   {for (j=1;j<=z_size;j++) 
  76:     printf("%10.6f",cmat[i][j]);
  77:    printf("\n");
  78:   }  
  79:  printf("\n");
  80:  Enter();
  81: }
  82: 
  83: int abso(int value)
  84: {if (value<0)  return(-value);
  85:         else   return(value);
  86: }
  87:   
  88: 
  89: void input(void)
  90: {int i,start=1;
  91:  float r1,r0;
  92:   printf(" i: P   Q  ckt x1   x0\n");      
  93:   while (1)
  94:    {nl++;
  95:     fscanf(infile,"%d %d %d %f %f %f %f",
  96:       kp+nl,kq+nl,ckt+nl,&r1,x1+nl,&r0,x0+nl);
  97:     if (kp[nl]==99) break;
  98:     printf("%2d: %2d %2d  %2d %5.2f %5.2f\n",
  99:            nl,kp[nl],kq[nl],ckt[nl],x1[nl],x0[nl]); 
 100:     if (kp[nl]>nb) nb=kp[nl]; if (kq[nl]>nb) nb=kq[nl];
 101:    }
 102:    nl--;
 103:   Enter(); 
 104:   printf(" i: P  Q  ckt- R  S  ckt  xm\n");
 105:   while (1)
 106:    {mc++;
 107:     fscanf(infile,"%d %d %d %d %d %d %f %f",
 108:                 p+mc,q+mc,ckt1+mc,r+mc,s+mc,ckt2+mc,&r1,xm+mc);
 109:     if (p[mc]==98) break;                            /*end of data*/
 110: /* Note that mg=1 initially.This is used for the mut_yn pointer to group=1*/
 111:     if (p[mc]==99) {GrStart[mg]=start;GrEnd[mg++]=mc-1;start=mc+1;}
 112: /*include group terminator 99 to simulate the database editor for
 113:   full blown short circuit software
 114: */
 115:     printf("%2d: %2d %2d %2d  %2d %2d %2d %5.2f\n",
 116:             mc,p[mc],q[mc],ckt1[mc],r[mc],s[mc],ckt2[mc],xm[mc]);
 117:     if (p[mc]!=99)  
 118:      for (i=1;i<=nl;i++)              /*code pointers vice-versa*/
 119:       {if (kp[i]==p[mc] && kq[i]==q[mc] && ckt[i]==ckt1[mc]) {pq[mc]=i;mut_yn[i]=mg;}
 120:        if (kp[i]==q[mc] && kq[i]==p[mc] && ckt[i]==ckt1[mc]) {pq[mc]=i;mut_yn[i]=mg;}
 121:        if (kp[i]==r[mc] && kq[i]==s[mc] && ckt[i]==ckt2[mc]) {rs[mc]=i;mut_yn[i]=mg;}
 122:        if (kp[i]==s[mc] && kq[i]==r[mc] && ckt[i]==ckt2[mc]) {rs[mc]=i;mut_yn[i]=mg;}
 123:       }
 124:    } /*while*/
 125:    mg--;mc--;
 126:    Enter(); 
 127:    printf("Group Start  End\n");
 128:    for (i=1;i<=mg;i++) 
 129:      printf("%2d   %2d    %2d\n",i,GrStart[i],GrEnd[i]);
 130:    printf(" i  pq   P   Q   R   S   rs\n");  
 131:    for (i=1;i<=mc;i++)
 132:     printf("%2d %2d  %2d  %2d  %2d  %2d  %2d\n",
 133:              i,pq[i],p[i],q[i],r[i],s[i],rs[i]);
 134: }
 135:  
 136: void primitive(int k)                  /*forms [z],inverts to [y]*/
 137: {int i,j,lr,lrr,ls,lss;                /*k-group number*/
 138:  float xm1;                            /*intialize [z] size*/
 139:  z_size=0;                             
 140:  for (i=GrStart[k];i<=GrEnd[k];i++)
 141:   {xm1=x0[pq[i]]; IncludeDiag(p[i],q[i],ckt1[i],xm1); /*include p-q in [z]*/
 142:    xm1=x0[rs[i]]; IncludeDiag(r[i],s[i],ckt2[i],xm1); /*include r-s in [z]*/
 143:   }
 144:  for (i=1;i<=z_size;i++)                 /*zero-out [z] off-diag*/
 145:    for (j=1;j<=z_size;j++) if (i!=j) cmat[i][j]=(float)0.0;
 146:  for (i=GrStart[k];i<=GrEnd[k];i++)      /*include [z] off-diag*/
 147: /*consider direction of coupling */
 148: /* lrr*lss<0 is the same as (lrr<0 ^ lss<0) -exclusive or */
 149:   {lrr=Find_pq(p[i],q[i],ckt1[i]); lr=abso(lrr);
 150:    lss=Find_pq(r[i],s[i],ckt2[i]); ls=abso(lss);
 151:    if (lrr*lss<0) {cmat[lr][ls]=-xm[i]; cmat[ls][lr]=-xm[i];}
 152:            else   {cmat[lr][ls]=xm[i]; cmat[ls][lr]=xm[i];}
 153:   }
 154:  if (debug) 
 155:    {printf("primitive [z]\n");print_zy();}
 156: /**********************/  
 157: /*Shipley_Inversion();*/
 158: /**********************/
 159:  for (k=1;k<=z_size;k++)                 
 160:   {cmat[k][k]= -1/cmat[k][k];           
 161:    for (i=1;i<=z_size;i++)              
 162:     if (i!=k) 
 163:         cmat[i][k]*=cmat[k][k];
 164:    for (i=1;i<=z_size;i++)
 165:      if (i!=k) 
 166:        for (j=1;j<=z_size;j++)
 167:           if (j!=k) 
 168:              cmat[i][j]+=cmat[i][k]*cmat[k][j];
 169:    for (i=1;i<=z_size;i++)
 170:      if (i!=k) 
 171:        cmat[k][i]*=cmat[k][k];
 172:   }
 173:   for (i=1;i<=z_size;i++)        /*reverse sign*/
 174:      for (j=1;j<=z_size;j++)
 175:         cmat[i][j]=-cmat[i][j];
 176:  if (debug) 
 177:    {printf("After Inversion [y]\n");print_zy();}
 178: }
 179: 
 180: /* find lp-lq in [z] whose size is z_size
 181:          zP   1     1     4
 182:          zQ   2     4     3
 183:       zP zQ +----+----+----+
 184:        1  2 | x  |    |    |
 185:             +----+----+----+
 186:  [z] = 1  4 |    |  x |    |
 187:             +----+----+----+
 188:        4  3 |    |    |  x |
 189:             +----+----+----+           */
 190: void IncludeDiag(int lp,int lq,int cktn,float xm1)
 191: {int i;                                                  
 192:  if (z_size)                                             
 193:    for (i=1;i<=z_size;i++)                               
 194:     { if (zP[i]==lp && zQ[i]==lq && zC[i]==cktn) return;
 195:       if (zP[i]==lq && zQ[i]==lp && zC[i]==cktn) return; 
 196:     }                                                    
 197:  zP[++z_size]=lp; zQ[z_size]=lq; zC[z_size]=cktn;        
 198:  cmat[z_size][z_size]=xm1;                               
 199: }                                                        
 200:                                            
 201: int Find_pq(int lp,int lq,int cktn)
 202: {int i;                     /*searches lp->lq in [z]. Negative if lp->lq*/
 203:  for (i=1;i<=z_size;i++)
 204:   {if (zP[i]==lp && zQ[i]==lq && zC[i]==cktn) {return(i);}
 205:    if (zP[i]==lq && zQ[i]==lp && zC[i]==cktn) {return(-i);}
 206:   }
 207:  printf("ERROR on [z] search.stop\n");
 208:  Enter();
 209:  return(0);
 210: }
 211:     
 212: /****************/
 213: /*    MAIN      */
 214: /****************/
 215: void main(void)
 216: {int i,j,k,lp,lq,lr,ls;
 217:  float a; 
 218:  infile=fopen("scm.dat","r");  
 219:  outfile=fopen("output.txt","w");
 220:  input(); 
 221:  printf("nl=%d nb=%d mg=%d mc=%d\n",nl,nb,mg,mc);
 222: //both Pos & Zer sequence Ybus will be formed
 223:  for (i=0;i<=nb;i++)
 224:    for (j=0;j<=nb;j++)       /*initialize Ybus,Zbus*/
 225:      {YP[i][j]=0;YZ[i][j]=0;}
 226: 
 227: /*form Ybus, Yij=-1/z(i,j),  Yii=sum 1/z(i,j), Yjj=sum 1/z(j,i) */
 228: /*Note: for parallel lines the equivalent y_eq=sum of y's
 229:   Note the meaning of YP[i][j]=-a  and YP[i][j]-=a.
 230:   The YP[i][j]=-a  is the previous definition of Yij, while
 231:   the YP[i][j]-=a  is y_eq=sum y's
 232: */
 233: 
 234:  for (k=1;k<=nl;k++)
 235:    {i=kp[k];j=kq[k];a=1/x1[k];  YP[i][i]+=a; YP[j][j]+=a; /*positive seq*/
 236:     if (i<j) 
 237:       YP[i][j]-=a; 
 238:     else 
 239:       YP[j][i]-=a;
 240:     if (!mut_yn[k]) 
 241:      {a=1/x0[k]; YZ[i][i]+=a; YZ[j][j]+=a; /*uncoupled zero seq*/
 242:       if (i<j) 
 243:         YZ[i][j]-=a; 
 244:       else 
 245:         YZ[j][i]-=a;
 246:      }
 247:    }
 248:  printf("Completed Pos Seq Ybus:\n"); 
 249:  print(YP);
 250:  printf("Partial Zero Seq Ybus:\n");     
 251:  print(YZ);
 252: 
 253: /*Include the effect of coupling*/
 254:  for (k=1;k<=mg;k++)        /*loop for coupled groups*/
 255:   {primitive(k);            /*for each group form [z],size=z_size*/
 256: /* Update Ybus. For full matrix programming,
 257:    it can be updated using,
 258:    Y(p,r)=Y(p,r) + y(pq,rs)   Y(p,s)=Y(p,s) - y(pq,rs)
 259:             ^ ^      ^  ^              ^ ^      ^   ^
 260:    Y(q,r)=Y(q,r) - y(pq,rs)   Y(q,s)=Y(q,s) + y(pq,rs)
 261:             ^ ^       ^ ^              ^ ^       ^  ^
 262:    Update only the upper triangle, i.e. i<=j.
 263:    Note that the entire [y] is used and test for i<=j is done on 
 264:    p<=r
 265:    p<=s
 266:    q<=r
 267:    q<=s
 268:    Otherwise, without these tests, even the lower triangle gets solved. 
 269: */
 270:    for (i=1;i<=z_size;i++)
 271:      {lp=zP[i];lq=zQ[i];
 272:       for (j=1;j<=z_size;j++)
 273:        {lr=zP[j];ls=zQ[j];
 274:         if (lp<=lr) YZ[lp][lr]+=cmat[i][j]; 
 275:         if (lp<=ls) YZ[lp][ls]-=cmat[i][j];
 276:         if (lq<=lr) YZ[lq][lr]-=cmat[i][j]; 
 277:         if (lq<=ls) YZ[lq][ls]+=cmat[i][j];
 278:        }
 279:       }
 280:   printf("After Updating YZ[][]\n"); 
 281:   print(YZ);
 282:   } /*for k*/
 283:  return;
 284: }

INPUT DATA
0  1  1 0  0.05  0  0.1
0  2  1 0  0.2   0  0.5
0  3  1 0  0.2   0  0.5
1  2  1 0  0.15  0  0.4
1  4  1 0  0.1   0  0.2
4  3  1 0  0.1   0  0.3
99 0  0 0  0     0  0
1  2  1    1  4  1   0  0.10
1  2  1    4  3  1   0  0.20
99 0  0    0  0  0   0  0 
98 0  0    0  0  0   0  0
  
 OUTPUT PART                                                   
 i: P   Q  ckt x1   x0                          
 1:  0  1   1  0.05  0.10                       
 2:  0  2   1  0.20  0.50                       
 3:  0  3   1  0.20  0.50                       
 4:  1  2   1  0.15  0.40                       
 5:  1  4   1  0.10  0.20                       
 6:  4  3   1  0.10  0.30                       
Enter to continue                               
                                                
 i: P  Q  ckt- R  S  ckt  xm                    
 1:  1  2  1   1  4  1  0.10                    
 2:  1  2  1   4  3  1  0.20                    
 3: 99  0  0   0  0  0  0.00                    
Enter to continue                               
                                                
Group Start  End                                
 1    1     2                                   
 i  pq   P   Q   R   S   rs                     
 1  4   1   2   1   4   5                       
 2  4   1   2   4   3   6                       
 3  0  99   0   0   0   0                       
nl=6 nb=4 mg=1 mc=3                             
Completed Pos Seq Ybus:                         
         1         2         3         4        
 36.666664 -6.666667  0.000000-10.000000        
  0.000000 11.666666  0.000000  0.000000        
  0.000000  0.000000 15.000000-10.000000        
  0.000000  0.000000  0.000000 20.000000        
                                               
Enter to continue                               
                                                
Partial Zero Seq Ybus:                          
         1         2         3         4        
 10.000000  0.000000  0.000000  0.000000        
  0.000000  2.000000  0.000000  0.000000        
  0.000000  0.000000  2.000000  0.000000        
  0.000000  0.000000  0.000000  0.000000        
                                               
Enter to continue                               
                                                
primitive [z]                                   
         1         1         4                  
         2         4         3                  
         1         1         1                  
  0.400000  0.100000  0.200000                  
  0.100000  0.200000  0.000000                  
  0.200000  0.000000  0.300000                  
                                               
Enter to continue                               
                                                
After Inversion [y]                             
         1         1         4                  
         2         4         3                  
         1         1         1                  
  4.615385 -2.307692 -3.076923                  
 -2.307692  6.153846  1.538461                  
 -3.076923  1.538461  5.384615                  
                                               
Enter to continue                               
                                                
After Updating YZ[][]                           
         1         2         3         4        
 16.153847 -2.307692  1.538461 -5.384615        
  0.000000  6.615385 -3.076923  0.769231        
  0.000000  0.000000  7.384615 -3.846153        
  0.000000  0.000000  0.000000  8.461538

HOME
CHAPTER ONE: INTRODUCTION
CHAPTER TWO: PROGRAMMING TOPICS
CHAPTER THREE: MATRIX INVERSION ROUTINE
CHAPTER FOUR: STEP-BY-STEP ALGORITHMS
CHAPTER FOUR: Part 2
CHAPTER FIVE: DIRECT METHODS IN Ybus FORMATION
CHAPTER SIX: SPARSITY
CHAPTER SEVEN: MATRIX FACTORIZATION AND OPTIMAL ORDERING
CHAPTER EIGHT: SPARSE Zbus FROM FACTORED SPARSE Ybus
CHAPTER NINE: FAULT CALCULATIONS
REFERENCES

i n n o v a t e o r s t a g n a t e