Solution of Load Flow by Newton Raphson method without Q-limit Violation in MATLAB | PSS LAB | EEE

 PROBLEM:

The load flow for sample power system are given below 1.04p.u and the minimum reactive power units of generator bus 2 are 0.35 and zero per unit reactive. Determine the set of flow equation at the end of first iteration by using Newton Raphson method.

BUS CODE	IMPEDANCE	LINE CHARGING ADMITTANCE
1-2	0.08+0.24j	0
1-3	0.02+0.06j	0
2-3	0.06+0.18j	0

 

BUS CODE	VOLTAGE	GENERATOR		LOAD
		P	Q	P	Q
1	1.06 ̷ ̲00	0	0	0	0
2	1.04 ̷ ̲00	0.2	0	0	0
3	1.0 ̷ ̲00	0	0	0.6	0.25

Theoretical Calculation:

MATLAB Program:

disp('          Newton - Raphson Method         ');

clear;

clear all;

n=input('No of buses:');

l=input('No of lines:');

s=input('Enter Impedance 1 or Admittance 2:');

for i=1:l

    a=input('Starting bus:');

    b=input('Ending bus:');

    t=input('Impedance/Admittance value:');

if s==1

    y(a,b)=1/t;

else

        y(a,b)=t;

end

y(b,a)=y(a,b);

end

ybus=zeros(n,n);

for i=1:n

for j=1:n

if i==j

for k=1:n

   ybus(i,j)=ybus(i,j)+y(i,k);

end

else

   ybus(i,j)=-y(i,j);

end

   ybus(j,i)=ybus(i,j);

end

end

disp(ybus);

y1=zeros(n,n);

y2=zeros(n,n);

y1=abs(ybus);

y2=angle(ybus);

disp(y1);

disp(y2);

v=zeros(1,n);

d=zeros(1,n);

pspec=zeros(1,n);

qspec=zeros(1,n);

for i=1:n

    v(i)=input('Enter V value:');

    d(i)=input('Enter deg value:');

    pspec(i)=input('Enter pspec value:');

    qspec(i)=input('Enter qspec value:');

end

disp(v);

disp(d);

pcal=zeros(n,1);

qcal=zeros(n,1);

for i=2:n

for j=1:n

       pcal(i)=pcal(i)+(v(i)*y1(i,j)*v(j)*cos(y2(i,j)+d(j)-d(i)));

       qcal(i)=qcal(i)+((-1)*(v(i)*y1(i,j)*v(j)*sin(y2(i,j)+d(j)-d(i))));

end

end

disp(pcal);

disp(qcal);

dp=zeros(n,1);

dq=zeros(n,1);

for i=2:n

    dp(i)=pspec(i)-pcal(i);

    dq(i)=qspec(i)-qcal(i);

end

disp(dp);

disp(dq);

 

D=zeros(n,1);

D=[dp(2);dp(3);dq(3)];

disp(D);

J1=zeros(n,n);

for i=2:n

for j=1:n

if j~=i

J1(i,i)=(v(i)*y1(i,j)*v(j)*sin((y2(i,j))+d(j)-d(i)));

for j=2:n

J1(i,j)=-(v(i)*y1(i,j)*v(j)*sin((y2(i,j))+d(j)-d(i)));

end

end

end

end

 disp(J1);

 J1(1,:)=[];

J1(:,1)=[];

disp(J1);

J2=zeros(n,n);

for i=3:n

for j=1:n

if j~=i

        a1=0;

        a2=0;

           a1=a1+(v(j)*y1(i,j)*cos((y2(i,j))+d(j)-d(i)));

           a2=a2+(2*v(i)^2)*y1(i,j)*cos(y2(i,j)+d(j)-d(i));

          J2(i,i)=a1+a2;

for i=2:n

for j=3:n

   J2(i,j)=(v(i)*v(j)*y1(i,j)*cos((y2(i,j))+d(j)-d(i)));

end

end

end

end

end

disp(J2);

J2(:,1)=[];

J2(1,:)=[];

J2(:,1)=[];

disp(J2);

J3=zeros(n,n);

for i=2:n

for j=1:n

if j~=i

     J3(i,i)=(v(i)*y1(i,j)*v(j)*cos(y2(i,j)+d(j)-d(i)));

end

for j=2:n

J3(i,j)=-(v(i)*y1(i,j)*v(j)*cos(y2(i,j)+d(j)-d(i)));

end

end

end

J3(1,:)=[];

J3(:,1)=[];

J3(1,:)=[];

disp(J3);

J4=zeros(n,n);

for i=3:n

for j=3:n

        b1=0;

b1=b1-(1*(2*(v(i)^2)*y1(i,i)*sin(y2(i,i)+d(j)-d(i))));

        b2=0;

for i=3:n

for j=1:n

if j~=i                      

b2=b2-(1*v(i)*v(j)*y1(i,j)*sin(y2(i,j))+d(j)-d(i));

end

 J4(i,i)=b1+b2;

end

end

end

end

disp(J4);

J4(1,:)=[];

J4(:,1)=[];

J4(1,:)=[];

J4(:,1)=[];

disp(J4);

J=zeros(n,n);

J=[J1 J2;J3 J4];

disp(J);

R=zeros(n,n);

R=inv(J);

disp(R);

E=zeros(n,1);

E=R*D;

disp(E);

dln=zeros();

for i=2:n

    dln(1)=0;

    dln(i)=d(i)+E(i-1,1);

end

disp(dln);

vnew=zeros();

for i=1:n

    vnew(i)=v(i);

end

for i=3:n

vnew(i)=v(i)+v(i)*E(i,1);

end

disp(vnew);

vnew1=zeros();

for i=1:n

vnew1(i)=vnew(i)*(cos(dln(i))+(1j*sin(dln(i))));

end

VoltageNew = ['Vnew = ', num2str(vnew1)];

disp(VoltageNew);

MATLAB Output:

YouTube Video of My Channel on the Problem you can use this also

  

    ( If you have any problem in doing this you can comment below, I will help as soon as possible)

💗  Thanks For Visiting  💗