Simulation of Newton-Raphson Load Flow analysis with matlab

1 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

v 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

q 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 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];

tolerence=1;indx=0;

% The Newton-Raphson iterations starts here

while tolerence >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];

tolerence =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

pwr=real(s);

qwr=imag(s);

q=qwr-c;

disp('#########################################################################################');

disp('-----------------------------------------------------------------------------------------');

disp(' Newton Raphson Loadflow Analysis ');

disp('-----------------------------------------------------------------------------------------');

disp('| Bus | V | Angle | Specified power | ');

disp('| No | pu | Degree | MW | MVar |');

for m = 1:nbus

disp('-----------------------------------------------------------------------------------------');

fprintf('%3g', m); fprintf(' %8.4f', mv(m)); fprintf(' %8.4f', th(m)*(180/pi));

fprintf(' %8.3f',100* p(m)); fprintf(' %8.3f',100* q(m));

fprintf('\n');

end

disp('-----------------------------------------------------------------------------------------');