BVM Engineering College (An Autonomous Institution)
Electrical Engineering Department
3EE12:POWER SYSTEMS-II LABORATORY
Practical No : Date :
AIM: To carry out load flow analysis using Newton Rapson Method
SOFTWARE: MATLAB.
Theory:
The Newton Raphson method of load flow analysis is an iterative inethod which approximates the set
of non-linear simultaneous cquations to a set of lincar simultanesus cquations using Taylor's series
expansion and the terms are limited to first order approximation. Tue load flow equations for Newton
Raphson method are non-linear equations in terms of rcal and imaginary part of bus voltages.
where, ep= Real part of Vp
fp= Imaginary part of Vp
Gpg, Bpq= Conductance and Susceptances of admittance Ypq respeclivcly.
Algorithm
Step1: Input the total number of buses. Input the details of series line impendence and line charging
admittance to calculate the Y-bus matrix.
Step2: Assume allbus voltage as 1 per unit except slack bus.
Step3: Set the iteration count as k=0 and bus count as p=1.
Step4: Calculate the real and reactive power pp and qp using the formula
P=EvpqYpq*cos(Qpq+ep-eq)
Qp=EVpqYpa*sin(qpqtep-ca)
Evalute pp*=psp-pp*
Step5: If the bus is gencrator (PV) bus, check the value of Qp*is within the limits. If it Violates
the limits, then equate the violated limit as reactivc power and treat' it as PQ bus. If limit is not
isolated then calculate,
|vpl^r=|vgpl^rspe-|vplr ;Qp*=qsp-qp*
Step6: Advance bus count by I and check if all the buses have been accounted if not go to step5.
Step7: Calculate the elements of Jacobcan matrix.
Step8: Calculate new bus voltage increment pk and fpk
Step9: Calculate new bus voltage ep*h+ cp*
Fp^k+1=fpK+fpK
Step10: Advance iteration count by I and go to step3.
Stepl1: Evaluate bus voltage and power flows through the line.
Procedure
> Enter the command window of the MATLAB.
Create a new M- file by selecting File - New -M-File.
Type and save the program in the editor Window.
Execute the program by pressing Tools - Run.
View the results.
Algorithm for the Example below:
STEP 1:Assume a flat profile l+j0 for all buses except the slack bus in the Specified voltage and
it is not modified in any iteration.
STEP 2: Assume a suitable value of [ called convergence criterion. Hence [ is A specified change
in the residue that used to compare the actual residue that is used to compare the actual
residue at the end of each iteration.
STEP 3: Set the iteration count K=0 and assunned voltage profile of the buses are denoted as
V10,V20.....Vn0.
STEP 4: Set the bus count p=1
STEP 5: Check for slack bus. If it is a slack bus then go to step 13. Otherwise Go to next step.
STEP 6: Calculate the real & reactive power of bus p using the following equation, Ppk=£k
q=l{epk(eqkGpq+fqkBpq)+fpk(fgkGpq-eqkBpq)}
Qpk-Ekq=1{fpk (eqk+Gpk+fqkBpq)+epk(qk Gpq- eqkBpq)}
STEP 7: Calculate the change in real power, APk=Pp spec -Ppk Where, Pp spec = specified real
power of bus p.
STEP 8: Check for generator bus. If it is a gencrator bus gob tc iext step othervise go to step 12.
STEP 9: Check for generator bus. If it is a generator reactivc poNcr limit Violation of generator
buses.For this compare the calculated reactive power Qwith specified limits. If the limits
are violated go to stepl1. otherwise go to next step.
STEP 10: Ifthecalculated reactive power is within the soecified limit then consider this bus as
generator bus. Now calculate the voltage residue using the equation |VPK|2 =|Vp spec2
|Vpk<2 ereh W |Vp specl= specified voltage magnitude for generation bus. Then go to
step 13.
STEP 11: Ifthe reactive power limit violated the treat this bus as a load bus. Now the specified
reactive power for this bus will correspond to Limit violated
STEP 12: Calculate the change in reactive power for load bus change in reactive power,
AQk=|4Qp spec,-Qpk
STEP 13: Repeat the stepSto 12 until all residues arc calculated for increment the bus count n.
by Ito 5 steps until the bus count isn.
STEP 14: Determine the largest of the absolute value of the residue (i.e.) Find the largest value
among APpk,AQpk or |4Vpk<2
STEP 15: Compare AE and E, if AE<E then go to step 20. If AE>E go to next Step.
STEP 16: Determine the element, the load flow equation using kth iteration Value.
STEP 17:Calculate the increment in real and reactive part of voltage Aepk and Afpk by solving the
matrix equation B=JC.
STEP 18: Calculate the new bus voltage.
STEP 19: Advance the iteration count k=k+1 and go to step 4. STEP 20: Calculate the line flows.
PROBLEM: Figure shows the one-line diagram of a three-bus power system with generators at
buses 1& 3. The magnitude of voltage at buslis adjusted to 1.05 per unit. The magnitude of voltage
at bus 3 is fixed at 1.04pu with a real power generation of 200 MW. Aload consists of 400 MW
and 250 MVAR is taken from bus2. Line impedances are marked in per unit on a 100MVA base
and the line charging susceptances are neglected.
BUS 1 0.02+ j0.04 BUS 2
400 MW
0.01+ j 0.03 0.0125+ j 0.025
250 MVAR
Slack bus
VË =1.05 Z0
BUS 3
200 MW Ivl= 1.04
i). Obtain the load flow solution by Newton-Raphson method.
ii). Write and execute a MATLAB program and also verify the output with the manual
calculation
results.
clear all;
clc;
v=[1.05;1.0;1.04);
d-[0;0;0];
ps-(-4;2.0];
qs--2.5;
n= input(Enter the number of buscs');
fprintf(Enter your choice');
p= input (1. impedance, 2. admittance');
if (p==1)
for q=1:n
for r-q+l:n
fprintf(Enter the impedance value between %d-%d'q,r);
z(q,r)=input(': );
if (z(q.r)=-0) y(q.r)-0;
else y(q.r)=inv(z(q.r));
end
y(r,q)= y(q,r);
end
end
elseif (p=-=2)
for a=1:n
for b=a+1:n
fprintf(Ernter the admittance value between %d-%d',a,b);
y(a,b)=input(':);
y(b,a)= y(a,b);
end
end
else
fprintf('enter the correct choice');
end
ybus=zeros(n,n);
for a= l:n
for b=1:n
if (a==b)
for c=l:n
ybus(a,a)= ybus(a,a)+ y(a,c);
end
else
ybus(a,b)=-y(b,.a):
end
end
end
ybus
y=abs(ybus);
tFangle(ybus);
iter-0:
pwracur-0.00025; % Power accuracy
dc=10; % Set the maximum power residualto a high value
while max(abs(dc))>pwracur
iter-iter+|
p-[v(2)*v(1)*y(2,1)*cos(t(2,1)
d(2)+d(1)+v(2y2*y(2,2)*cos(t(2,2)) +v(2)*v(3)*y(2,3) *cos((2,3)-d(2)+d(3);
v(3)*v(1)*y(3,1)*cos(t(3,1)
d(3)+d(1)+v(3)^2*y(3,3)*cos(t(3,3)}+v(3)*v(2)*y(3,2)*cos(1(3,2)-d(3)+d(2)):
q-v(2)°v(1)"y(2,1)"sin(t(2, I)-d(2)+d(1)-v2)2*y(2,2)*sin(2,2))-v(2)*V(3)*"y(2,3) "sin(t(2,3)
d(2)+d(3);
j(1,1)=v(2)*v(|)"y(2, 1)*sin(t(2,1)-d(2)+d(1) +v(2)*v(3)*y(2,3)*sin((2,3)-d(2)+d(3):
j(1,2)=-v(2)*v(3)*y(2,3)*sin(t(2,3)-d(2)+d(3));
i(1,3)=v(I)*y(2,1)*cos(t(2, I)-d(2)+d(1))+2*v(2)*y(2,2)*cos(t(2,2)) +v(3)"y(2,3)*cos((2,3)
d(2)+d(3):
j(2,I)=-v(3)*v(2)*y(3,2)*sin(3,2)-d(3)+d(2);
j(2,2)=v(3)*v(1)"y(3, I) "sin(3,2)-d(3)+d(1)+ v(3)*v(2)*y(3,2)* sin((3,2)-d(3)+d(2);
j(2,3)=v(3)*y(2,3)*cos(t(3,2)-d(3) +d(2));
j(3,I)=v(2)*v(1)"y(2, 1)*cos(I(2, 1)-d(2)+d(1)+v(2) *v(3)*y(2,3)*cos(t(2,3)-d(2)+d(3);
(3,2)=-v(2)*v(3)*y(2,3)*cos((3,2)-d(2)+d(3);
i(3,3)-v(I)"y(2,)"sin(2, 1)-d(2)+d(1)-2*v(2)*y(2,2)*sin(2,2))-v(3)*y(2,3)"sin(1(2,3)
d(2)+d(3));
dp-ps-p;
dq=qs-q;
de-[dp;dql
dx=j\dc
d(2)=d(2)+dx(1);
d(3)=d(3)+dx(2);
v(2)=v(2)+dx(3);
v,d,delta-180/pi*d;
end
pl=v(1y2*y(1,1)*cost(1,1)+v(1)*v(2)*y(1,2)*cos(t(1,2)
d(l)+d(2)) +v(1)*v(3)*y(1,3)*cos(t(1,3)-d(1)+d(3);
ql=-v(1)^2*y(1,1)*sin(t(1,1))-v(1)*v(2)*y(1,2)*sin(t(1,2)-d(1) +d(2))
v(l)*v(3) *y(1,3)*sin((1,3)-d(1)+d(3));
q3=-v(3) *v(1)*y(3, 1)*sin(t(3,1)-d(3)+d(1))-v(3)*v(2) *y(3,2)*sin((3,2)-d(3)+ d(2))
V(3^2*y(3,3)*sin(t(3,3);
Results :
