Newton-Raphson

CHAPTER 4: Newton-Raphson Loadflow (NRLF)

Unlike Gauss-Siedel, (GS) which updates the node voltage one at a time, Newton Raphson, (NR) solves a voltage correction for all the nodes and update them. A simultaneous linear equation in the form of Ax=b is evaluated each iteration. Look again on "Program Listing_1: DC Circuit Flow" on how to solve Ax=b by Gauss-Jordan Reduction. Comparing NR with GS, GS has problems when the system becomes large. One reason is the presence of negative impedance as a result of transformer 3-winding representation. GS tends to increase in iteration count and is slow in computer time.

Let us solve a sample 2 equations in 2 unknowns and non-linear.
NRLF_Eq. (1)
Set an initial estimate, get a correction, update the estimate. By Taylor Series Expansion;
NRLF_Eq. (2)
can be expanded to
NRLF_Eq. (3)
Eq. (3) can be expanded to 3,4,etc as the exponent. The last 2 terms of the 2 equations have x in exponent 2 and the derivatives are in second order. By neglecting these last 2 terms, we get only the linear partial derivatives. This is just an ordinary derivative.
NRLF_Eq. (4)
Eq. (4) can be written in matrix form as shown on the left with the red arrow. This is in the form of Ax=b. Instead of solving the simultaneous non-linear equations, we are now solving simultaneous linear equations.
NRLF_Eq. (5)

Performing the derivatives on Eq. (1), yields the elements of the Jacobian matrix. The error vector is the difference between the scheduled value and the computed value when the state variables are substituted to the set of equations. The unknown is the correction to the state variables.
NRLF_Eq. (6)

Eq. (7) is the matrix representation of Eq. (4).
NRLF_Eq. (7)

Iteration=1:
Let x(1)=2,x(2)=2 be the initial estimate. Substituting to Eq. (7), and solving for the correction and updating x(1),x(2) yields;
NRLF_Eq. (8)
Iteration=2:
Substitute again the updated x(1),x(2) to Eq. (7), and solving for the correction and updating x(1),x(2) yields;
NRLF_Eq. (9)
This is to be continued until the error vector is less than some tolerance.

Derivation of NRLF Formula
This equation was taken from GSLF_Eq. (20). From the above 2 equation sample, this is similar to Eq. (1).
NRLF_Eq. (10) bus injected power
Comparing to Eq. (1), the numbers 2 and -16 are the actual values when x(1) and x(2) are substituted on the equations. Similarly, when the computed voltages are applied to Eq. (10), it should result to the scheduled value of generated power less load. But because there is some degree of error, then there is some mismatch. The error vector then becomes the scheduled injected power less than the computed injected power. The state variables are the voltage and the voltage angle, so the correction vector will be delta theta and delta e. If there are 2 nodes, there will be 4 state variables. So the size of the Jacobian matrix is 4x4.
NRLF_Eq. (11)
Eq. (11) can be re-written as Eq. (12)
NRLF_Eq. (12)
We now derive the elements of each submatrices H,J,N and L.

NRLF_Eq. (13)
Eq. (10) includes when j=i. It means that the diagonal of the Ybus is included. When taking the derivative of variable-i with respect to variable-j, there is only one node-j in the summation symbol. So the resulting derivative has no summation symbol. Whereas taking the derivative of i with respect to i, the summation symbol is included. There exists a product term of i to j. At the bottom of the summation, it says j is not equal to i. Taking a derivative with respect to theta when j=i simply is non existant because theta is zero. However, when taking the derivative with respect to E, there is a value even if theta is zero because cosine of zero is 1, while its sine is zero.

Some Improvement to NRLF
By multiplying N and L with E and dividing delta E with E, computing time can be reduced. The code is simplified too.
NRLF_Eq. (14)
The red numbers from 1 to 9 are for explanation purposes only. In line 2 the one term of 2EiEiGii is transferred inside the summation to satisfy when j=i. The summation portion now becomes the computed injected real power. Also the summation portion is equal to N'ij as shown in line 5. This is also true to line 7. One part of -2EiEiBii is thrown inside the summation which is equal to the computed reactive power injection.
NRLF_Eq. (15)

SUMMARY OF FORMULA FOR H,J,N',L' The prime in N and L is not Hij=EiEj(GijsinOij - BijcosOij) shown in the code. Hii=-sumHij, j.ne.i note: O is for theta Nij=EiEj(GijcosOij + BijsinOij) Nii=2EiEiGii + sum Nij, j.ne.i Jij=-Nij Jii=Nii - 2EiEiGii Lij=Hij Lii=-Hii - 2EiEiBii deltaPi=Pi(sch)-Nii +EiEiGii delta is for error deltaQi=Qi(sch)+Hii +EiEiBii sch=generated-load COMPUTER IMPLEMENTATION loop i loop j 1. Hij=EiEj(GijsinOij - BijcosOij) 2. Nij=EiEj(GijcosOij + BijsinOij) 3. Jij=-Nij 4. Lij=Hij 5. Hii-=Hij this is a C/C++ code 6. Nii+=Nij next j 7. Jii=Nii //Nii is not complete yet 8. Nii+=2EiEiGii 9. Lii= -Hii -2EiEiBii 10. deltaPi=Pi(sch)-Nii +EiEiGii 11. deltaQi=Qi(sch)+Hii +EiEiBii next i IMPROVED PSEUDO CODE This will require a single sweep on line data and a single loop on bus data. Gij,Bij can be stored along np[],nq[] loop k=1,lines i=np[k]; j=nq[k]; EiEj=E[i]*E[j]; Oij=O[i]-O[j]; Oji=-Oij; Hij=EiEj*(Gij*sin(Oij)-Bij*cos(Oij)); Not symmetrical Hji=EiEj*(Gij*sin(Oji)-Bij*cos(Oji)); Hii-=Hij; Hjj-=Hji; Nij=EiEj*(Gij*cos(Oij)+Bij*sin(Oij)); Nji=EiEj*(Gij*cos(Oji)+Bij*sin(Oji)); Nii+=Nij; Njj+=Nji; Jij=-Nij; Jji=-Nji; Lij=Hij; Lji=Hji; next k loop i=1,nbus EiEi=E[i]*E[i]; Jii=Nii; //Nii is not complete yet Nii+=2EiEi*Gii; Lii=-Hii-2EiEi*Bii; deltaP[i]=Pi(sch) -Nii +EiEi*Gii; deltaQ[i]=Qi(sch) +Hii +EiEi*Bii; next i

Program Listing_1: NEWTON RAPHSON LOADFLOW - with row/col delete to solve Ax=b

