task2dwithPG

clear
clc

V1=1.01;
PG=0.6;
PSOP=0;
QSOP2=0;
QSOP4=0;

ZL12=0.05+i*0.01;
ZL13=0.04;
ZL34=0.01;
Zload3=2.21;
Zload4=2.08+i*0.52;

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

Y=[YL12+YL13,-YL12,-YL13,0;
-YL12,YL12,0,0;
-YL13,0,YL13+YL34+Yload3,-YL34;
0,0,-YL34,YL34+Yload4];

G12=real(Y(1,2));
B12=imag(Y(1,2));

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

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

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

for PSOP=0:0.001:1
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)); V3*V4*(G34*cos(theta4-theta3)+B34*sin(theta4-theta3))+V4*V4*G44-PSOP; V1*V2*(G12*sin(theta2)-B12*cos(theta2))-V2*V2*B22+QSOP2; V1*V3*(G13*sin(theta3)-B13*cos(theta3))-V3*V3*B33+V3*V4*(G34*sin(theta3-theta4)-B34*cos(theta3-theta4)); V3*V4*(G34*sin(theta4-theta3)-B34*cos(theta4-theta3))-V4*V4*B44-QSOP4];

fp=@(x) f(x(1),x(2),x(3),x(4),x(5),x(6));

[x, fval, info] = fsolve (fp, [1.01;1.01;1.01;0;0;0]);

disp(x);

V2=x(1);
V3=x(2);
V4=x(3);
theta2=x(4);
theta3=x(5);
theta4=x(6);

if x(1)>0.98&&x(2)>0.98&&x(3)>0.98&&x(1)<1.02&&x(2)<1.02&&x(3)<1.02
    disp(PSOP);
    break;
else
    continue;
end
end
for PSOP=0.3960:0.001:2
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)); V3*V4*(G34*cos(theta4-theta3)+B34*sin(theta4-theta3))+V4*V4*G44-PSOP; V1*V2*(G12*sin(theta2)-B12*cos(theta2))-V2*V2*B22+QSOP2; V1*V3*(G13*sin(theta3)-B13*cos(theta3))-V3*V3*B33+V3*V4*(G34*sin(theta3-theta4)-B34*cos(theta3-theta4)); V3*V4*(G34*sin(theta4-theta3)-B34*cos(theta4-theta3))-V4*V4*B44-QSOP4];

fp=@(x) f(x(1),x(2),x(3),x(4),x(5),x(6));

[x, fval, info] = fsolve (fp, [1.01;1.01;1.01;0;0;0]);

disp(x);

V2=x(1);
V3=x(2);
V4=x(3);
theta2=x(4);
theta3=x(5);
theta4=x(6);

if x(1)<0.98||x(2)<0.98||x(3)<0.98||x(1)>1.02||x(2)>1.02||x(3)>1.02
    disp(PSOP);
    break;
else
    continue;
end
end