obtaining p-v curves using Matlab

hi, everyone, I am trying to Obtain P-V curves using Matlab can anyone help me through it, please

[Shahabullah Amin]

clc;
clear all
syms X
z=0.1+0.5*j;
Vs=1;
A=1;
a1=real(A); a2=imag(A);
A=a1+a2*j;
B=z;
b1=real(B); b2=imag(B);
C=0;
D=A;
fi=acos(1);
K1=a1*(b2-b1*tan(fi))+a2*(b1+b2*tan(fi));
K2=a1*(b1+b2*tan(fi))+a2*(b1-b2*tan(fi));
deltarcrit=(pi/4)+0.5*atan(K2/-K1);
Vrcrit=Vs/(2*(a1*cos(deltarcrit)+a2*sin(deltarcrit)));
K3=b1*cos(deltarcrit)+b2*sin(deltarcrit);
K4=a1*cos(deltarcrit)+a2*sin(deltarcrit);
Prcrit=((Vs^2)*(2*K3*K4-(a1*b1+a2*b2)))/((b1^2+b2^2)*4*K4)
Vr=[];
for P=0.1:0.01:1

      Qr=P*tan(fi);
      P1=a1^2+a2^2;
      P2=2*P*(a1*b1+a2*b2)+2*Qr*(a1*b2+a2*b1)-Vs^2;
     P3=((b1+b2).^2)*(P^2+Qr^2);
     equation=P1*(X^2)+P2*X+P3;

      Vr=[Vr roots([P1 P2 P3])];

end

      Pr=(0.1:0.01:1);

plot(Pr,Vr(1,:))
hold on
plot(Pr,Vr(2,:))
Prcrit

this is the code which I have written but I have made some mistake in this Power and voltage curves obtaining Matlab code if anyone can help me through it, I will be thankful

[Walter Roberson]

What difficulty are you observing?

[Shahabullah Amin]

it's plotting the result in the not solving area which I have shown in the picture with the circle

[Shahabullah Amin]

Answer 1

[Walter Roberson]

clc;
clear all
syms X
z=0.1+0.5*1j;
Vs=1;
A=1;
a1=real(A); a2=imag(A);
A=a1+a2*1j;
B=z;
b1=real(B); b2=imag(B);
C=0;
D=A;
fi=acos(1);
K1=a1*(b2-b1*tan(fi))+a2*(b1+b2*tan(fi));
K2=a1*(b1+b2*tan(fi))+a2*(b1-b2*tan(fi));
deltarcrit=(pi/4)+0.5*atan(K2/-K1);
Vrcrit=Vs/(2*(a1*cos(deltarcrit)+a2*sin(deltarcrit)));
K3=b1*cos(deltarcrit)+b2*sin(deltarcrit);
K4=a1*cos(deltarcrit)+a2*sin(deltarcrit);
Prcrit=((Vs^2)*(2*K3*K4-(a1*b1+a2*b2)))/((b1^2+b2^2)*4*K4);
Vr=[];
for P=0.1:0.01:1

        Qr=P*tan(fi);
        P1=a1^2+a2^2;
        P2=2*P*(a1*b1+a2*b2)+2*Qr*(a1*b2+a2*b1)-Vs^2;
       P3=((b1+b2).^2)*(P^2+Qr^2);
       equation=P1*(X^2)+P2*X+P3;
        these_roots = roots([P1 P2 P3]);
        mask = any(imag(these_roots) ~= 0,2);
        these_roots(mask,:) = nan;
        Vr=[Vr these_roots];
  end

        Pr=(0.1:0.01:1);

plot(Pr,Vr.')
display(Prcrit)

[Walter Roberson]

The area that it does not draw is the area where the roots go complex.

If you change the P loop to

        syms P
        Qr=P*tan(fi);
        P1=a1^2+a2^2;
        P2=2*P*(a1*b1+a2*b2)+2*Qr*(a1*b2+a2*b1)-Vs^2;
       P3=((b1+b2).^2)*(P^2+Qr^2);
       equation=P1*(X^2)+P2*X+P3;
       Vr = solve(equation, X);

then because you do not change anything other than P in the loop, you can get the general form, which is

equation = (9*P^2)/25 + X^2 + X*(P/5 - 1)

and then Vr is

   1/2 - (5^(1/2)*(-(7*P - 5)*(P + 1))^(1/2))/10 - P/10
 (5^(1/2)*(-(7*P - 5)*(P + 1))^(1/2))/10 - P/10 + 1/2

That has a term

(-(7*P - 5)*(P + 1))^(1/2)

so the equation is real-valued if -(7*P - 5)*(P + 1) is positive. When P is positive (as is the case in your for loop), P+1 is always positive. So real or imaginary is going for the root is going to have a boundary when 7*P - 5 becomes 0, which is P = 5/7 which is about 0.714285714285 . Below that you have real roots; above that you have only imaginary roots.

Answer 2

[Shahabullah Amin]

the plotted signal should be like this

[Walter Roberson]

Use a higher resolution on P.

Or use the symbolic form I showed, and then

fplot(Vr, [0 0.75])

The code you posted certainly does not have real-valued solutions as far out as 2.7