1: /* Program  NRLF1 -powerflow by Newton-Raphson
   2:    1. contiguous bus numbering
   3:    2. full matrix               
   4:    3. single sweep on line data to form the Jacobian
   5:    4. Slack, PV, randomly located
   6:    5. Ybus, linear storage
   7:    Date: October 15, 1999  Vancouver,BC, Canada
   8:    file: nrlfc1.cpp
   9:    Sparsity Technology, Joesel Pante
  10:    http://www.geocities.com/siliconvalley/lab/4223   for details
  11: */
  12: #define LMAX 20
  13: #define BMAX  15
  14: #include <stdio.h>
  15: #include <math.h>
  16: 
  17: int nl,nb,size,np[LMAX],nq[LMAX];  //p,q bus    
  18: float base=100;  //base Mva 
  19: float gii[BMAX],bii[BMAX],gij[LMAX],bij[LMAX];   //Ybus,linear storage
  20: int bus[BMAX],type[BMAX],tcol[2*BMAX],tnum[2*BMAX];     
  21: //type[] is used to guide lufact on Ax=b, as A contains Slack and PV bus.
  22: //without having to delete row-column.
  23: 
  24: float Pg[BMAX],Qg[BMAX],Pd[BMAX],Qd[BMAX],Qgmin[BMAX],Qgmax[BMAX]; 
  25: float yc[LMAX],Pij[LMAX],Qij[LMAX]; //line flow 
  26: float Es[BMAX],E[BMAX],d[BMAX];   //Es[]=voltage for PV    
  27: float A[20][21],x[2*BMAX];  //jacobian                        
  28: FILE *file;
  29: 
  30: void Enter(void)
  31: {char *ch; printf("Enter to continue\n");gets(ch);}
  32:           
  33: 
  34: void input(void)
  35: { int k,m,i;     
  36:   float r,x,ch,a1;     
  37:   nb=0,nl=0; //initialize counters nl=lines, nb=bus 
  38:   for (i=0;i<BMAX;i++) 
  39:    {gii[i]=0;       //initialize Yii
  40:     bii[i]=0; 
  41:    } 
  42:   printf(" count  p   q     r      x      MVAR\n");
  43:   while (1)
  44:    { nl++; 
  45:      fscanf(file,"%d %d %f %f %f",&k,&m,&r,&x,&ch);
  46:      if (k==99) break;
  47:      printf("%5d %3d %3d %7.2f %7.2f %7.2f\n",nl,k,m,r,x,ch);
  48:      if (k>nb) nb=k; if (m>nb) nb=m;  /*get the largest bus*/    
  49:      np[nl]=k; nq[nl]=m;
  50:      a1=r*r+x*x;
  51:      gij[nl]=-r/a1;   //off-diagonal
  52:      bij[nl]=x/a1; 
  53:      yc[nl]=ch/200;    //charging in admittance is +j
  54:      gii[k]-=gij[nl]; //diagonal
  55:      gii[m]-=gij[nl];   
  56:      bii[k]=bii[k]-bij[nl]+yc[nl];    //charging to pi and pu
  57:      bii[m]=bii[m]-bij[nl]+yc[nl];
  58:       
  59:     }
  60:    nl--;     //correct nl=lines read  
  61:    printf("number of lines:%2d  number of nodes:%2d\n",nl,nb);  
  62:    nb=0; 
  63:    printf("count Bus Type    E      Pg      Qg      Pd      Qd     Qgmin    Qgmax\n");
  64:    while (1)
  65:    { nb++; 
  66:      fscanf(file,"%d %d %f %f %f %f %f %f %f",
  67:      bus+nb,type+nb,E+nb,Pg+nb,Qg+nb,Pd+nb,Qd+nb,Qgmin+nb,Qgmax+nb); 
  68:      Es[nb]=E[nb];
  69:      if (bus[nb]==99) break;
  70:      printf("%5d %3d %3d %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f\n",
  71:      nb,bus[nb],type[nb],E[nb],Pg[nb],Qg[nb],Pd[nb],Qd[nb],Qgmin[nb],Qgmax[nb]); 
  72:      Pg[nb]/=base;Qg[nb]/=base;Pd[nb]/=base;Qd[nb]/=base;Qgmin[nb]/=base;Qgmax[nb]/=base;
  73:    }   
  74:    nb--;  
  75: }
  76:  
  77:  
  78: void printY(void)  
  79: {int j; 
  80:  printf("count np  nq  gij[]  bij[]\n");  
  81:  for (j=1;j<=nl;j++)
  82:   printf("%4d %4d %4d %9.4f %9.4f\n",j,np[j],nq[j],gij[j],bij[j]); 
  83:  printf("count gii[]  bii[]\n");  
  84:  for (j=1;j<=nb;j++)
  85:   printf("%4d  %9.4f %9.4f\n",j,gii[j],bii[j]); 
  86:   
  87:  Enter();
  88: }     
  89: 
  90: void printA(void)  
  91: {int i,j; 
  92:  printf("The Jacobian:\n");  
  93:  for (i=1;i<=size;i++) 
  94:    printf("%9d",tcol[i]);printf("\n"); 
  95:  for (i=1;i<=size;i++) 
  96:    printf("%9d",tnum[i]);printf("\n");      
  97:  for (i=1;i<=size;i++)
  98:    {printf("%2d ",tcol[i]); 
  99:     printf("%2d ",tnum[i]);  
 100:     for (j=1;j<=size;j++)
 101:      printf("%9.4f",A[i][j]); 
 102:     printf("%9.4f\n",x[i]);  
 103:    }   
 104:  printf("\n"); 
 105:  Enter();
 106: }               
 107:  
 108: //Solution of simultaneous equations Ax=b
 109: //by unsymmetric LU Factorization              
 110: void lufact(void)
 111: {int i,j,k;
 112:  for (k=1; k<=size-1; k++)
 113:   {for (i=k+1; i<=size; i++)  //solve for L     
 114:       A[i][k]/=A[k][k];
 115:    for (i=k+1; i<=size; i++)   //block  
 116:       for (j=k+1; j<=size; j++) 
 117:         A[i][j]=A[i][j]-A[k][j]*A[i][k];
 118:         
 119:    }
 120:   printf("After factoring\n");  
 121:   printA();
 122: //forward course
 123:  for (k=1;k<=size-1; k++) 
 124:    for (i=k+1; i<=size; i++)  
 125:       x[i]=x[i]-A[i][k]*x[k]; 
 126: //backward course
 127:  x[size]/=A[size][size];
 128:  for (k=size-1;k>=1; k--)   
 129:    {for (j=k+1;j<=size;j++)  
 130:       x[k]=x[k]-A[k][j]*x[j];   
 131:     x[k]/=A[k][k];  
 132:    } 
 133: //
 134:   printf("After forward and backward course\n");                      
 135:   printA();
 136: }             
 137:   
 138: void delRC(int rc)
 139: {int ii,j1,i1,jj;        
 140: //delete column
 141:  for (i1=1;i1<=size;i1++) 
 142:    for (j1=rc+1;j1<=size+1;j1++)  
 143:      {ii=i1;jj=j1-1;
 144:       A[ii][jj]=A[i1][j1];
 145:      }   
 146:      
 147: // printf("After delete column\n");  printA();        
 148: //delete row  
 149:  for (i1=rc+1;i1<=size;i1++) 
 150:    for (j1=1;j1<=size;j1++)  
 151:      {ii=i1-1;jj=j1;
 152:       A[ii][jj]=A[i1][j1];
 153:      }                            
 154: //delete pointer too
 155:   for (i1=rc+1;i1<=size;i1++)  
 156:    {ii=i1-1;
 157:     tcol[ii]=tcol[i1]; tnum[ii]=tnum[i1];   
 158:     x[ii]=x[i1];                         //delete deltaPQ!
 159:    }
 160:   size--;   
 161:  //printf("After delete row\n");  printA();                
 162: 
 163: }
 164:   
 165: void deleRC(void)
 166: {int i=0; 
 167:  while (i<=size)
 168:    {i++;
 169:     if (tcol[i]) 
 170:       {delRC(i);
 171:        i--;     //move back pointer
 172:       }
 173:    }   
 174: }
 175: 
 176: void nrlfc1(void)
 177: {int iteration=0,itermax=20,i,j,k;
 178:  float maxPQ=0.01,delPQ;
 179:  float sinOij,cosOij,EiEj;                     
 180: //form the Jacobian matrix by single sweep on the line data  
 181:  do
 182:    {delPQ=0;   //set max error
 183:     iteration++; 
 184:      size=nb*2;  
 185: //replicate type[1->nb] to type[nb+1->2*nb]  
 186:   for (i=1;i<=nb;i++)
 187:    {tcol[i]=type[i];tcol[i+nb]=type[i];    
 188:     tnum[i]=i;tnum[i+nb]=i;
 189:    }   
 190: //zero out the square Jacobian matrix first 
 191:     for (i=1; i<=size; i++)
 192:       for (j=1; j<=size;j++)
 193:          A[i][j]=0;
 194:     for (k=1;k<=nl;k++)     //single sweep on line data to form HJNL
 195:      {i=np[k];j=nq[k];
 196: //Hij=EiEj(Gij.sinOij -Bij.cosOij)
 197:       sinOij=(float)sin(d[i]-d[j]);
 198:       cosOij=(float)cos(d[i]-d[j]); 
 199:       EiEj=E[i]*E[j];
 200:       A[i][j]=EiEj*(gij[k]*sinOij-bij[k]*cosOij);    //Hij   
 201: //Nij=EiEj(Gij.cosOij +Bij.sinOij)   
 202:       A[i][j+nb]=EiEj*(gij[k]*cosOij+bij[k]*sinOij); //Nij  
 203:       sinOij=-sinOij;
 204:       A[j][i]=EiEj*(gij[k]*sinOij-bij[k]*cosOij);    //Hji
 205:       A[j][i+nb]=EiEj*(gij[k]*cosOij+bij[k]*sinOij); //Nji
 206: //Jij= -Nij
 207:       A[i+nb][j]= -A[i][j+nb];
 208:       A[j+nb][i]= -A[j][i+nb];
 209: //Lij=Hij
 210:       A[i+nb][j+nb]= A[i][j];
 211:       A[j+nb][i+nb]= A[j][i];
 212: //Hii= sum -Hij
 213:       A[i][i]-=A[i][j];
 214:       A[j][j]-=A[j][i];
 215: //Nii= sum Nij
 216:       A[i][i+nb]+=A[i][j+nb];
 217:       A[j][j+nb]+=A[j][i+nb];
 218:      }
 219:   
 220: //include j=i. Loop on bus data
 221:     for (i=1;i<=nb;i++)
 222:      {A[i+nb][i]=A[i][i+nb];           //Jii=Nii  
 223:       A[i][i+nb]+=2*E[i]*E[i]*gii[i];  //Nii= sum Nij +2EiEiGii  
 224:       A[i+nb][i+nb]= -A[i][i] -2*E[i]*E[i]*bii[i]; //Lii=-Hii-2EiEiBii   
 225: //delta P=Pi(sch) - Nii + EiEiGii    
 226:       x[i]=Pg[i]-Pd[i] -A[i][i+nb] + E[i]*E[i]*gii[i];  
 227:       if (tcol[i]!=2)
 228:         if((float)fabs(x[i])>delPQ) delPQ=(float)fabs(x[i]);    
 229: //delta Q=Qi(sch) + Hii + EiEiBii  
 230:       x[i+nb]=Qg[i]-Qd[i] +A[i][i] +E[i]*E[i]*bii[i]; 
 231:       if (!tcol[i+nb])
 232:          if((float)fabs(x[i+nb])>delPQ) delPQ=(float)fabs(x[i+nb]);    
 233:      }            
 234:      printf("iter:%2d max delPQ:%8.6f\n",iteration,delPQ);
 235:      printA(); 
 236:      deleRC();
 237:      lufact(); 
 238: //update voltage angle and magnitude      
 239:      k=size/2;  
 240:      for (i=1;i<=k;i++)
 241:        {j=tnum[i];
 242:         d[j]+=x[i];
 243:         E[j]=E[j]*(1+x[i+k]);
 244:         printf("E[%2d]=%9.6f d=%9.6f \n",j,E[j],d[j]);
 245:        }                                             
 246:   }while (iteration<=itermax && delPQ>maxPQ); 
 247: }                    
 248: 
 249: void pflows(void)
 250: {int i,k,m;    
 251:  float angle,co,si,Gij,Bij;
 252:  for (i=1; i<=nl; i++) 
 253:   {k=np[i];m=nq[i]; 
 254:    angle=d[k]-d[m]; 
 255:    co=(float)cos(angle);
 256:    si=(float)sin(angle); 
 257:    Gij=-gij[i];Bij=-bij[i];  //from Ybus
 258:    Pij[i]=Gij*E[k]*E[k]-E[k]*E[m]*(Gij*co+Bij*si);
 259:    Qij[i]=-(yc[i]+Bij)*E[k]*E[k]-E[k]*E[m]*(Gij*si-Bij*co); 
 260:    printf("i:%2d j:%2d Pij:%7.2f Qij:%7.2f\n",k,m,Pij[i]*base,Qij[i]*base);
 261:   }
 262:  }     
 263:   
 264: //Ybus is not a square matrix. Doesn`t negate Gij,Bij
 265: void gener(void)  
 266: {int i,j,m;     
 267:  float co,si;
 268:  for (i=1; i<=nb; i++)
 269:   if (type[i]==2)
 270:     {Pg[i]=Qg[i]=0; //solve Pg,Qg   
 271: //line data, off diagonals only    
 272:      for (m=1;m<=nl; m++)             
 273:        {if (np[m]==i)
 274:           {j=nq[m];
 275:            co=(float)cos(d[i]-d[j]);
 276:            si=(float)sin(d[i]-d[j]);  
 277:            Pg[i]=Pg[i]+E[j]*(gij[m]*co+bij[m]*si);
 278:            Qg[i]=Qg[i]+E[j]*(gij[m]*si-bij[m]*co); 
 279:           }  
 280:         else if (nq[m]==i)
 281:           {j=np[m];
 282:            co=(float)cos(d[i]-d[j]);
 283:            si=(float)sin(d[i]-d[j]);  
 284:            Pg[i]=Pg[i]+E[j]*(gij[m]*co+bij[m]*si);
 285:            Qg[i]=Qg[i]+E[j]*(gij[m]*si-bij[m]*co); 
 286:           } 
 287:        }                                 
 288: //include diagonal, j=i       
 289:      Pg[i]+=E[i]*gii[i];
 290:      Qg[i]-=E[i]*bii[i]; 
 291: //Ei before the summation formula       
 292:      Pg[i]*=E[i];
 293:      Qg[i]*=E[i];  
 294: //add load
 295:      Pg[i]+=Pd[i];
 296:      Qg[i]+=Qd[i];     
 297:      printf("Slack Pg[%2d]=%7.2f Qg=%7.2f\n",i,Pg[i]*base,Qg[i]*base);  
 298:     }
 299:    else if (type[i]==1) 
 300:     {Qg[i]=0;  //solve Qg
 301:      for (m=1;m<=nl; m++)          //includes j=i
 302:         {if (np[m]==i)
 303:           {j=nq[m];
 304:            co=(float)cos(d[i]-d[j]);
 305:            si=(float)sin(d[i]-d[j]); 
 306:            Qg[i]=Qg[i]+E[j]*(gij[m]*si-bij[m]*co); 
 307:           }   
 308:          else if (nq[m]==i)
 309:           {j=np[m];
 310:            co=(float)cos(d[i]-d[j]);
 311:            si=(float)sin(d[i]-d[j]); 
 312:            Qg[i]=Qg[i]+E[j]*(gij[m]*si-bij[m]*co); 
 313:           }  
 314:         }  
 315:  //include diagonal, j=i       
 316:      Qg[i]-=E[i]*bii[i]; 
 317: //Ei before the summation formula         
 318:      Qg[i]*=E[i]; 
 319:      Qg[i]+=Qd[i];
 320:      printf("PV Qg[%2d]=%7.2f\n",i,Qg[i]*base);       
 321:     }
 322: }
 323:         
 324: void main(void)
 325: {int i;
 326:  file=fopen("nrlf.dat","r");
 327:  input();            //read line,bus data
 328:  printY();           //print Ybus=g[][]+jb[][];
 329:  nrlfc1();           //iterate voltage
 330:  pflows();           //solve power flows
 331:  gener();            //solve generation, Pg,Qg=SL, Qg=PV
 332:  i=fclose(file);
 333: }

OUTPUT:
 count  p   q     r      x      MVAR                                   
    1   1   2    0.08    0.24    5.00                                  
    2   1   3    0.02    0.06    6.00                                  
    3   2   3    0.06    0.18    4.00                                  
number of lines: 3  number of nodes: 3                                 
count Bus Type    E      Pg      Qg      Pd      Qd     Qgmin    Qgmax 
    1   1   2    1.05    0.00    0.00    0.00    0.00    0.00    0.00  
    2   2   0    1.00   20.00   18.93   50.00   20.00    0.00    0.00  
    3   3   0    1.00    0.00    0.00   60.00   25.00    0.00    0.00  
count np  nq  gij[]  bij[]                                             
   1    1    2   -1.2500    3.7500                                     
   2    1    3   -5.0000   15.0000                                     
   3    2    3   -1.6667    5.0000                                     
count gii[]  bii[]                                                     
   1     6.2500  -18.6950                                              
   2     2.9167   -8.7050                                              
   3     6.6667  -19.9500                                              
Enter to continue
The Jacobian:                                                        
        2        0        0        2        0        0               
        1        2        3        1        2        3               
 2  1   19.6875  -3.9375 -15.7500   7.2187  -1.3125  -5.2500  -0.3281
 0  2   -3.9375   8.9375  -5.0000  -1.3125   2.8542  -1.6667  -0.2375
 0  3  -15.7500  -5.0000  20.7500  -5.2500  -1.6667   6.4167  -0.3500
 2  1   -6.5625   1.3125   5.2500  21.5350  -3.9375 -15.7500  -0.9237
 0  2    1.3125  -2.9792   1.6667  -3.9375   8.4725  -5.0000   0.2218
 0  3    5.2500   1.6667  -6.9167 -15.7500  -5.0000  19.1500   0.5500
                                                                    
Enter to continue                                                    
                                                                     
After factoring                                                      
The Jacobian:                                                        
        0        0        0        0                                 
        2        3        2        3                                 
 0  2    8.9375  -5.0000   2.8542  -1.6667  -0.2375                  
 0  3   -0.5594  17.9528  -0.0699   5.4843  -0.3500                  
 0  2   -0.3333  -0.0000   9.4239  -5.5556   0.2218                  
 0  3    0.1865  -0.3333  -0.5895  18.0138   0.5500                  
                                                                    
Enter to continue                                                    
After forward and backward course                                     
The Jacobian:                                                         
        0        0        0        0                                  
        2        3        2        3                                  
 0  2    8.9375  -5.0000   2.8542  -1.6667  -0.0513                   
 0  3   -0.5594  17.9528  -0.0699   5.4843  -0.0355                   
 0  2   -0.3333  -0.0000   9.4239  -5.5556   0.0321                   
 0  3    0.1865  -0.3333  -0.5895  18.0138   0.0287                   
                                                                     
Enter to continue                                                     
                                                                      
E[ 2]= 1.032068 d=-0.051344                                           
E[ 3]= 1.028723 d=-0.035546
--------------------------------------                                        
iter: 2 max delPQ:0.026448                                            
The Jacobian:                                                         
        2        0        0        2        0        0                
        1        2        3        1        2        3                
 2  1   20.5120  -4.1279 -16.3841   7.8154  -1.1442  -4.8216  -0.9248 
 0  2   -3.9889   9.2688  -5.2799  -1.5614   2.7989  -1.8532   0.0078 
 0  3  -16.0002  -5.3359  21.3361  -5.9732  -1.6854   6.4517   0.0035 
 2  1   -5.9658   1.1442   4.8216  20.7104  -4.1279 -16.3841  -0.0992 
 0  2    1.5614  -3.4145   1.8532  -3.9889   9.2757  -5.2799  -0.0141 
 0  3    5.9732   1.6854  -7.6586 -16.0002  -5.3359  20.8890  -0.0264 
                                                                     
Enter to continue
After factoring                                                        
The Jacobian:                                                          
        0        0        0        0                                   
        2        3        2        3                                   
 0  2    9.2688  -5.2799   2.7989  -1.8532   0.0078                    
 0  3   -0.5757  18.2965  -0.0742   5.3848   0.0035                    
 0  2   -0.3684  -0.0050  10.3064  -5.9356  -0.0141                    
 0  3    0.1818  -0.3661  -0.5697  19.8157  -0.0264                    
                                                                      
Enter to continue                                                      
                                                                       
After forward and backward course                                      
The Jacobian:                                                          
        0        0        0        0                                   
        2        3        2        3                                   
 0  2    9.2688  -5.2799   2.7989  -1.8532   0.0016                    
 0  3   -0.5757  18.2965  -0.0742   5.3848   0.0009                    
 0  2   -0.3684  -0.0050  10.3064  -5.9356  -0.0020                    
 0  3    0.1818  -0.3661  -0.5697  19.8157  -0.0016                    
                                                                      
Enter to continue                                                      
                                                                       
E[ 2]= 1.030006 d=-0.049707                                            
E[ 3]= 1.027097 d=-0.034653                                            
iter: 3 max delPQ:0.000048                                             
                                  
Enter to continue                                                      
E[ 2]= 1.030001 d=-0.049704          
E[ 3]= 1.027094 d=-0.034652          
i: 1 j: 2 Pij:  22.94 Qij:  -1.10    
i: 1 j: 3 Pij:  68.39 Qij:  15.06    
i: 2 j: 3 Pij:  -7.44 Qij:   2.09    
Slack Pg[ 1]=  91.33 Qg=  13.96

Code Explanation:
There is no need to store the Ybus matrix as a square. The formation of HJNL uses a single sweep on the line data and the Ybus elements are accessed in that order. The diagonal is stored with the bus data. The Jacobian matrix is coded with bus type in vector tcol[2*nb] and bus number in vector tnum[2*nb]. Remember that the matrix size is twice the bus data. Column and row entries for Slack and PV bus are deleted to simulate how to do it using a calculator. The bus number code tnum[] is also deleted for the entry of Slack and PV type. Update of voltage angle and magnitude are only for those nodes in tnum[].

Program Listing_2: NEWTON RAPHSON LOADFLOW - no delete to solve Ax=b

