Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- %===========================================%
- %Power Electronics in Distribution Networks %
- %Coursework Task 3 %
- %===========================================%
- %Clear Data
- clear
- clc
- %System Parameters
- V1=1.01;
- ZL12=0.05+1i*0.01;
- ZL13=0.04;
- ZL34=0.01;
- Zload3=2.21;
- Zload4=2.08+1i*0.52;
- %PG=0.6;
- PG=0;
- PSOP=0;
- QSOP2=0;
- QSOP4=0;
- %Admittance Matrix
- 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));
- G14=real(Y(1,4));
- B14=imag(Y(1,4));
- G22=real(Y(2,2));
- B22=imag(Y(2,2));
- G23=real(Y(2,3));
- B23=imag(Y(2,3));
- G24=real(Y(2,4));
- B24=imag(Y(2,4));
- 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));
- voltage_at_2=[];
- %Change QSOP2 when PG=0
- for n= 1:10000
- QSOP2=-60.01+0.01*n;
- f=@(V2,V3,V4,theta2,theta3,theta4) [
- V1*V2*(G12*cos(theta2)+B12*sin(theta2))+V2*V2*G22+V2*V3*(G23*cos(theta2-theta3)+B23*sin(theta2-theta3))+V4*V2*(G24*cos(theta2-theta4)+B24*sin(theta2-theta4))-PG+PSOP;
- V1*V3*(G13*cos(theta3)+B13*sin(theta3))+V2*V3*(G23*cos(theta3-theta2)+B23*sin(theta3-theta2))+V3*V3*G33+V4*V3*(G34*cos(theta3-theta4)+B34*sin(theta3-theta4));
- V1*V4*(G14*cos(theta4)+B14*sin(theta4))+V2*V4*(G24*cos(theta4-theta2)+B24*sin(theta4-theta2))+V4*V3*(G34*cos(theta4-theta3)+B34*sin(theta4-theta3))+V4*V4*G44-PSOP;
- V1*V2*(G12*sin(theta2)-B12*cos(theta2))-V2*V2*B22+V2*V3*(G23*sin(theta2-theta3)-B23*cos(theta2-theta3))+V2*V4*(G24*sin(theta2-theta4)-B24*cos(theta2-theta4))+QSOP2;
- V1*V3*(G13*sin(theta3)-B13*cos(theta3))+V2*V3*(G23*sin(theta3-theta2)-B23*cos(theta3-theta2))-V3*V3*B33+V3*V4*(G34*sin(theta3-theta4)-B34*cos(theta3-theta4));
- V1*V4*(G14*sin(theta4)-B14*cos(theta4))+V2*V4*(G24*sin(theta4-theta2)-B24*cos(theta4-theta2))+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);
- voltage_at_2(n)=x(1);
- check(n)=info;%%%check if the answer is reliable
- end
- %Change QSOP2 when PG=0.6
- for n= 1:10000
- PG=0.6;
- QSOP2=-60.01+0.01*n;
- f=@(V2,V3,V4,theta2,theta3,theta4) [
- V1*V2*(G12*cos(theta2)+B12*sin(theta2))+V2*V2*G22+V2*V3*(G23*cos(theta2-theta3)+B23*sin(theta2-theta3))+V4*V2*(G24*cos(theta2-theta4)+B24*sin(theta2-theta4))-PG+PSOP;
- V1*V3*(G13*cos(theta3)+B13*sin(theta3))+V2*V3*(G23*cos(theta3-theta2)+B23*sin(theta3-theta2))+V3*V3*G33+V4*V3*(G34*cos(theta3-theta4)+B34*sin(theta3-theta4));
- V1*V4*(G14*cos(theta4)+B14*sin(theta4))+V2*V4*(G24*cos(theta4-theta2)+B24*sin(theta4-theta2))+V4*V3*(G34*cos(theta4-theta3)+B34*sin(theta4-theta3))+V4*V4*G44-PSOP;
- V1*V2*(G12*sin(theta2)-B12*cos(theta2))-V2*V2*B22+V2*V3*(G23*sin(theta2-theta3)-B23*cos(theta2-theta3))+V2*V4*(G24*sin(theta2-theta4)-B24*cos(theta2-theta4))+QSOP2;
- V1*V3*(G13*sin(theta3)-B13*cos(theta3))+V2*V3*(G23*sin(theta3-theta2)-B23*cos(theta3-theta2))-V3*V3*B33+V3*V4*(G34*sin(theta3-theta4)-B34*cos(theta3-theta4));
- V1*V4*(G14*sin(theta4)-B14*cos(theta4))+V2*V4*(G24*sin(theta4-theta2)-B24*cos(theta4-theta2))+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);
- voltage_at_22(n)=x(1);
- end
- for n=1:10000
- ub(n)=1.02;
- lb(n)=0.98;
- end
- figure(1)
- index=-60:0.01:39.99;
- plot(index,voltage_at_2)
- hold on
- plot(index,voltage_at_22)
- hold on
- plot(index,ub)
- hold on
- plot(index,lb)
- hold off
- %Chnage QSOP4 with an interval of 0.01
- for n= 1:10000
- QSOP4=-60.01+0.01*n;
- f=@(V2,V3,V4,theta2,theta3,theta4) [
- V1*V2*(G12*cos(theta2)+B12*sin(theta2))+V2*V2*G22+V2*V3*(G23*cos(theta2-theta3)+B23*sin(theta2-theta3))+V4*V2*(G24*cos(theta2-theta4)+B24*sin(theta2-theta4))-PG+PSOP;
- V1*V3*(G13*cos(theta3)+B13*sin(theta3))+V2*V3*(G23*cos(theta3-theta2)+B23*sin(theta3-theta2))+V3*V3*G33+V4*V3*(G34*cos(theta3-theta4)+B34*sin(theta3-theta4));
- V1*V4*(G14*cos(theta4)+B14*sin(theta4))+V2*V4*(G24*cos(theta4-theta2)+B24*sin(theta4-theta2))+V4*V3*(G34*cos(theta4-theta3)+B34*sin(theta4-theta3))+V4*V4*G44-PSOP;
- V1*V2*(G12*sin(theta2)-B12*cos(theta2))-V2*V2*B22+V2*V3*(G23*sin(theta2-theta3)-B23*cos(theta2-theta3))+V2*V4*(G24*sin(theta2-theta4)-B24*cos(theta2-theta4))+QSOP2;
- V1*V3*(G13*sin(theta3)-B13*cos(theta3))+V2*V3*(G23*sin(theta3-theta2)-B23*cos(theta3-theta2))-V3*V3*B33+V3*V4*(G34*sin(theta3-theta4)-B34*cos(theta3-theta4));
- V1*V4*(G14*sin(theta4)-B14*cos(theta4))+V2*V4*(G24*sin(theta4-theta2)-B24*cos(theta4-theta2))+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);
- voltage_at_3(n)=x(2);
- voltage_at_4(n)=x(3);
- check(n)=info;
- end
- for n=1:10000
- ub(n)=1.02;
- lb(n)=0.98;
- end
- figure(2)
- index=-60:0.01:39.99;
- plot(index,voltage_at_3)
- hold on
- plot(index,voltage_at_4)
- hold on
- plot(index,ub)
- hold on
- plot(index,lb)
- hold off
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement