Chapter 4: Load Flow Studies
Forming Ybus Matrix
This is a function that can be called by various programs. The function can be invoked by the statement
[yb,ych]=ybus;
where 'yb' and 'ych' are respectively the Ybus matrix and a matrix containing the line charging admittances. It
is assumed that the system data of Table 4.1 are given in matrix form and the matrix that contains line
impedances is 'zz', while 'ych' contains the line charging information. This program is stored in the file
ybus.m. The program listing is given below.
% Function ybus
% THIS IS THE PROGRAM FOR CREATING Ybus MATRIX.
function [yb,ych]=ybus
% The line impedances are
zz=[0 0.02+0.1i 0 0 0.05+0.25i
0.02+0.1i 0 0.04+0.2i 0 0.05+0.25i
0 0.04+0.2i 0 0.05+0.25i 0.08+0.4i
0 0 0.05+0.25i 0 0.1+0.5i
0.05+0.25i 0.05+0.25i 0.08+0.4i 0.1+0.5i 0];
% The line chargings are
ych=j*[0 0.03 0 0 0.02
0.03 0 0.025 0 0.020
0 0.025 0 0.02 0.01
0 0 0.02 0 0.075
0.02 0.02 0.01 0.075 0];
% The Ybus matrix is formed here
for i=1:5
for j=1:5
if zz(i,j) == 0
yb(i,j)=0;
else
yb(i,j)=-1/zz(i,j);
end
end
end
for i=1:5
ysum=0;
csum=0;
for j=1:5
ysum=ysum+yb(i,j);
csum=csum+ych(i,j);
end
yb(i,i)=csum-ysum;
end
Chapter 4: Load Flow Studies
Gauss-Seidel Load Flow
The Gauss-Seidel program is stored in the file loadflow_gs.m. This calls the ybus.m function
discussed above. The program allows the selection of the acceleration factor. The program lists the
number of iterations required to converge, bus voltages and their magnitudes and real and reactive
power. The program listing is given below.
% Program loadflow_gs
% THIS IS A GAUSS-SEIDEL POWER FLOW PROGRAM
clear all
d2r=pi/180;w=100*pi;
% The Y_bus matrix is
[ybus,ych]=ybus;
g=real(ybus);b=imag(ybus);
% The given parameters and initial conditions are
p=[0;-0.96;-0.35;-0.16;0.24];
q=[0;-0.62;-0.14;-0.08;-0.35];
mv=[1.05;1;1;1;1.02];
th=[0;0;0;0;0];
v=[mv(1);mv(2);mv(3);mv(4);mv(5)];
acc=input('Enter the acceleration constant: ');
del=1;indx=0;
% The Gauss-Seidel iterations starts here
while del>1e-6
% P-Q buses
for i=2:4
tmp1=(p(i)-j*q(i))/conj(v(i));
tmp2=0;
for k=1:5
if (i==k)
tmp2=tmp2+0;
else
tmp2=tmp2+ybus(i,k)*v(k);
end
end
vt=(tmp1-tmp2)/ybus(i,i);
v(i)=v(i)+acc*(vt-v(i));
end
% P-V bus
q5=0;
for i=1:5
q5=q5+ybus(5,i)*v(i);
end
q5=-imag(conj(v(5))*q5);
tmp1=(p(5)-j*q5)/conj(v(5));
tmp2=0;
for k=1:4
tmp2=tmp2+ybus(5,k)*v(k);
end
vt=(tmp1-tmp2)/ybus(5,5);
v(5)=abs(v(5))*vt/abs(vt);
% Calculate P and Q
for i=1:5
sm=0;
for k=1:5
sm=sm+ybus(i,k)*v(k);
end
s(i)=conj(v(i))*sm;
end
% The mismatch
delp=p-real(s)';
delq=q+imag(s)';
delpq=[delp(2:5);delq(2:4)];
del=max(abs(delpq));
indx=indx+1;
if indx==1
pause
end
end
'GS LOAD FLOW CONVERGES IN ITERATIONS',indx,pause
'FINAL VOLTAGE MAGNITUDES ARE',abs(v)',pause
'FINAL ANGLES IN DEGREE ARE',angle(v)'/d2r,pause
'THE REAL POWERS IN EACH BUS IN MW ARE',(real(s)+[0 0 0 0
0.24])*100,pause
'THE REACTIVE POWERS IN EACH BUS IN MVar ARE',(-imag(s)+[0 0 0 0
0.11])*100
Chapter 4: Load Flow Studies
Solving Nonlinear Equations using Newton-Raphson
This program gives the solution of the nonlinear equations of Example 4.1 using the NewtonRaphson method. The equations and the Jacobian matrix are explicitly entered in the program itself.
The program gives the number of iterations and the final values of x1 , x2 and x3 . The program listing
is given below.
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
Program nwtraph
THIS IS A NEWTON-RAPHSON PROGRAM
We have to solve three nonlinear equations given by
g1=x1^2-x2^2+x3^2-11=0
g2=x1*x2+x2^2-3x3-3=0
g3=x1-x1*x3+x2*x3-6=0
Let us assume the initial conditions of x1=x2=x3=1
The Jacobian matrix is
J=[2x1 -2x2
x2
x1+2x2
1-x3
x3
2x3
-3
-x1+x2];
clear all
x=[1;1;1];
% The Newton-Raphson iterations starts here
del=1;
indx=0;
while del>1e-6
g=[x(1)^2-x(2)^2+x(3)^2-11;x(1)*x(2)+x(2)^2-3*x(3)-3;
x(1)-x(1)*x(3)+x(2)*x(3)-6];
J=[2*x(1) -2*x(2) 2*x(3);x(2) x(1)+2*x(2) -3;1-x(3) x(3)
-x(1)+x(2)];
delx=-inv(J)*g;
x=x+delx;
del=max(abs(g));
indx=indx+1;
end
'NEWTON-RAPHSON SOLUTION CONVERGES IN
ITERATIONS',indx,pause
'FINAL VALUES OF x ARE',x
Chapter 4: Load Flow Studies
Newton-Raphson Load Flow
The Newton-Raphson load flow program is stored in the files loadflow_nr.m. The outputs of the
program can be checked by typing
indx
the number of iterations
bus voltages in Cartesian form
abs(v)
magnitude of bus voltages
angle(v)/d2r
angle of bus voltage in degree
preal
real power in MW
preac
reactive power in MVAr
pwr
power flow in the various line segments
qwr
reactive power flow in the various line segments
reactive power entering or leaving a bus
pl
real power losses in various line segments
ql
reactive drops in various line segments
It is to be noted that in calculating the power and reactive power the conventions that the power
entering a node is positive and leaving it is negative are maintained. The program listing for the
Newton-Raphson load flow is given below.
% Program loadflow_nr
% THIS IS THE NEWTON-RAPHSON POWER FLOW PROGRAM
clear all
d2r=pi/180;w=100*pi;
% The Ybus matrix is
[ybus,ych]=ybus;
g=real(ybus);b=imag(ybus);
% The given parameters and initial conditions are
p=[0;-0.96;-0.35;-0.16;0.24];
q=[0;-0.62;-0.14;-0.08;-0.35];
mv=[1.05;1;1;1;1.02];
th=[0;0;0;0;0];
del=1;indx=0;
% The Newton-Raphson iterations starts here
while del>1e-6
for i=1:5
temp=0;
for k=1:5
temp=temp+mv(i)*mv(k)*(g(i,k)-j*b(i,k))*exp(j*(th(i)th(k)));
end
pcal(i)=real(temp);qcal(i)=imag(temp);
end
% The mismatches
delp=p-pcal';
delq=q-qcal';
% The Jacobian matrix
for i=1:4
ii=i+1;
for k=1:4
kk=k+1;
j11(i,k)=mv(ii)*mv(kk)*(g(ii,kk)*sin(th(ii)-th(kk))b(ii,kk)*cos(th(ii)-th(kk)));
end
j11(i,i)=-qcal(ii)-b(ii,ii)*mv(ii)^2;
end
for i=1:4
ii=i+1;
for k=1:4
kk=k+1;
j211(i,k)=-mv(ii)*mv(kk)*(g(ii,kk)*cos(th(ii)-th(kk))b(ii,kk)*sin(th(ii)-th(kk)));
end
j211(i,i)=pcal(ii)-g(ii,ii)*mv(ii)^2;
end
j21=j211(1:3,1:4);
j12=-j211(1:4,1:3);
for i=1:3
j12(i,i)=pcal(i+1)+g(i+1,i+1)*mv(i+1)^2;
end
j22=j11(1:3,1:3);
for i=1:3
j22(i,i)=qcal(i+1)-b(i+1,i+1)*mv(i+1)^2;
end
jacob=[j11 j12;j21 j22];
delpq=[delp(2:5);delq(2:4)];
corr=inv(jacob)*delpq;
th=th+[0;corr(1:4)];
mv=mv+[0;mv(2:4).*corr(5:7);0];
del=max(abs(delpq));
indx=indx+1;
end
preal=(pcal+[0 0 0 0 0.24])*100;
preac=(qcal+[0 0 0 0 0.11])*100;
% Power flow calculations
for i=1:5
v(i)=mv(i)*exp(j*th(i));
end
for i=1:4
for k=i+1:5
if (ybus(i,k)==0)
s(i,k)=0;s(k,i)=0;
c(i,k)=0;c(k,i)=0;
q(i,k)=0;q(k,i)=0;
cur(i,k)=0;cur(k,i)=0;
else
cu=-(v(i)-v(k))*ybus(i,k);
s(i,k)=-v(i)*cu'*100;
s(k,i)=v(k)*cu'*100;
c(i,k)=100*abs(ych(i,k))*abs(v(i))^2;
c(k,i)=100*abs(ych(k,i))*abs(v(k))^2;
cur(i,k)=cu;cur(k,i)=-cur(i,k);
end
end
end
pwr=real(s);
qwr=imag(s);
q=qwr-c;
% Power loss
ilin=abs(cur);
for i=1:4
for k=i+1:5
if (ybus(i,k)==0)
pl(i,k)=0;pl(k,i)=0;
ql(i,k)=0;ql(k,i)=0;
else
z=-1/ybus(i,k);
r=real(z);
x=imag(z);
pl(i,k)=100*r*ilin(i,k)^2;pl(k,i)=pl(i,k);
ql(i,k)=100*x*ilin(i,k)^2;ql(k,i)=ql(i,k);
end
end
end