1: /* Program  NRLF2 -powerflow by Newton-Raphson
   2:    1. contiguous bus numbering
   3:    2. full matrix               
   4:    3. single sweep on line data to form the Jacobian
   5:    4. Slack, PV, randomly located
   6:    5. Ybus, linear storage   
   7:    6. No delete of row/col on the Jacobian to perform Ax=b
   8:    Date: October 25, 1999  Vancouver,BC, Canada
   9:    file: nrlfc1.cpp
  10:    Sparsity Technology, Joesel Pante
  11:    http://www.geocities.com/siliconvalley/lab/4223   for details
  12: */
  13: #define LMAX 20
  14: #define BMAX  15
  15: #include <stdio.h>
  16: #include <math.h>
  17: 
  18: int nl,nb,size,np[LMAX],nq[LMAX];  //p,q bus    
  19: float base=100;  //base Mva 
  20: float gii[BMAX],bii[BMAX],gij[LMAX],bij[LMAX];   //Ybus,linear storage
  21: int bus[BMAX],type[BMAX],tcol[2*BMAX],tnum[2*BMAX];     
  22: //type[] is used to guide lufact on Ax=b, as A contains Slack and PV bus.
  23: //without having to delete row-column.
  24: 
  25: float Pg[BMAX],Qg[BMAX],Pd[BMAX],Qd[BMAX],Qgmin[BMAX],Qgmax[BMAX]; 
  26: float yc[LMAX],Pij[LMAX],Qij[LMAX]; //line flow 
  27: float Es[BMAX],E[BMAX],d[BMAX];   //Es[]=voltage for PV    
  28: float A[20][21],x[2*BMAX];  //jacobian                        
  29: FILE *file;
  30: 
  31: void Enter(void)
  32: {char *ch; printf("Enter to continue\n");gets(ch);}
  33:           
  34: 
  35: void input(void)
  36: { int k,m,i;     
  37:   float r,x,ch,a1;     
  38:   nb=0,nl=0; //initialize counters nl=lines, nb=bus 
  39:   for (i=0;i<BMAX;i++) 
  40:    {gii[i]=0;       //initialize Yii
  41:     bii[i]=0; 
  42:    } 
  43:   printf(" count  p   q     r      x      MVAR\n");
  44:   while (1)
  45:    { nl++; 
  46:      fscanf(file,"%d %d %f %f %f",&k,&m,&r,&x,&ch);
  47:      if (k==99) break;
  48:      printf("%5d %3d %3d %7.2f %7.2f %7.2f\n",nl,k,m,r,x,ch);
  49:      if (k>nb) nb=k; if (m>nb) nb=m;  /*get the largest bus*/    
  50:      np[nl]=k; nq[nl]=m;
  51:      a1=r*r+x*x;
  52:      gij[nl]=-r/a1;   //off-diagonal
  53:      bij[nl]=x/a1; 
  54:      yc[nl]=ch/200;    //charging in admittance is +j
  55:      gii[k]-=gij[nl]; //diagonal
  56:      gii[m]-=gij[nl];   
  57:      bii[k]=bii[k]-bij[nl]+yc[nl];    //charging to pi and pu
  58:      bii[m]=bii[m]-bij[nl]+yc[nl];
  59:       
  60:     }
  61:    nl--;     //correct nl=lines read  
  62:    printf("number of lines:%2d  number of nodes:%2d\n",nl,nb);  
  63:    nb=0; 
  64:    printf("count Bus Type    E      Pg      Qg      Pd      Qd     Qgmin    Qgmax\n");
  65:    while (1)
  66:    { nb++; 
  67:      fscanf(file,"%d %d %f %f %f %f %f %f %f",
  68:      bus+nb,type+nb,E+nb,Pg+nb,Qg+nb,Pd+nb,Qd+nb,Qgmin+nb,Qgmax+nb); 
  69:      Es[nb]=E[nb];
  70:      if (bus[nb]==99) break;
  71:      printf("%5d %3d %3d %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f\n",
  72:      nb,bus[nb],type[nb],E[nb],Pg[nb],Qg[nb],Pd[nb],Qd[nb],Qgmin[nb],Qgmax[nb]); 
  73:      Pg[nb]/=base;Qg[nb]/=base;Pd[nb]/=base;Qd[nb]/=base; Qgmin[nb]/=base;Qgmax[nb]/=base;
  74:    }   
  75:    nb--;  
  76: }
  77:  
  78:  
  79: void printY(void)  
  80: {int j; 
  81:  printf("count np  nq  gij[]  bij[]\n");  
  82:  for (j=1;j<=nl;j++)
  83:   printf("%4d %4d %4d %9.4f %9.4f\n",j,np[j],nq[j],gij[j],bij[j]); 
  84:  printf("count gii[]  bii[]\n");  
  85:  for (j=1;j<=nb;j++)
  86:   printf("%4d  %9.4f %9.4f\n",j,gii[j],bii[j]); 
  87:   
  88:  Enter();
  89: }     
  90: 
  91: void printA(void)  
  92: {int i,j; 
  93:  printf("The Jacobian:\n");  
  94:  for (i=1;i<=size;i++) 
  95:    printf("%9d",tcol[i]);printf("\n"); 
  96:  for (i=1;i<=size;i++) 
  97:    printf("%9d",tnum[i]);printf("\n");      
  98:  for (i=1;i<=size;i++)
  99:    {printf("%2d ",tcol[i]); 
 100:     printf("%2d ",tnum[i]);  
 101:     for (j=1;j<=size;j++)
 102:      printf("%9.4f",A[i][j]); 
 103:     printf("%9.4f\n",x[i]);  
 104:    }   
 105:  printf("\n"); 
 106:  Enter();
 107: }               
 108:  
 109: //Solution of simultaneous equations Ax=b
 110: //by unsymmetric LU Factorization              
 111: void lufact(void)
 112: {int i,j,k;
 113:  for (k=1; k<size; k++)  
 114:    if (tcol[k]!=2)
 115:      {for (i=k+1; i<=size; i++)  //solve for L     
 116:         A[i][k]/=A[k][k];
 117:       for (i=k+1; i<=size; i++)   //block    
 118:         if (tcol[i]!=2)
 119:           for (j=k+1; j<=size; j++)    
 120:              if (tcol[k]!=2)
 121:                 A[i][j]=A[i][j]-A[k][j]*A[i][k];
 122:         
 123:      }
 124:   printf("After factoring\n");  
 125:   printA();
 126: //forward course
 127:  for (k=1;k<size; k++) 
 128:    if (tcol[k]!=2)
 129:     for (i=k+1; i<=size; i++)      //does not contain tcol[]
 130:       x[i]=x[i]-A[i][k]*x[k]; 
 131: //backward course
 132:  x[size]/=A[size][size];
 133:  for (k=size-1;k>=1; k--)   
 134:    {if (tcol[k]!=2)
 135:     for (j=k+1;j<=size;j++)     
 136:       if (tcol[j]!=2)
 137:         x[k]=x[k]-A[k][j]*x[j];   
 138:     x[k]/=A[k][k];         
 139:    }                    
 140: //
 141:   printf("After forward and backward course\n");                      
 142:   printA();
 143: }             
 144:   
 145: void nrlfc1(void)
 146: {int iteration=0,itermax=20,i,j,k;
 147:  float maxPQ=0.01,delPQ;
 148:  float sinOij,cosOij,EiEj;                     
 149: //form the Jacobian matrix by single sweep on the line data 
 150:  size=nb*2;  //outside the loop, set only once  
 151: //replicate type[1->nb] to type[nb+1->2*nb]  
 152:   for (i=1;i<=nb;i++)
 153:    {tcol[i]=type[i];tcol[i+nb]=type[i];    
 154:     tnum[i]=i;tnum[i+nb]=i;
 155:    }   
 156:  do
 157:    {delPQ=0;   //set max error
 158:     iteration++; 
 159: //zero out the square Jacobian matrix first 
 160:     for (i=1; i<=size; i++)
 161:       for (j=1; j<=size;j++)
 162:          A[i][j]=0;
 163:     for (k=1;k<=nl;k++)     //single sweep on line data to form HJNL
 164:      {i=np[k];j=nq[k];
 165: //Hij=EiEj(Gij.sinOij -Bij.cosOij)
 166:       sinOij=(float)sin(d[i]-d[j]);
 167:       cosOij=(float)cos(d[i]-d[j]); 
 168:       EiEj=E[i]*E[j];
 169:       A[i][j]=EiEj*(gij[k]*sinOij-bij[k]*cosOij);    //Hij   
 170: //Nij=EiEj(Gij.cosOij +Bij.sinOij)   
 171:       A[i][j+nb]=EiEj*(gij[k]*cosOij+bij[k]*sinOij); //Nij  
 172:       sinOij=-sinOij;
 173:       A[j][i]=EiEj*(gij[k]*sinOij-bij[k]*cosOij);    //Hji
 174:       A[j][i+nb]=EiEj*(gij[k]*cosOij+bij[k]*sinOij); //Nji
 175: //Jij= -Nij
 176:       A[i+nb][j]= -A[i][j+nb];
 177:       A[j+nb][i]= -A[j][i+nb];
 178: //Lij=Hij
 179:       A[i+nb][j+nb]= A[i][j];
 180:       A[j+nb][i+nb]= A[j][i];
 181: //Hii= sum -Hij
 182:       A[i][i]-=A[i][j];
 183:       A[j][j]-=A[j][i];
 184: //Nii= sum Nij
 185:       A[i][i+nb]+=A[i][j+nb];
 186:       A[j][j+nb]+=A[j][i+nb];
 187:      }
 188:   
 189: //include j=i. Loop on bus data
 190:     for (i=1;i<=nb;i++)
 191:      {A[i+nb][i]=A[i][i+nb];           //Jii=Nii  
 192:       A[i][i+nb]+=2*E[i]*E[i]*gii[i];  //Nii= sum Nij +2EiEiGii  
 193:       A[i+nb][i+nb]= -A[i][i] -2*E[i]*E[i]*bii[i]; //Lii=-Hii-2EiEiBii   
 194: //delta P=Pi(sch) - Nii + EiEiGii    
 195:       x[i]=Pg[i]-Pd[i] -A[i][i+nb] + E[i]*E[i]*gii[i];  
 196:       if (tcol[i]!=2)
 197:         if((float)fabs(x[i])>delPQ) delPQ=(float)fabs(x[i]);    
 198: //delta Q=Qi(sch) + Hii + EiEiBii  
 199:       x[i+nb]=Qg[i]-Qd[i] +A[i][i] +E[i]*E[i]*bii[i]; 
 200:       if (!tcol[i+nb])
 201:          if((float)fabs(x[i+nb])>delPQ) delPQ=(float)fabs(x[i+nb]);    
 202:      }            
 203:      printf("iter:%2d max delPQ:%8.6f\n",iteration,delPQ);
 204:      printA(); 
 205:      lufact(); 
 206: //update voltage angle and magnitude      
 207:      for (i=1;i<=nb;i++) 
 208:        if (tcol[i]!=2)
 209:        {j=tnum[i];
 210:         d[j]+=x[i];
 211:         E[j]=E[j]*(1+x[i+nb]);
 212:         printf("E[%2d]=%9.6f d=%9.6f \n",j,E[j],d[j]);
 213:        }                                             
 214:   }while (iteration<=itermax && delPQ>maxPQ); 
 215: }                    
 216: 
 217: void pflows(void)
 218: {int i,k,m;    
 219:  float angle,co,si,Gij,Bij;
 220:  for (i=1; i<=nl; i++) 
 221:   {k=np[i];m=nq[i]; 
 222:    angle=d[k]-d[m]; 
 223:    co=(float)cos(angle);
 224:    si=(float)sin(angle); 
 225:    Gij=-gij[i];Bij=-bij[i];  //from Ybus
 226:    Pij[i]=Gij*E[k]*E[k]-E[k]*E[m]*(Gij*co+Bij*si);
 227:    Qij[i]=-(yc[i]+Bij)*E[k]*E[k]-E[k]*E[m]*(Gij*si-Bij*co); 
 228:    printf("i:%2d j:%2d Pij:%7.2f Qij:%7.2f\n",k,m,Pij[i]*base,Qij[i]*base);
 229:   }
 230:  }     
 231:   
 232: //Ybus is not a square matrix. Doesn`t negate Gij,Bij
 233: void gener(void)  
 234: {int i,j,m;     
 235:  float co,si;
 236:  for (i=1; i<=nb; i++)
 237:   if (type[i]==2)
 238:     {Pg[i]=Qg[i]=0; //solve Pg,Qg   
 239: //line data, off diagonals only    
 240:      for (m=1;m<=nl; m++)             
 241:        {if (np[m]==i)
 242:           {j=nq[m];
 243:            co=(float)cos(d[i]-d[j]);
 244:            si=(float)sin(d[i]-d[j]);  
 245:            Pg[i]=Pg[i]+E[j]*(gij[m]*co+bij[m]*si);
 246:            Qg[i]=Qg[i]+E[j]*(gij[m]*si-bij[m]*co); 
 247:           }  
 248:         else if (nq[m]==i)
 249:           {j=np[m];
 250:            co=(float)cos(d[i]-d[j]);
 251:            si=(float)sin(d[i]-d[j]);  
 252:            Pg[i]=Pg[i]+E[j]*(gij[m]*co+bij[m]*si);
 253:            Qg[i]=Qg[i]+E[j]*(gij[m]*si-bij[m]*co); 
 254:           } 
 255:        }                                 
 256: //include diagonal, j=i       
 257:      Pg[i]+=E[i]*gii[i];
 258:      Qg[i]-=E[i]*bii[i]; 
 259: //Ei before the summation formula       
 260:      Pg[i]*=E[i];
 261:      Qg[i]*=E[i];  
 262: //add load
 263:      Pg[i]+=Pd[i];
 264:      Qg[i]+=Qd[i];     
 265:      printf("Slack Pg[%2d]=%7.2f Qg=%7.2f\n",i,Pg[i]*base,Qg[i]*base);  
 266:     }
 267:    else if (type[i]==1) 
 268:     {Qg[i]=0;  //solve Qg
 269:      for (m=1;m<=nl; m++)          //includes j=i
 270:         {if (np[m]==i)
 271:           {j=nq[m];
 272:            co=(float)cos(d[i]-d[j]);
 273:            si=(float)sin(d[i]-d[j]); 
 274:            Qg[i]=Qg[i]+E[j]*(gij[m]*si-bij[m]*co); 
 275:           }   
 276:          else if (nq[m]==i)
 277:           {j=np[m];
 278:            co=(float)cos(d[i]-d[j]);
 279:            si=(float)sin(d[i]-d[j]); 
 280:            Qg[i]=Qg[i]+E[j]*(gij[m]*si-bij[m]*co); 
 281:           }  
 282:         }  
 283:  //include diagonal, j=i       
 284:      Qg[i]-=E[i]*bii[i]; 
 285: //Ei before the summation formula         
 286:      Qg[i]*=E[i]; 
 287:      Qg[i]+=Qd[i];
 288:      printf("PV Qg[%2d]=%7.2f\n",i,Qg[i]*base);       
 289:     }
 290: }
 291:         
 292: void main(void)
 293: {int i;
 294:  file=fopen("nrlf.dat","r");
 295:  input();            //read line,bus data
 296:  printY();           //print Ybus=g[][]+jb[][];
 297:  nrlfc1();           //iterate voltage
 298:  pflows();           //solve power flows
 299:  gener();            //solve generation, Pg,Qg=SL, Qg=PV
 300:  i=fclose(file);
 301: }

