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
[Link] 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 :