Coursework Task3: The SOP extracts or injects reactive power (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=2.21; % load impedances (pu)
Zload4=2.08+1i*0.52;
PG=0.6; %PV power generation (pu)
Qsop2=-2;
Qsop4=0.1;

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));

k=100
V2=zeros(1,k)
V3=zeros(1,k)
V4=zeros(1,k)
Q=zeros(1,k)
Lower_Limit=zeros(1,k)
for i=0:k
    Qsop4=Qsop4-0.1;
    disp(['Qsop4 = ',num2str(Qsop4),' pu'])


    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-Qsop2; 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-Qsop4];
    %Load flow equation [P2-PG;P3;P4;Q2-Qsop2;Q3;Q4-Qsop4]=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(10,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

    for n=2:length(x)
        x(n,:)=x(n-1,:)-((Jp(x(n-1,:))^(-1))*fp(x(n-1,:)))';
        disp(norm(fp(x(n,:))));
    end
    disp('Solution:');
    disp(x(n,:));
    V3(1,i+1)=x(n,2)
    V4(1,i+1)=x(n,3)
    Q(1,i+1)=-i/10
    Lower_Limit(1,i+1)=0.98
end
disp(V3);
disp(V4);



figure
plot(Q,V3,Q,V4,Q,Lower_Limit)

ylabel('Voltage(pu)')
xlabel('Qsop4(pu)')
legend("V3","V4","Lower limit(0.98 pu)")