Code Explanation: The solution of Ax=b does not actually delete row/column due to slack bus. In order to visualize the process of factorization, we have included Listing-1. Shown below will be the modification in Listing-1 and -2 if we have to include PV busses. In Listing-2, just ignore the row and column after nb (nb=number of nodes) for PV busses. PV bus does not get a voltage update.

//For Listing-1 modification
void deleRC(void)
{int i=0,j=0; 
 while (i<=size)
   {i++;j++;        //j keeps track of nb
    if (tcol[i]==2 || tcol[i]==1 && j>nb) 
      {delRC(i);
       i--;     //move back pointer
      }
   }   
}
//---------------------
//update voltage angle from 1 to nb-1, exclude Slack     
     for (i=1;i<nb;i++)
       {j=tnum[i];
        d[j]+=x[i]; 
       }
//update voltage magnitude, from nb to size, excludes Slack and PV
     for (i=nb;i<=size;i++)        
       {j=tnum[i];
        E[j]=E[j]*(1+x[i]); 
       }
     for (i=1;i<=nb;i++)   
        printf("E[%2d]=%9.6f d=%9.6f \n",i,E[i],d[i]);
//For Listing-2 modification-------------------------
//Solution of simultaneous equations Ax=b
//by unsymmetric LU Factorization              
void lufact(void)
{int i,j,k;
 for (k=1; k<size; k++)  
   if (tcol[k]!=2)                    //exclude Slack
       {if (k<nb ||!tcol[k])          //exclude PV
          {for (i=k+1; i<=size; i++)  //solve for L     
              A[i][k]/=A[k][k];  
           for (i=k+1; i<=size; i++)   //block    
              for (j=k+1; j<=size; j++)    
                A[i][j]=A[i][j]-A[k][j]*A[i][k];
          }      
       } 
     
  printf("After factoring\n");  
  printA();
//forward course
 for (k=1;k<size; k++) 
   if (tcol[k]!=2)            //excludeSlack
     {if (k<nb || !tcol[k])    //exclude PV
       for (i=k+1; i<=size; i++)      //does not contain tcol[]
         x[i]=x[i]-A[i][k]*x[k]; 
     }    
//backward course
 x[size]/=A[size][size];
 for (k=size-1;k>=1; k--)   
   if (tcol[k]!=2)                   //exclude Slack
     {if (k<nb || !tcol[k])          //exclude PV
        for (j=k+1;j<=size;j++)     
          if (tcol[j]!=2)            //exclude Slack
           {if (j<nb || !tcol[j])    //exclude PV
               x[k]=x[k]-A[k][j]*x[j];
           }    
      x[k]/=A[k][k];         
     }                    
//
  printf("After forward and backward course\n");                      
  printA();
}             
///--------on the voltage update portion
//update voltage angle and magnitude      
     for (i=1;i<=nb;i++) 
       if (tcol[i]!=2)
       {j=tnum[i];
        d[j]+=x[i]; 
        if (!tcol[i+nb])                //exclude PV voltage update
          E[j]=E[j]*(1+x[i+nb]);
        printf("E[%2d]=%9.6f d=%9.6f \n",j,E[j],d[j]);
       }

Newton Raphson Using the Algorithm by Van Ness

The Jacobian matrix is in rectangular and updating voltage and angle is in polar. Let
NRLF_Eq. (16)
The equation above shows the polar form of the Ybus. This is needed in the derivation of formula which is tricky.
NRLF_Eq. (17)
Rewriting GSLF_Eq (5) and taking out the term when i=j is shown above. The summation simply means that we are substituting the non-zero elements of row=i from the Ybus. The first summation does not contain (i not equal to j), so it includes the diagonal of the Ybus. The next 2 summations excludes when i=j. The dot above the variable refers that it is a complex number. The asterisk means that we take the conjugate of the complex number. We can now separate the real part and equate it to the real part of the right hand side term. Yii is composed of Gii and jBii.
NRLF_Eq. (18)

Transferring the negative sign,
NRLF_Eq. (19)
From Eqs (18) and (19) we can simplify the derivation process, let A and B be substituted to as shown below,
NRLF_Eq. (20)
The purpose of the above Eqs (18) and (19) is to equate Pi and Qi in terms of and be able to get the derivatives of Pi and Qi with respect to and conform with the definition of the Jacobian matrix. The tricky part is that we will not use the sine and cosine terms. Watch on.

Now let us get the rectangular parts of

NRLF_Eq. (21)
NRLF_Eq. (22)
For simplification, aij, bij is substituted,
NRLF_Eq. (23)
We can now equate real part to real part, and imaginary part to imaginary part.
NRLF_Eq. (24)
NRLF_Eq. (25)
Eqs (18) and (19) can now be replaced with these new derived values. We will need Eq (26) in the programming code.
NRLF_Eq. (26)
The method of Van Ness is different from the ordinary Newton Raphson, because he group individual H,N,J,L for branch i-j near each other instead of 4 separate sub-matrices.
NRLF_Eq. (27)
NRLF_Eq. (31)
We now get the derivatives. I hope you know how to. The aspect of this website is to show the reader how to create the industry grade programming aspect, rather than the details of the mathematics involved.
NRLF_Eq. (32)
NRLF_Eq. (33)
NRLF_Eq. (34)
NRLF_Eq. (35)

Eqs (24) and (25) will be used instead of the sine and cosine terms to form the H,N,J,L.

We are now ready to show the programming code.

COMPUTER IMPLEMENTATION loop k=1, Lines i=pbus[k]; j=qbus[k]; 1. aij=e[j]Gij - f[j]Bij; aji=e[i]Gij - f[i]Bij; 2. bij=e[j]Bij + f[j]Gij; bji=e[i]Bij + f[i]Gij; 3. ReEYEij=e[i]aij + f[i]bij; ReEYEji=e[j]aji + f[j]bji; 4. ImEYEij=e[i]bij - f[i]aij; ImEYEji=e[j]bji - f[j]aji; 5. Hij= -ImEYEij; Hji= -ImEYEji; 6. Nij= ReEYEij; Nji= ReEYEji; 7. Jij= -ReEYEij; Jji= -ReEYEji; 8. Lij= -ImEYEij; Hji= -ImEYEji; collect Pi, Qi 9. Pi = Pi + ReEYEij; Pj = Pj + ReEYEji; 10.Qi = Qi + ImEYEij; Qj = Qj + ImEYEji; next k loop i=1, Nodes 11. Hii=Qi //Qi is not complete yet 12. Jii=Pi //Pi is not complete yet E[i] squared = e[i] squared + f[i] squared Let E2=e[i]e[i]+ f[i]f[i] 13. Pi= Pi + (E2)Gii 14. Qi= -Qi - (E2)Bii 15. Nii= Pi + (E2)Gii 16. Lii= Qi - (E2)Bii deltaP[i]=Pi(sch) -Pi; deltaQ[i]=Qi(sch) -Qi; next i

Program Listing_3: NRLF by Van Ness Algorithm - no delete to solve Ax=b

1: /* Program  NRLF3 -powerflow by Newton-Raphson using Van Ness Algorithm
   2:    1. contiguous bus numbering
   3:    2. full matrix               
   4:    3. single sweep on line data to form the Jacobian
   5:    4. Slack, PV, randomly located
   6:    5. Ybus, linear storage   
   7:    6. No delete of row/col on the Jacobian to perform Ax=b
   8:    Date: Nov 11, 1999  Vancouver,BC, Canada
   9:    file: nrlfc3.cpp
  10:    Sparsity Technology, Joesel Pante
  11:    http://www.geocities.com/siliconvalley/lab/4223   for details
  12: */
  13: #define LMAX 20
  14: #define BMAX  15  
  15: #define BMAX2 30
  16: #include <stdio.h>
  17: #include <math.h>
  18: 
  19: int nl,nb,size,np[LMAX],nq[LMAX];  //p,q bus    
  20: float base=100;  //base Mva 
  21: float gii[BMAX],bii[BMAX],gij[LMAX],bij[LMAX];   //Ybus,linear storage
  22: int bus[BMAX],type[BMAX],tcol[BMAX2],tnum[BMAX2];     
  23: //type[] is used to guide lufact on Ax=b, as A contains Slack and PV bus.
  24: //without having to delete row-column.
  25: 
  26: float Pg[BMAX],Qg[BMAX],Pd[BMAX],Qd[BMAX],Qgmin[BMAX],Qgmax[BMAX];    
  27: float P[BMAX],Q[BMAX];
  28: float yc[LMAX],Pij[LMAX],Qij[LMAX]; //line flow 
  29: float Es[BMAX],e[BMAX],f[BMAX];   //Es[]=voltage for PV    
  30: float A[BMAX2][BMAX2],x[BMAX2];  //jacobian                        
  31: FILE *file;
  32: 
  33: void Enter(void)
  34: {char *ch; printf("Enter to continue\n");gets(ch);}
  35:           
  36: 
  37: void input(void)
  38: { int k,m,i;     
  39:   float r,x,ch,a1;     
  40:   nb=0,nl=0; //initialize counters nl=lines, nb=bus 
  41:   for (i=0;i<BMAX;i++) 
  42:    {gii[i]=0;       //initialize Yii
  43:     bii[i]=0; 
  44:    } 
  45:   printf(" count  p   q     r      x      MVAR\n");
  46:   while (1)
  47:    { nl++; 
  48:      fscanf(file,"%d %d %f %f %f",&k,&m,&r,&x,&ch);
  49:      if (k==99) break;
  50:      printf("%5d %3d %3d %7.2f %7.2f %7.2f\n",nl,k,m,r,x,ch);
  51:      if (k>nb) nb=k; if (m>nb) nb=m;  /*get the largest bus*/    
  52:      np[nl]=k; nq[nl]=m;
  53:      a1=r*r+x*x;
  54:      gij[nl]=-r/a1;   //off-diagonal
  55:      bij[nl]=x/a1; 
  56:      yc[nl]=ch/200;    //charging in admittance is +j
  57:      gii[k]-=gij[nl]; //diagonal
  58:      gii[m]-=gij[nl];   
  59:      bii[k]=bii[k]-bij[nl]+yc[nl];    //charging to pi and pu
  60:      bii[m]=bii[m]-bij[nl]+yc[nl];
  61:       
  62:     }
  63:    nl--;     //correct nl=lines read  
  64:    printf("number of lines:%2d  number of nodes:%2d\n",nl,nb);  
  65:    nb=0; 
  66:    printf("count Bus Type    E      Pg      Qg      Pd      Qd     Qgmin    Qgmax\n");
  67:    while (1)
  68:    { nb++; 
  69:      fscanf(file,"%d %d %f %f %f %f %f %f %f",
  70:      bus+nb,type+nb,e+nb,Pg+nb,Qg+nb,Pd+nb,Qd+nb,Qgmin+nb,Qgmax+nb); 
  71:      Es[nb]=e[nb];    f[nb]=0;  //zero the j part of E
  72:      if (bus[nb]==99) break;
  73:      printf("%5d %3d %3d %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f\n",
  74:      nb,bus[nb],type[nb],e[nb],Pg[nb],Qg[nb],Pd[nb],Qd[nb],Qgmin[nb],Qgmax[nb]); 
  75:      Pg[nb]/=base;Qg[nb]/=base;Pd[nb]/=base;Qd[nb]/=base;Qgmin[nb]/=base;Qgmax[nb]/=base;
  76:    }   
  77:    nb--;  
  78: }
  79:  
  80:  
  81: void printY(void)  
  82: {int j; 
  83:  printf("count np  nq  gij[]  bij[]\n");  
  84:  for (j=1;j<=nl;j++)
  85:   printf("%4d %4d %4d %9.4f %9.4f\n",j,np[j],nq[j],gij[j],bij[j]); 
  86:  printf("count gii[]  bii[]\n");  
  87:  for (j=1;j<=nb;j++)
  88:   printf("%4d  %9.4f %9.4f\n",j,gii[j],bii[j]); 
  89:   
  90:  Enter();
  91: }     
  92: 
  93: void printA(void)  
  94: {int i,j; 
  95:  printf("The Jacobian:\n");  
  96:  for (i=1;i<=size;i++) 
  97:    printf("%9d",tcol[i]);printf("\n"); 
  98:  for (i=1;i<=size;i++) 
  99:    printf("%9d",tnum[i]);printf("\n");      
 100:  for (i=1;i<=size;i++)
 101:    {printf("%2d ",tcol[i]); 
 102:     printf("%2d ",tnum[i]);  
 103:     for (j=1;j<=size;j++)
 104:      printf("%9.4f",A[i][j]); 
 105:     printf("%9.4f\n",x[i]);  
 106:    }   
 107:  printf("\n"); 
 108:  Enter();
 109: }               
 110:  
 111: //Solution of simultaneous equations Ax=b
 112: //by unsymmetric LU Factorization              
 113: void lufact(void)
 114: {int i,j,k;
 115:  for (k=1; k<size; k++)  
 116:    if (tcol[k]!=2)          
 117: //place here check for PV on even number only
 118: //if (tcol[k]==1 && !(k%2))  exit   
 119: // k%2=0 for even number, =1 for odd number   
 120:      {for (i=k+1; i<=size; i++)  //solve for L     
 121:         A[i][k]/=A[k][k];
 122:       for (i=k+1; i<=size; i++)   //block    
 123:         if (tcol[i]!=2)
 124:           for (j=k+1; j<=size; j++)    
 125:              if (tcol[k]!=2)
 126:                 A[i][j]=A[i][j]-A[k][j]*A[i][k];
 127:         
 128:      }
 129:   printf("After factoring\n");  
 130:   printA();
 131: //forward course
 132:  for (k=1;k<size; k++) 
 133:    if (tcol[k]!=2)
 134:     for (i=k+1; i<=size; i++)      //does not contain tcol[]
 135:       x[i]=x[i]-A[i][k]*x[k]; 
 136: //backward course
 137:  x[size]/=A[size][size];
 138:  for (k=size-1;k>=1; k--)   
 139:    {if (tcol[k]!=2)
 140:     for (j=k+1;j<=size;j++)     
 141:       if (tcol[j]!=2)
 142:         x[k]=x[k]-A[k][j]*x[j];   
 143:     x[k]/=A[k][k];         
 144:    }                    
 145: //
 146:   printf("After forward and backward course\n");                      
 147:   printA();
 148: }             
 149:   
 150: void nrlfc1(void)
 151: {int iteration=0,itermax=1,i,j,k;
 152:  float maxPQ=0.01,delPQ;
 153:  float AIJ,AJI,BIJ,BJI,ReEYEij,ReEYEji,ImEYEij,ImEYEji,E2,E1,d1,deltaE,deltaO;                    
 154: //form the Jacobian matrix by single sweep on the line data 
 155:  size=nb*2;  //outside the loop, set only once  
 156: //replicate type[1->nb] to type[nb+1->2*nb]  
 157:   for (i=1;i<=nb;i++)
 158:    {tcol[i*2]=type[i];tcol[i*2-1]=type[i];    
 159:     tnum[i*2]=i;tnum[i*2-1]=i;
 160:    }   
 161:  do
 162:    {delPQ=0;   //set max error
 163:     iteration++; 
 164: //zero out the square Jacobian matrix first 
 165:     for (i=1; i<=size; i++)
 166:       for (j=1; j<=size;j++)
 167:          A[i][j]=0;     
 168:     for (i=1;i<=nb;i++)  //P,Q accumulator
 169:       {P[i]=0;
 170:        Q[i]=0;
 171:       }
 172:            
 173:     for (k=1;k<=nl;k++)     //single sweep on line data to form HJNL
 174:      {i=np[k];j=nq[k];
 175: //Hij= 
 176:       AIJ=e[j]*gij[k]-f[j]*bij[k]; AJI=e[i]*gij[k]-f[i]*bij[k];   
 177:       
 178:       BIJ=e[j]*bij[k]+f[j]*gij[k]; BJI=e[i]*bij[k]+f[i]*gij[k];  
 179:       
 180:       ReEYEij=e[i]*AIJ+f[i]*BIJ; ReEYEji=e[j]*AJI+f[j]*BJI;      
 181:       
 182:       ImEYEij=e[i]*BIJ-f[i]*AIJ; ImEYEji=e[j]*BJI-f[j]*AJI;      
 183:       
 184:       P[i]+=ReEYEij;P[j]+=ReEYEji;    
 185:       Q[i]+=ImEYEij;Q[j]+=ImEYEji;
 186:       
 187:       A[i*2-1][j*2-1]= -ImEYEij;    //Hij    
 188:       A[j*2-1][i*2-1]= -ImEYEji;    //Hji
 189:       A[i*2-1][j*2  ]=  ReEYEij;    //Nij  
 190:       A[j*2-1][i*2  ]=  ReEYEji;    //Nji    
 191:       A[i*2  ][j*2-1]= -ReEYEij;    //Jij
 192:       A[j*2  ][i*2-1]= -ReEYEji;    //Jji
 193:       A[i*2  ][j*2  ]= -ImEYEij;    //Lij=Hij
 194:       A[j*2  ][i*2  ]= -ImEYEji;    //Lji=Hji
 195:      }
 196:   
 197: //include j=i. Loop on bus data, 
 198:     for (i=1;i<=nb;i++)
 199:      {A[i*2-1][i*2-1]=Q[i];       //Hii, sum Hij, sum ImEYEij  
 200:       A[i*2  ][i*2-1]=P[i];       //Jii= sum ReEYEij
 201:       E2=e[i]*e[i]+f[i]*f[i];     //E*E  
 202:       P[i]+=E2*gii[i];
 203:       Q[i]= -E2*bii[i] -Q[i];
 204:       A[i*2-1][i*2]=P[i] +E2*gii[i];        //Nii
 205:       A[i*2  ][i*2]=Q[i] -E2*bii[i];        //Lii
 206: //Delta P/Q      
 207:       x[i*2-1]=Pg[i]-Pd[i]-P[i];   
 208:       x[i*2  ]=Qg[i]-Qd[i]-Q[i];  
 209:       if (type[i]!=2) 
 210:         {if((float)fabs(x[i*2-1])>delPQ) delPQ=(float)fabs(x[i*2-1]);  //0,1     
 211:          if (!type[i]) 
 212:           if((float)fabs(x[i*2-1])>delPQ) delPQ=(float)fabs(x[i*2-1]); //0    
 213:         }  
 214:      }            
 215:      printf("iter:%2d max delPQ:%8.6f\n",iteration,delPQ);
 216:      printA(); 
 217:      lufact(); 
 218: //update voltage angle and magnitude      
 219:     for (i=1;i<=nb;i++)     
 220: //convert rectangular to polar for polar update   //there must be a better way   
 221:      {E1=(float)(sqrt(e[i]*e[i]+f[i]*f[i]));  
 222:       d1=(float)(atan(f[i]/e[i]));
 223:       e[i]=E1;
 224:       f[i]=d1;
 225:       deltaE=x[i*2]*e[i];
 226:       deltaO=x[i*2-1];
 227:       if (!type[i]) e[i]+=deltaE;    //update voltage for bus PQ only
 228:       if (type[i]!=2) f[i]+=deltaO;  //update theta for PQ and PV
 229: //convert from polar to rectangular again   
 230:       printf("i:%4d E=%8.5f deltaE=%8.5f O=%8.5f deltaO=%8.5f \n",
 231:       i,e[i],deltaE,f[i],deltaO);  
 232:       E1=e[i];
 233:       d1=f[i];
 234:       e[i]=E1*(float)(cos(d1));
 235:       f[i]=E1*(float)(sin(d1));
 236:      } 
 237:       
 238:   }while (iteration<=itermax && delPQ>maxPQ);   
 239: //We can get the Pg/Qg of slack, Qg of PV from P[i]/Q[i] here,        
 240: //the problem is that the voltages are not updated yet
 241: 
 242: //convert to polar Voltages
 243:  for (i=1;i<=nb;i++) 
 244:    { E1=(float)(sqrt(e[i]*e[i]+f[i]*f[i]));  
 245:     d1=(float)(atan(f[i]/e[i]));
 246:     e[i]=E1;                          //use e[] instead of E[]
 247:     f[i]=d1;                          //use f[] instead of d[]
 248:    } 
 249:    
 250: }                    
 251: 
 252: 
 253: void pflows(void)
 254: {int i,k,m;    
 255:  float angle,co,si,Gij,Bij;
 256:  for (i=1; i<=nl; i++) 
 257:   {k=np[i];m=nq[i]; 
 258:    angle=f[k]-f[m];       //d[i]
 259:    co=(float)cos(angle);
 260:    si=(float)sin(angle); 
 261:    Gij=-gij[i];Bij=-bij[i];  //from Ybus
 262:    Pij[i]=Gij*e[k]*e[k]-e[k]*e[m]*(Gij*co+Bij*si);
 263:    Qij[i]=-(yc[i]+Bij)*e[k]*e[k]-e[k]*e[m]*(Gij*si-Bij*co);  //E[]
 264:    printf("i:%2d j:%2d Pij:%7.2f Qij:%7.2f\n",k,m,Pij[i]*base,Qij[i]*base);
 265:   }
 266:  }  
 267:  
 268:  
 269: //Ybus is not a square matrix. Doesn`t negate Gij,Bij
 270: void gener(void)  
 271: {int i,j,m;     
 272:  float co,si;
 273:  for (i=1; i<=nb; i++)
 274:   if (type[i]==2)
 275:     {Pg[i]=Qg[i]=0; //solve Pg,Qg   
 276: //line data, off diagonals only    
 277:      for (m=1;m<=nl; m++)             
 278:        {if (np[m]==i)
 279:           {j=nq[m];
 280:            co=(float)cos(f[i]-f[j]);
 281:            si=(float)sin(f[i]-f[j]);  
 282:            Pg[i]=Pg[i]+e[j]*(gij[m]*co+bij[m]*si);
 283:            Qg[i]=Qg[i]+e[j]*(gij[m]*si-bij[m]*co); 
 284:           }  
 285:         else if (nq[m]==i)
 286:           {j=np[m];
 287:            co=(float)cos(f[i]-f[j]);
 288:            si=(float)sin(f[i]-f[j]);  
 289:            Pg[i]=Pg[i]+e[j]*(gij[m]*co+bij[m]*si);
 290:            Qg[i]=Qg[i]+e[j]*(gij[m]*si-bij[m]*co); 
 291:           } 
 292:        }                                 
 293: //include diagonal, j=i       
 294:      Pg[i]+=e[i]*gii[i];
 295:      Qg[i]-=e[i]*bii[i]; 
 296: //Ei before the summation formula       
 297:      Pg[i]*=e[i];
 298:      Qg[i]*=e[i];  
 299: //add load
 300:      Pg[i]+=Pd[i];
 301:      Qg[i]+=Qd[i];     
 302:      printf("Slack Pg[%2d]=%7.2f Qg=%7.2f\n",i,Pg[i]*base,Qg[i]*base);  
 303:     }
 304:    else if (type[i]==1) 
 305:     {Qg[i]=0;  //solve Qg
 306:      for (m=1;m<=nl; m++)          //includes j=i
 307:         {if (np[m]==i)
 308:           {j=nq[m];
 309:            co=(float)cos(f[i]-f[j]);
 310:            si=(float)sin(f[i]-f[j]); 
 311:            Qg[i]=Qg[i]+e[j]*(gij[m]*si-bij[m]*co); 
 312:           }   
 313:          else if (nq[m]==i)
 314:           {j=np[m];
 315:            co=(float)cos(f[i]-f[j]);
 316:            si=(float)sin(f[i]-f[j]); 
 317:            Qg[i]=Qg[i]+e[j]*(gij[m]*si-bij[m]*co); 
 318:           }  
 319:         }  
 320:  //include diagonal, j=i       
 321:      Qg[i]-=e[i]*bii[i]; 
 322: //Ei before the summation formula         
 323:      Qg[i]*=e[i]; 
 324:      Qg[i]+=Qd[i];
 325:      printf("PV Qg[%2d]=%7.2f\n",i,Qg[i]*base);       
 326:     }
 327: }
 328:                   
 329: void main(void)
 330: {int i;
 331:  file=fopen("nrlf.dat","r");
 332:  input();            //read line,bus data
 333:  printY();           //print Ybus=g[][]+jb[][];
 334:  nrlfc1();           //iterate voltage
 335:  pflows();           //solve power flows
 336:  gener();            //solve generation, Pg,Qg=SL, Qg=PV
 337:  i=fclose(file);
 338: }

Code Explanation: This code simulates the algorithm used by a full blown Newton Raphson made by PECo. We can learn plenty of things by analyzing how a successful program implements area control, tap-changing under load transformer, PV busses hitting its limits and others.

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