How Can I Do a Power Flow Analysis of a Simple Circuit?
The other day our professor was teaching nodal analysis and the admittance matrix in the power system during the lecture.
The mathematics behind was largely linear algebra, which is itself an interesting topic to learn.
But these are just useless jargon unless I could really write an algorithm that can run a power flow analysis with them. So I decided to model a simple circuit step by step. The circuit was given at class as an example:
The following content contains some codes of R after the theorem is explained. These sections of codes can be skipped if you are only interested in knowing the theorem behind the algorithm.
If you code with other languages, the codes provided are still easy to adapt because they require no additional library in R.
The Topology (Network) Matrix
Before moving on, we need to first know some basics of topology.
A node is where at least three lines meet each other.
A line, on the other hand, always connects a node with another node.
To describe the relationship between the lines and the nodes, the topology matrix is introduced.
Here is how the topology matrix [N] is constructed.
For every line, we define its direction. Then, entry [N]_ij of the matrix [N] would be 1 if the i-th arrow starts from the j-th node, -1 if it ends at the j-th node, 0 if otherwise. The direction of the arrows are arbitrary.
This implies that for a system with m lines and n nodes, the dimension of [N] will be m*n.
For example, the topology matrix [N] of our circuit would be:
The Line Admittance Matrix
In a power system, we can assign each line a value of its total admittance. We can then form a m*m diagonal matrix [Y_l], where [Y_l]_ii is the total admittance of the i-th line.
The lines which does not connect two known nodes in the system or does not have an assigned admittance, we do not include them in our system.
For our circuit, the line admittance matrix would be:
The Additional Current Vector
For each node, there might be additional current coming from outside the system. These additional currents are represented in the additional current vector {I_0}.
To tell whether a current is additional or not, see if the line it comes from (or leaves) is inside the system.
The direction of the additional must follow the same denotation as the topology matrix. In our case, a positive additional current means the current is moving out from the node, and a negative additional current means the current is moving into the node.
In our circuit above there is no additional current. Our additional current vector is thus:
The Additional Voltage Vector
For each line, there might be an additional voltage change between two nodes other than the one caused by impedance. This is usually the case when there is a generator or generator on the line. We represent it in the additional voltage vector {U_0}.
The sign of this additional voltage depends on the direction of arrow we have assigned earlier to the line. If the positive pole of the battery or generator heads toward the direction of the arrow of a line, then the voltage change is positive.
In our circuit the opposite happens: all three additional voltage are thus negative. So our additional voltage vector is:
In nodal analysis, only nodes can have assigned voltage values, so additional voltages should be turned into additional currents via Norton-Thevenin theorem.
If you don’t do so, you will still be able to do the power flow analysis, but you will not be able to calculate the current of each line via the matrix form of Ohm’s Law.
The Current and Voltage Laws in Matrix Form
After applying the Norton-Thevenin theorem, the Kirchhoff’s current law of each node can be represented as:
Meanwhile, the Ohm’s Law for every line can also be written as a matrix form:
Substituting the {I} in equation(1) into (2), we will quickly obtain:
In equation(3), only the voltage of the nodes are unknown, so this forms a classic linear algebra problem with the type of Ax=b.
But the set of equations is not solvable. It can be shown that there are only n-1 degrees of freedom in the system. Therefore we conveniently set the voltage of the last node to be zero.
The Node (Bus) Admittance Matrix
Equation(3) is sometimes written as
Where
and
[Y_eq] is sometimes called the bus (node) admittance matrix, and {I_eq} is sometimes called the equivalent current.
The node admittance matrix actually has many beautiful properties. For example, in our case the bus admittance matrix is:
(Only the first 2*2 block is shown, since it is the only relevant part in solving the equations.)
At first glimpse we can see that the matrix is symmetric. It turns out to be the basic property of the admittance matrix.
Two other important properties of this matrix:
- The i-th diagonal entry is the summation of the admittance of the lines that has an end at the i-th node.
- The #ij entry off the diagonal is minus one times the admittance of the line which connects the i-th and j-th nodes.
Test of Algorithm: Nodal Analysis
Now all the building blocks of the power flow analysis are ready, we can test if it really works.
The following code executes the algorithm:
##Functions
#Gaussian Seidel Iteration
GS<-function(A,b){
error=1
x_0=c()
for(i in 1:length(b)){
x_0[i]=0
}
x=x_0
while(error>10^(-6)){
sum=c()
for(i in 1:(length(b)-1)){
sum=0
for(j in 1:(length(b)-1)){
if(i!=j){
sum=sum+A[i,j]*x[j]
}
}
x[i]=(1/diag(A)[i])*(b[i]-sum)
}
error=sum(abs((A%*%x-b)*Conj(b)))/sum((abs(b))^2)
}
return(x)
}#Solve for U
U<-function(N,Y,I_0,U_0){
Y_eq=t(N)%*%Y%*%N
I_eq=-t(N)%*%Y%*%U_0-I_0
U=GS(Y_eq,I_eq)
return(U)
}
The below is the set parameters:
##Parameters
N=matrix(c(1,0,1,1,1,-1,1,0,0,-1,0,-1,-1,-1,0),ncol=3)
Y=complex(imaginary=1)*diag(c(2,1,1,-.5,.5))
I_0=c(0,0,0)
U_0=sample(0,size=nrow(N),replace=T)
U_0[2:4]=c(2,3,5)
U_sol=U(N,Y,I_0,U_0)The resulting voltages on node 1 and 2 are {-1.59,-1.71} (pu).
From Nodal Analysis to Power Flow Analysis
In our example circuit, we can calculate directly the voltage vector {U} when the parameters in[Y_eq] and {U_0} are given by solving the linear system.
This means that we must know the additional voltage and impedance of every line. However, it is usually the case that we only know the following details on the nodes:
Condition one: the apparent power output of each node {S}={P}+j{Q} (a load bus)
Condition two: the real power {P} and the magnitudes of the voltage {U} of each node (a generator bus)
This means that instead of solving the system explicitly, we have to implicitly find a set of voltage such that the given conditions are met.
For the first condition, we
- first guess a voltage vector, {U_guess}.
- can get {I_eq,guess} with equation(4).
- compute the guess power using equation(8).
- see how much this guess value deviates from the exact value and use Newton-Raphson method to get the next iteration of {U_guess}. That is, to solve equation(9).
Some Extra Notes on the Jacobian
In equation(9), [J] is the Jacobian matrix:
The partial derivatives of real and reactive power deviations with respect to magnitude of voltage and phase angle are (the Y here is actually the bus admittance matrix, denoted as Y_eq previously):
Note that the two sets of iterative variables in this algorithm are the magnitudes and phase angles of the voltage at each node.
For the second condition, it is tempting to think that we can directly plug in {U} to solve the system of equations.
However, we do not know the relative phase angle of the voltages, so we still have to take a guess of the phase angles and construct our {U_guess}. We also need to guess the reactive power {Q_guess}.
Then we get to use Newton-Raphson method again, this time taking the phase angles and reactive power as iterating variables. The partial derivative of the deviation of reactive power and Q_guess is -1, so the Jacobian matrix of equation (10) can be transformed quite easily.
Test of Algorithm: Power Flow Analysis
Consider our circuit with known power output at node 1 and 2.
The following code is given for Newton-Raphson algorithm:
#Gaussian Seidel Iteration for Power
GS<-function(A,b){
error=1
x_0=c()
for(i in 1:length(b)){
x_0[i]=0
}
x=x_0
counter=0
while(error>10^(-6)){
sum=c()
for(i in 1:length(b)){
sum=0
for(j in 1:length(b)){
if(i!=j){
sum=sum+A[i,j]*x[j]
}
}
x[i]=(1/diag(A)[i])*(b[i]-sum)
}
error=sum(abs((A%*%x-b)*Conj(b)))/sum((abs(b))²)
}
return(x)
}The iterative functions for the first condition is the following.
#Jacobian, condition 1
Jacobian<-function(Y_eq,U_old,phi_old){
J=matrix(nrow=2*length(U_old),ncol=2*length(U_old))
for(i in 1:length(U_old)){
#partial derivatives with respect to phi
#P
J[i,1:(ncol(J)/2)]=U_old[i]*U_old*(Re(Y_eq[i,])*sin(phi_old[i]-phi_old)-Im(Y_eq[i,])*cos(phi_old[i]-phi_old))
#Q
J[i+(ncol(J)/2),1:(ncol(J)/2)]=U_old[i]*U_old*(-Re(Y_eq[i,])*cos(phi_old[i]-phi_old)-Im(Y_eq[i,])*sin(phi_old[i]-phi_old))
#partial derivatives with respect to U
#P
J[i,(ncol(J)/2+1):ncol(J)]=U_old[i]*(Re(Y_eq[i,])*cos(phi_old[i]-phi_old)+Im(Y_eq[i,])*sin(phi_old[i]-phi_old))
J[i,(ncol(J)/2)+i]=2*U_old[i]*Re(Y_eq[i,i])
#Q
J[i+(ncol(J)/2),(ncol(J)/2+1):ncol(J)]=U_old[i]*(Re(Y_eq[i,])*sin(phi_old[i]-phi_old)-Im(Y_eq[i,])*cos(phi_old[i]-phi_old))
J[i+(ncol(J)/2),(ncol(J)/2)+i]=-2*U_old[i]*Im(Y_eq[i,i])
}
return(J)
}#Solve for U, condition 1
Var<-function(N,Y,S,U_old,phi_old){
U=U_old*complex(real=cos(phi_old),imaginary=sin(phi_old))
Y_eq=t(N)%*%Y%*%N
I_eq=Y_eq%*%U
delta_S=c()
for(i in 1:nrow(Y_eq)){
delta_S[i]=U[i]*Conj(I_eq[i])-S[i]
}
if(sum(delta_S^2)==0){
var=c(U_old,phi_old)
}
else{
J=Jacobian(Y_eq[1:(nrow(Y_eq)-1),1:(ncol(Y_eq)-1)],U_old[1:(nrow(Y_eq)-1)],phi_old[1:(nrow(Y_eq)-1)])
delta_var=GS(-J,c(Re(delta_S)[1:(nrow(Y_eq)-1)],Im(delta_S)[1:(nrow(Y_eq)-1)]))
phi_new=phi_old+c(delta_var[1:(nrow(Y_eq)-1)],0)
U_new=U_old+c(delta_var[nrow(Y_eq):length(delta_var)],0)
var=c(U_new,phi_new)
}return(var)
}
The iterative function for the second method is:
#Jacobian, condition 2
Jacobian<-function(Y_eq,U_old,phi_old){
J=matrix(nrow=2*length(U_old),ncol=2*length(U_old))
for(i in 1:length(U_old)){
#partial derivatives with respect to phi
#P
J[i,1:(ncol(J)/2)]=U_old[i]*U_old*(Re(Y_eq[i,])*sin(phi_old[i]-phi_old)-Im(Y_eq[i,])*cos(phi_old[i]-phi_old))
#Q
J[i+(ncol(J)/2),1:(ncol(J)/2)]=U_old[i]*U_old*(-Re(Y_eq[i,])*cos(phi_old[i]-phi_old)-Im(Y_eq[i,])*sin(phi_old[i]-phi_old))
#partial derivatives with respect to Q
#P
J[i,(ncol(J)/2+1):ncol(J)]=0
J[i,(ncol(J)/2)+i]=-1
#Q
J[i+(ncol(J)/2),(ncol(J)/2+1):ncol(J)]=0
J[i+(ncol(J)/2),(ncol(J)/2)+i]=-1
}
return(J)
}#Solve for U, condition 2
Var<-function(N,Y,P,V,Q_old,phi_old){
U_0=V*complex(real=cos(phi_old),imaginary=sin(phi_old))
Y_eq=t(N)%*%Y%*%N
I_eq=Y_eq%*%U_0
S=complex(real=P,imaginary=Q_old)
delta_S=c()
for(i in 1:nrow(Y_eq)){
delta_S[i]=U_0[i]*Conj(I_eq[i])-S[i]
}
if(sum(delta_S^2)==0){
var=c(Q_old,phi_old)
}
else{
J=Jacobian(Y_eq[1:(nrow(Y_eq)-1),1:(ncol(Y_eq)-1)],V[1:(nrow(Y_eq)-1)],phi_old[1:(nrow(Y_eq)-1)])
delta_var=GS(-J,c(Re(delta_S)[1:(nrow(Y_eq)-1)],Im(delta_S)[1:(nrow(Y_eq)-1)]))
phi_new=phi_old+c(delta_var[1:(nrow(Y_eq)-1)],0)
Q_new=Q_old+c(delta_var[nrow(Y_eq):length(delta_var)],0)
var=c(Q_new,phi_new)
}
return(var)
}
For the first condition, the following parameters is computed:
##Parameters for Power Flow, condition 1
P=c(0,0,0)
Q=c(-.1,-.1,0)
S=complex(real=P,imaginary=Q)
N=matrix(c(1,0,1,1,1,-1,1,0,0,-1,0,-1,-1,-1,0),ncol=3)
Y=complex(imaginary=1)*diag(c(1,1,1,-.5,.5))The iterative codes:
TIME=1000
U_guess=matrix(nrow=TIME,ncol=ncol(N))
U_guess[1,]=c(1,1,0)
phi_guess=matrix(nrow=TIME,ncol=ncol(N))
phi_guess[1,]=c(0,0,0)
var=c()
for(i in 2:TIME){
var=Var(N,Y,S,U_guess[i-1,],phi_guess[i-1,])
U_guess[i,]=var[1:(length(var)/2)]
phi_guess[i,]=var[(length(var)/2+1):length(var)]
}The resulting voltages for node 1 and 2 are {0.3897342, 0.3485889} (pu).
For the second condition, the following parameters is computed:
##Parameters for Power Flow, condition 2
P=c(0,0,0)
U=c(.38,.35,0)
N=matrix(c(1,0,1,1,1,-1,1,0,0,-1,0,-1,-1,-1,0),ncol=3)
Y=complex(imaginary=1)*diag(c(1,1,1,-.5,.5))TIME=1000
Q_guess=matrix(nrow=TIME,ncol=ncol(N))
Q_guess[1,]=c(1,1,0)
phi_guess=matrix(nrow=TIME,ncol=ncol(N))
phi_guess[1,]=c(0,0,0)
var=c()
for(i in 2:TIME){
var=Var(N,Y,P,U,Q_guess[i-1,],phi_guess[i-1,])
Q_guess[i,]=var[1:(length(var)/2)]
phi_guess[i,]=var[(length(var)/2+1):length(var)]
}
The resulting reactive powers for node 1 and node 2 are {-0.08930, -0.10675} (pu). The convergence rate is much faster than condition 1; it only took one iteration to get to the exact solution.
Ending Notes
- In real power systems, both P-Q and P-U conditions co-exist, so the Jacobian matrix, the initial parameters, and the iterative variables must be set according.
- If the calculated reactive power of a P-U bus exceed the capacity of the bus, we set the reactive power to the limit and turn it into a P-Q bus.
- The solutions may not be unique; in the first case of the power flow analysis, there probably are four non-trivial solutions. However, only stable solutions can be found via iterative methods.
- There are other methods to compute the power flow directly, the most promising one is perhaps the Holomorphic embedding load flow method.
About This Series…
This article is part of the series My Coding Diary with R. To know how the series has been developing, check the intro article out:
