SHARE
TWEET

Coursework Task 3

Ailili1997 Nov 11th, 2019 96 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1.  
  2. %===========================================%
  3. %Power Electronics in Distribution Networks %
  4. %Coursework Task 3                          %
  5. %===========================================%
  6. %Clear Data
  7. clear
  8. clc
  9. %System Parameters
  10. V1=1.01;
  11. ZL12=0.05+1i*0.01;
  12. ZL13=0.04;
  13. ZL34=0.01;
  14. Zload3=2.21;
  15. Zload4=2.08+1i*0.52;
  16. %PG=0.6;
  17. PG=0;
  18. PSOP=0;
  19. QSOP2=0;
  20. QSOP4=0;
  21. %Admittance Matrix
  22. YL12=1/ZL12;
  23. YL13=1/ZL13;
  24. YL34=1/ZL34;
  25. Yload3=1/Zload3;
  26. Yload4=1/Zload4;
  27.  
  28. Y=[YL12+YL13,-YL12,-YL13,0;-YL12,YL12,0,0; -YL13,0,YL13+YL34+Yload3,-YL34;0,0,-YL34,YL34+Yload4];
  29.  
  30. G12=real(Y(1,2));
  31. B12=imag(Y(1,2));
  32.  
  33. G13=real(Y(1,3));
  34. B13=imag(Y(1,3));
  35.  
  36. G14=real(Y(1,4));
  37. B14=imag(Y(1,4));
  38.  
  39. G22=real(Y(2,2));
  40. B22=imag(Y(2,2));
  41.  
  42. G23=real(Y(2,3));
  43. B23=imag(Y(2,3));
  44.  
  45. G24=real(Y(2,4));
  46. B24=imag(Y(2,4));
  47.  
  48. G33=real(Y(3,3));
  49. B33=imag(Y(3,3));
  50.  
  51. G34=real(Y(3,4));
  52. B34=imag(Y(3,4));
  53.  
  54. G44=real(Y(4,4));
  55. B44=imag(Y(4,4));
  56.  
  57. voltage_at_2=[];
  58.  
  59. %Change QSOP2 when PG=0
  60. for n= 1:10000
  61. QSOP2=-60.01+0.01*n;
  62. f=@(V2,V3,V4,theta2,theta3,theta4) [
  63.     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;
  64.     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));
  65.     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;
  66.     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;
  67.     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));
  68.     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;];
  69.  
  70. fp=@(x) f(x(1),x(2),x(3),x(4),x(5),x(6));
  71.  
  72. [x, fval, info] = fsolve (fp, [1.01;1.01;1.01;0;0;0]);
  73. disp(x);
  74. V2=x(1);
  75. V3=x(2);
  76. V4=x(3);
  77. theta2=x(4);
  78. theta3=x(5);
  79. theta4=x(6);
  80. voltage_at_2(n)=x(1);
  81. check(n)=info;%%%check if the answer is reliable
  82. end
  83.  
  84.  
  85. %Change QSOP2 when PG=0.6
  86. for n= 1:10000
  87. PG=0.6;
  88. QSOP2=-60.01+0.01*n;
  89. f=@(V2,V3,V4,theta2,theta3,theta4) [
  90.     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;
  91.     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));
  92.     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;
  93.     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;
  94.     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));
  95.     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;];
  96.  
  97. fp=@(x) f(x(1),x(2),x(3),x(4),x(5),x(6));
  98.  
  99. [x, fval, info] = fsolve (fp, [1.01;1.01;1.01;0;0;0]);
  100.  
  101. disp(x);
  102.  
  103. V2=x(1);
  104. V3=x(2);
  105. V4=x(3);
  106. theta2=x(4);
  107. theta3=x(5);
  108. theta4=x(6);
  109. voltage_at_22(n)=x(1);
  110.  
  111. end
  112.  
  113. for n=1:10000
  114. ub(n)=1.02;
  115. lb(n)=0.98;
  116. end
  117. figure(1)
  118. index=-60:0.01:39.99;
  119. plot(index,voltage_at_2)
  120. hold on
  121. plot(index,voltage_at_22)
  122. hold on
  123. plot(index,ub)
  124. hold on
  125. plot(index,lb)
  126. hold off
  127.  
  128.  
  129. %Chnage QSOP4 with an interval of 0.01
  130. for n= 1:10000
  131. QSOP4=-60.01+0.01*n;
  132. f=@(V2,V3,V4,theta2,theta3,theta4) [
  133.     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;
  134.     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));
  135.     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;
  136.     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;
  137.     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));
  138.     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;];
  139.  
  140. fp=@(x) f(x(1),x(2),x(3),x(4),x(5),x(6));
  141.  
  142. [x, fval, info] = fsolve (fp, [1.01;1.01;1.01;0;0;0]);
  143. disp(x);
  144.  
  145. V2=x(1);
  146. V3=x(2);
  147. V4=x(3);
  148. theta2=x(4);
  149. theta3=x(5);
  150. theta4=x(6);
  151. voltage_at_3(n)=x(2);
  152. voltage_at_4(n)=x(3);
  153. check(n)=info;
  154. end
  155.  
  156. for n=1:10000
  157. ub(n)=1.02;
  158. lb(n)=0.98;
  159. end
  160. figure(2)
  161. index=-60:0.01:39.99;
  162. plot(index,voltage_at_3)
  163. hold on
  164. plot(index,voltage_at_4)
  165. hold on
  166. plot(index,ub)
  167. hold on
  168. plot(index,lb)
  169. hold off
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Top