Coursework Task1: The power-flow problem using the Newton-Raphson method (MATLAB Script)

clear
clc


V1=1.01;

ZL12=0.05+1i*0.01; %impedances between bus i,j (pu)
ZL13=0.04;
ZL34=0.01;
Zload3=1; % load impedances (pu)
Zload4=2.08+1i*0.52;
PG=0.6; %PV power generation (pu)

YL12=1/ZL12; %admittance = 1/impedance
YL13=1/ZL13;
YL34=1/ZL34;
Yload3=1/Zload3;
Yload4=1/Zload4;

Y=[YL12+YL13,-YL12,-YL13,0; %admittance matrix 4*4 (4nodes)
-YL12,YL12,0,0;
-YL13,0,YL13+YL34+Yload3,-YL34;
0,0,-YL34,YL34+Yload4];

G12=real(Y(1,2)); %Real and Imaginary values of each component in admittance matrix Y
B12=imag(Y(1,2));

G22=real(Y(2,2));
B22=imag(Y(2,2));

G13=real(Y(1,3));
B13=imag(Y(1,3));

G33=real(Y(3,3));
B33=imag(Y(3,3));

G34=real(Y(3,4));
B34=imag(Y(3,4));

G44=real(Y(4,4));
B44=imag(Y(4,4));


f=@(V2,V3,V4,theta2,theta3,theta4) [V1*V2*(G12*cos(theta2)+B12*sin(theta2))+V2*V2*G22-PG; V1*V3*(G13*cos(theta3)+B13*sin(theta3))+V3*V3*G33+V3*V4*(G34*cos(theta3-theta4)+B34*sin(theta3-theta4)); V4*V3*(G34*cos(theta4-theta3)+B34*sin(theta4-theta3))+V4*V4*G44;V2*V1*(G12*sin(theta2)-B12*cos(theta2))-V2*V2*B22; V1*V3*(G13*sin(theta3)-B13*cos(theta3))-V3*V3*B33+V4*V3*(G34*sin(theta3-theta4)-B34*cos(theta3-theta4)); V4*V3*(G34*sin(theta4-theta3)-B34*cos(theta4-theta3))-V4*V4*B44];
%Load flow equation [P2-PG;P3;P4;Q2;Q3;Q4]=0 to solve unknown parameters (V2,V3,V4,theta2,theta3,theta4)

J=@(V2,V3,V4,theta2,theta3,theta4) [2*V2*G22+V1*(G12*cos(theta2)+B12*sin(theta2)),0,0,V1*V2*(-G12*sin(theta2)+B12*cos(theta2)),0,0;0,2*V3*G33+V1*(G13*cos(theta3)+B13*sin(theta3))+V4*(G34*cos(theta3-theta4)+B34*sin(theta3-theta4)),V3*(G34*cos(theta3-theta4)+B34*sin(theta3-theta4)),0,V3*V1*(-G13*sin(theta3)+B13*cos(theta3))+V3*V4*(-G34*sin(theta3-theta4)+B34*cos(theta3-theta4)),V3*V4*(G34*sin(theta3-theta4)-B34*cos(theta3-theta4));0,V4*(G34*cos(theta4-theta3)+B34*sin(theta4-theta3)),2*V4*G44+V3*(G34*cos(theta4-theta3)+B34*sin(theta4-theta3)),0,V4*V3*(G34*sin(theta4-theta3)-B34*cos(theta4-theta3)),V4*V3*(-G34*sin(theta4-theta3)+B34*cos(theta4-theta3));-2*V2*B22+V1*(G12*sin(theta2)-B12*cos(theta2)),0,0,V2*V1*(G12*cos(theta2)+B12*sin(theta2)),0,0;0,-2*V3*B33+V1*(G13*sin(theta3)-B13*cos(theta3))+V4*(G34*sin(theta3-theta4)-B34*cos(theta3-theta4)),V3*(G34*sin(theta3-theta4)-B34*cos(theta3-theta4)),0,V3*V1*(G13*cos(theta3)+B13*sin(theta3))+V3*V4*(G34*cos(theta3-theta4)+B34*sin(theta3-theta4)),V3*V4*(-G34*cos(theta3-theta4)-B34*sin(theta3-theta4));0,V4*(G34*sin(theta4-theta3)-B34*cos(theta4-theta3)),-2*V4*B44+V3*(G34*sin(theta4-theta3)-B34*cos(theta4-theta3)),0,V4*V3*(-G34*cos(theta4-theta3)-B34*sin(theta4-theta3)),V4*V3*(G34*cos(theta4-theta3)+B34*sin(theta4-theta3))];

fp=@(x) f(x(1),x(2),x(3),x(4),x(5),x(6)); %fp is a function to input the x value in to f function
Jp=@(x) J(x(1),x(2),x(3),x(4),x(5),x(6));

x=zeros(3,6); %zero matrix x with 3*6 size
x(1,:)=[1.01;1.01;1.01;0;0;0]; %x(1,:) is an initial value for x
disp(x(1,:));
disp('Error:');
disp(norm(fp(x(1,:))));
for n=2:length(x)
    x(n,:)=x(n-1,:)-((Jp(x(n-1,:))^(-1))*fp(x(n-1,:)))';
    disp(['Iteration',num2str(n)]);
    disp(x(n,:));
    disp('Error:');
    disp(norm(fp(x(n,:))));
end
disp('Solution:');
disp(x(n,:));
disp(['V2 = ',num2str(x(n,1)),' pu']);
disp(['V3 = ',num2str(x(n,2)),' pu']);
disp(['V4 = ',num2str(x(n,3)),' pu']);
disp(['Theta2 = ',num2str(x(n,4)),' rad']);
disp(['Theta3 = ',num2str(x(n,5)),' rad']);
disp(['Theta4 = ',num2str(x(n,6)),' rad']);