NEWTON RAPHSON MATLAB CODE :-

clc
clear all;

data=[1 2 0.01008 0.05040 0.05125;

1 3 0.00744 0.03720 0.03875;
2 4 0.00744 0.03720 0.03875;
3 4 0.01272 0.06360 0.06375];

fb=data(:,1);
tb=data(:,2);
r=data(:,3);
x=data(:,4);
y2=data(:,5);
z=r+(1i*x);
y=1./z;

nbus=max(max(fb),max(tb));
nbranch=length(fb);
Y=zeros(nbus,nbus);

% Ybus OFF Diagonal Elements

for k=1:nbranch
Y(fb(k),tb(k)) = Y(fb(k),tb(k))-y(k);
Y(tb(k),fb(k)) = Y(fb(k),tb(k));
end
% Ybus Diagonal Elements
for m =1:nbus
for n =1:nbranch
if fb(n) == m
Y(m,m) = Y(m,m) + y(n)+(1i*y2(n));
elseif tb(n) == m
Y(m,m) = Y(m,m) + y(n)+(1i*y2(n));
end
end
end
disp('Y-bus matrix of given system = ')
disp(Y)

% bus no.| V | del| Pg | Qg | Pd | Qd

busdata= [1 1.00 0 0 0 50 30.99;
2 1.00 0 0 0 170 105.35;
3 1.00 0 0 0 200 123.94;
4 1.02 0 318 0 80 49.58];
BMVA=100;
busno=busdata(:,1);
v=busdata(:,2);
Vb=v;
Del=busdata(:,3);
Pg=busdata(:,4);
Qg=busdata(:,5);
Pd=busdata(:,6);
Qd=busdata(:,7);
Psch=(Pg-Pd)./BMVA;
Qsch=(Qg-Qd)./BMVA;
n=input('enter the no. of iteration = ' );

%to find the type of bus

%type-- 1-slack 2-PQ 3-PV
type=zeros(nbus,1);
type(1)=1;
for c=2:nbus
if Pg(c)==0
if Pd(c)~=0
if Qd(c)~=0
type(c)=2;
end
end
end
if Pg(c)~=0
type(c)=3;
end
end
type;

pq=find(type==2);
npq=length(pq);
pv=find(type==3);
npv=length(pv);
nou=(2*npq)+npv;

for iter=1:n
V=Vb.*exp(1i*Del);
Pcal=zeros(nbus,1);
Qcal=zeros(nbus,1);
for k=2:nbus
soi=0;
for j=1:nbus;
if j~=k
soi=soi+(V(j,1)*Y(k,j));
end
soi;
if type(k)==2 || type(k)==3
Pcal(k,1) = real(conj(V(k,1))*(soi + Y(k,k)*V(k,1)));
end
if type(k)==2
Qcal(k,1) = -imag(conj(V(k,1))*(soi + Y(k,k)*V(k,1)));
end

end
end
Pcal;
Qcal;

for x=1:npq
y=pq(x);
QschM(x,1)=Qsch(y);
QcalM(x,1)=Qcal(y);
end
delP=Psch(2:nbus)-Pcal(2:nbus);
delQ=QschM-QcalM;
MM=[delP;delQ]; %Mismatch Matrix

%to find J1
J1=zeros(nbus-1);
for a=1:(nbus-1)
m=a+1;
for b=1:(nbus-1)
l=b+1;
if m==l
J1(a,b)=-Qcal(m)-imag(conj(V(m))*V(l)*Y(m,l));
else
J1(a,b)=imag(conj(V(m)*V(l)*Y(m,l)));
end
end
end
J1;
%to find J2
j2=zeros(nbus-1,nbus-1-npv);
for a=1:(nbus-1)
m=a+1;
for b=1:(nbus-1-npv)
l=pq(b);
jjj=0;
J2(a,b)=real(V(m)*V(l)*(Y(m,l)))/abs(V(l));
if l==m
J2(a,b)=(2*J2(a,b))+real(soi);
end
end
end
J2;
%to find J3
J3=zeros(nbus-1-npv,nbus-1);
for a=1:(nbus-1-npv)
m=pq(a);
for b=1:nbus-1
l=b+1;
if l==m
for l=1:nbus
J3(a,b)=J3(a,b)+real(conj(V(m))*V(l)*Y(m,l));
end
J3(a,b)=J3(a,b)-real(conj(V(m))*V(m)*Y(m,m));
else
J3(a,b)=-real(conj(V(m))*V(l)*Y(m,l));
end
end
end
J3;

%to find J4
J4=zeros(nbus-1-npv,nbus-1-npv);
for a=1:(nbus-1-npv)
m=pq(a);
for b=1:(nbus-1-npv)
l=pq(b);
if l==m
J4(a,b)=Qcal(m)-((abs(V(m)))^2)*imag(Y(m,m));
else
J4(a,b)=(-imag(conj(V(m))*V(l)*Y(m,l)))/V(l);
%(J1(m,l))/V(l);
end

end
end
J4;
J=[J1,J2;J3,J4];
IM=(inv(J))*MM; %Improvement Matrix
delV=zeros(nbus,1);
delD=zeros(nbus,1);
delD(2:nbus)=IM(1:(nbus-1));
for t=1:npq
r=pq(t);
delV(r)=IM((nbus-1)+t);
end
Vb=Vb+delV;
Del=Del+(delD.*(180/pi));
fprintf('iteration no. = ')
disp(iter)
disp('jacobian matrix = ')
disp(J)
disp('Mismatch Matrix = ')
disp(MM)
disp('Improvement Matrix = ')
disp(IM)
disp('bus voltages = ')
disp(Vb)
disp('voltage angle del (in degree) = ')
disp(Del)
end

OUTPUT OF NR METHOD :-
Y-bus matrix of given system =
8.9852 -44.8360i -3.8156 +19.0781i -5.1696 +25.8478i 0.0000 + 0.0000i
-3.8156 +19.0781i 8.9852 -44.8360i 0.0000 + 0.0000i -5.1696 +25.8478i
-5.1696 +25.8478i 0.0000 + 0.0000i 8.1933 -40.8638i -3.0237 +15.1185i
0.0000 + 0.0000i -5.1696 +25.8478i -3.0237 +15.1185i 8.1933 -40.8638i

enter the no. of iteration = 2

iteration no. = 1

jacobian matrix =
45.4429 0 -26.3648 9.7771 0
0 41.2687 -15.4209 0 8.1933
-26.3648 -15.4209 42.5147 -5.2730 -3.0842
-9.0886 0 5.2730 44.2290 0
0 -8.2537 3.0842 0 40.4590

bus voltages =
1.0000
0.9834
0.9710
1.0200

voltage angle del (in degree) =

0
-0.9290
-1.7959
1.5162

iteration no. = 2

jacobian matrix =
-14.9906 0 -18.7126 -57.2964 0
0 -24.7429 -15.2179 0 52.4889
-18.7126 -15.2179 42.5147 -18.9957 1.2937
-29.9292 0 12.6543 101.7131 0
0 -17.8524 -5.4912 0 101.7930

bus voltages =
1.0000
0.7630
-0.0428
1.0200

voltage angle del (in degree) =

0
45.2801
-108.0489
-56.6632