Eacpower.m
上传用户:eighthdate
上传日期:2014-05-24
资源大小:270k
文件大小:2k
源码类别:

其他行业

开发平台:

Matlab

  1. function eacpower(P0, E, V, X)
  2. % This program obtains the power angle curve for a one-machine system
  3. % during normal operation. Using equal area criterion the maximum input
  4. % power that can be suddenly applied for the machine to remain critically
  5. % stable is obtained.
  6. %
  7. % Copyright (c) 1998 H.  Saadat
  8. if exist('P0')~=1
  9. P0 = input('Generator initial power in p.u. P0 = '); else, end
  10. if exist('E')~=1
  11. E = input('Generator e.m.f. in p.u. E = '); else, end
  12. if exist('V')~=1
  13. V = input('Infinite bus-bar voltage in p.u. V = '); else, end
  14. if exist('X')~=1
  15. X = input('Reactance between internal emf and infinite bus in p.u. X = '); else, end
  16. Pemax= E*V/X;
  17. if P0 >= Pemax
  18.   fprintf('nP0 must be less than the peak electrical power Pemax = %5.3f p.u. n', Pemax)
  19.   fprintf('Try again. nn')
  20. return, end
  21. d0=asin(P0/Pemax);
  22. delta = 0:.01:pi;
  23. Pe = Pemax*sin(delta);
  24. dmax=pi;
  25. Ddmax=1;
  26. while abs(Ddmax) > 0.00001
  27. Df = cos(d0) - (sin(dmax)*(dmax-d0)+cos(dmax));
  28. J=cos(dmax)*(dmax-d0);
  29. Ddmax=Df/J;
  30. dmax=dmax+Ddmax;
  31. end
  32. dc=pi-dmax;
  33. Pm2=Pemax*sin(dc);
  34. Pmx =[0  pi-d0]*180/pi; Pmy=[P0 P0];
  35. Pm2x=[0  dmax]*180/pi; Pm2y=[Pm2 Pm2];
  36. x0=[d0 d0]*180/pi; y0=[0 Pm2]; xc=[dc dc]*180/pi; yc=[0 Pemax*sin(dc)];
  37. xm=[dmax dmax]*180/pi; ym=[0 Pemax*sin(dmax)];
  38. d0=d0*180/pi; dmax=dmax*180/pi; dc=dc*180/pi;
  39. x=(d0:.1:dc);
  40. y=Pemax*sin(x*pi/180);
  41. %y1=Pe2max*sin(d0*pi/180);
  42. %y2=Pe2max*sin(dc*pi/180);
  43. x=[d0 x dc];
  44. y=[Pm2 y Pm2];
  45. xx=dc:.1:dmax;
  46. h=Pemax*sin(xx*pi/180);
  47. xx=[dc xx  dmax];
  48. hh=[Pm2 h  Pm2];
  49. delta=delta*180/pi;
  50. %clc
  51. fprintf('nInitial power                      =%7.3f p.u.n', P0)
  52. fprintf('Initial power angle                =%7.3f degrees n', d0)
  53. fprintf('Sudden additional power            =%7.3f p.u.n', Pm2-P0)
  54. fprintf('Total power for critical stability =%7.3f p.u.n', Pm2)
  55. fprintf('Maximum angle swing                =%7.3f degrees n', dmax)
  56. fprintf('New operating angle                =%7.3f degrees nnn', dc)
  57. fill(x,y,'m')
  58. hold;
  59. fill(xx,hh,'c')
  60. plot(delta, Pe,'-', Pmx, Pmy,'g', Pm2x,Pm2y,'g', x0,y0,'c', xc,yc, xm,ym,'r'), grid
  61. Title('Equal-area criterion applied to the sudden change in power')
  62. xlabel('Power angle, degree'), ylabel(' Power, per unit')
  63. axis([0 180  0 1.1*Pemax])
  64. hold off;