Advertisement
Ailili1997

Coursework Task 3

Nov 11th, 2019
179
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 5.35 KB | None | 0 0
  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
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement