Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- clear
- clc
- % Fixed voltage magnitude at the slack bus
- V1=1.01;
- Pg=0;
- PSOP=0;
- % Load parameters of the two feeder system (p.u.)
- ZL12=0.05+i*0.01;
- ZL13=0.04;
- ZL34=0.01;
- Zload3=2.21;
- Zload4=2.08+i*0.52;
- % Admittance calculation
- YL12=1/ZL12;
- YL13=1/ZL13;
- YL34=1/ZL34;
- Yload3=1/Zload3;
- Yload4=1/Zload4;
- % Admittane matrix
- Y=[YL12+YL13,-YL12,-YL13,0;
- -YL12,YL12,0,0;
- -YL13,0,YL13+YL34+Yload3,-YL34;
- 0,0,-YL34,YL34+Yload4];
- % Obtaining real part of the complex number for Y matrix
- G11=real(Y(1,1));
- % Obtaining imaginary part of the complex number for Y matrix
- B11=imag(Y(1,1));
- 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));
- PSOPArray = [];
- for PSOP = 0:0.01:1
- % PSOP has been added for Q2a to extract 0.8 pu of active power from (2) to (4)
- % PSOP not added for Q1a
- % 6 unknowns variables - V2, V3, V4, theta2, theta3 and theta 4 are
- % calculated below by creating the function
- f=@(V2,V3,V4,theta2,theta3,theta4) [V1*V2*(G12*cos(theta2)+B12*sin(theta2))+V2*V2*G22 + PSOP;
- V1*V3*(G13*cos(theta3)+B13*sin(theta3))+V3*V3*G33+V3*V4*(G34*cos(theta3-theta4)+B34*sin(theta3-theta4));
- % PSOP not added for Q1a
- % PSOP has been subtracted for Q2a to extract 0.8 pu of active power from (2) to (4)
- V4*V4*G44+V3*V4*(G34*cos(theta4-theta3)+B34*sin(theta4-theta3))- PSOP;
- V1*V2*(G12*sin(theta2)-B12*cos(theta2))-V2*V2*B22;
- V1*V3*(G13*sin(theta3)-B13*cos(theta3))-V3*V3*B33+V3*V4*(G34*sin(theta3-theta4)-B34*cos(theta3-theta4));
- -V4*V4*B44+V3*V4*(G34*sin(theta4-theta3)-B34*cos(theta4-theta3))];
- % derivative of the defined function to build Jacabian, J.
- J=@(V2,V3,V4,theta2,theta3,theta4) [V1*(G12*cos(theta2)+B12*sin(theta2))+2*V2*G22, 0, 0, V1*V2*(-G12*sin(theta2)+B12*cos(theta2)), 0, 0;
- 0, V1*(G13*cos(theta3)+B13*sin(theta3))+2*V3*G33+V4*(G34*cos(theta3-theta4)+B34*sin(theta3-theta4)), V3*(G34*cos(theta3-theta4)+B34*sin(theta3-theta4)), 0, V1*V3*(-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)), V3*(G34*cos(theta4-theta3)+B34*sin(theta4-theta3))+2*V4*G44, 0, V3*V4*(G34*sin(theta4-theta3)-B34*cos(theta4-theta3)), V3*V4*(-G34*sin(theta4-theta3)+B34*cos(theta4-theta3));
- V1*(G12*sin(theta2)-B12*cos(theta2))-2*V2*B22, 0, 0, V1*V2*(G12*cos(theta2)+B12*sin(theta2)),0,0;
- 0, V1*(G13*sin(theta3)-B13*cos(theta3))-2*V3*B33+V4*(G34*sin(theta3-theta4)-B34*cos(theta3-theta4)), V3*(G34*sin(theta3-theta4)-B34*cos(theta3-theta4)), 0, V1*V3*(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, V3*V4*(-G34*cos(theta4-theta3)-B34*sin(theta4-theta3)), V3*V4*(G34*cos(theta4-theta3)+B34*sin(theta4-theta3))];
- fp=@(x) f(x(1),x(2),x(3),x(4),x(5),x(6));
- Jp=@(x) J(x(1),x(2),x(3),x(4),x(5),x(6));
- % Newton-Raphson
- x=zeros(6,6);
- x(1,:)=[1.01,1.01,1.01,0,0,0];
- 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(norm(fp(x(n,:))));
- end
- disp('Solution:');
- disp(x(n,:));
- PSOPArray = [PSOPArray [x(n,1); x(n,2) ; x(n,3)]];
- end
- clf
- hold on
- plot(0:0.01:1,PSOPArray')
- figure(1)
- hold on
- % plot graph of Voltages against PSOP for Q2d
- disp('Intersection 1:');
- P = InterX([0:0.01:1;PSOPArray(1,:)],[0:0.01:1;PSOPArray(2,:)])
- plot(P(1),P(2),'ko');
- disp('Intersection 2:');
- P = InterX([0:0.01:1;PSOPArray(1,:)],[0:0.01:1;PSOPArray(3,:)])
- plot(P(1),P(2),'bo');
- disp('Intersection 3:');
- P = InterX([0:0.01:1;PSOPArray(2,:)],[0:0.01:1;PSOPArray(3,:)])
- plot(P(1),P(2),'ro');
- hold off
- xlabel ('PSOP (p.u.)')
- ylabel ('Voltage (p.u.)')
- title ('Voltages against PSOP range')
- legend("V2","V3","V4")
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement