Coursework Task2d: Finding the largest Psop (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; %Without PV power generation (pu)
Psop=0.8; %A Soft Open Point extracts active power from (2) and inject it into (4)

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

V2=0.96867
V3=1.0053
V4=1.0086

while (V2 < 0.98 || V2 > 1.02 || V3 < 0.98 || V3 > 1.02 || V4 < 0.98 || V4 > 1.02)
    Psop = Psop - 0.001;
    disp(['Psop = ',num2str(Psop),' pu']);
    f=@(V2,V3,V4,theta2,theta3,theta4) [V1*V2*(G12*cos(theta2)+B12*sin(theta2))+V2*V2*G22-PG+Psop;
        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-Psop;
        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+Psop;P3;P4-Psop;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

    for n=2:length(x)
        x(n,:)=x(n-1,:)-((Jp(x(n-1,:))^(-1))*fp(x(n-1,:)))';
    end
    disp('Solution:');
    disp(x(n,:));
    V2=x(n,1)
    V3=x(n,2)
    V4=x(n,3)
end
disp(['The largest Psop = ',num2str(Psop),' pu']);