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

其他行业

开发平台:

Matlab

  1. %   Fast Decoupled Power Flow Solution
  2. %   Copyright (c)1998 by Hadi Saadat.
  3. ns=0; Vm=0; delta=0; yload=0; deltad=0;
  4. nbus = length(busdata(:,1));
  5. for k=1:nbus
  6. n=busdata(k,1);
  7. kb(n)=busdata(k,2); Vm(n)=busdata(k,3); delta(n)=busdata(k, 4);
  8. Pd(n)=busdata(k,5); Qd(n)=busdata(k,6); Pg(n)=busdata(k,7); Qg(n) = busdata(k,8);
  9. Qmin(n)=busdata(k, 9); Qmax(n)=busdata(k, 10);
  10. Qsh(n)=busdata(k, 11);
  11.     if Vm(n) <= 0  Vm(n) = 1.0; V(n) = 1 + j*0;
  12.     else delta(n) = pi/180*delta(n);
  13.          V(n) = Vm(n)*(cos(delta(n)) + j*sin(delta(n)));
  14.          P(n)=(Pg(n)-Pd(n))/basemva;
  15.          Q(n)=(Qg(n)-Qd(n)+ Qsh(n))/basemva;
  16.          S(n) = P(n) + j*Q(n);
  17.     end
  18. if kb(n) == 1, ns = ns+1; else, end
  19. nss(n) = ns;
  20. end
  21. Ym = abs(Ybus); t = angle(Ybus);
  22. ii=0;
  23. for ib=1:nbus
  24.      if kb(ib) == 0 | kb(ib) == 2
  25.      ii = ii+1;
  26.       jj=0;
  27.          for jb=1:nbus
  28.              if kb(jb) == 0 | kb(jb) == 2
  29.              jj = jj+1;
  30.              B1(ii,jj)=imag(Ybus(ib,jb));
  31.              else,end
  32.          end
  33.      else, end
  34. end
  35. ii=0;
  36. for ib=1:nbus
  37.      if kb(ib) == 0
  38.      ii = ii+1;
  39.       jj=0;
  40.          for jb=1:nbus
  41.              if kb(jb) == 0
  42.              jj = jj+1;
  43.              B2(ii,jj)=imag(Ybus(ib,jb));
  44.              else,end
  45.          end
  46.      else, end
  47. end
  48. B1inv=inv(B1); B2inv = inv(B2);
  49. maxerror = 1; converge = 1; 
  50. iter = 0;
  51. % Start of iterations
  52. while maxerror >= accuracy & iter <= maxiter % Test for max. power mismatch
  53. iter = iter+1;
  54. id=0; iv=0;
  55. for n=1:nbus
  56. nn=n-nss(n);
  57. J11=0;  J33=0;
  58.    for i=1:nbr
  59.      if nl(i) == n | nr(i) == n
  60.         if nl(i) == n,  l = nr(i); end
  61.         if nr(i) == n,  l = nl(i); end
  62.         J11=J11+ Vm(n)*Vm(l)*Ym(n,l)*sin(t(n,l)- delta(n) + delta(l));
  63.         J33=J33+ Vm(n)*Vm(l)*Ym(n,l)*cos(t(n,l)- delta(n) + delta(l));
  64.      else , end
  65.    end
  66.    Pk = Vm(n)^2*Ym(n,n)*cos(t(n,n))+J33;
  67.    Qk = -Vm(n)^2*Ym(n,n)*sin(t(n,n))-J11;
  68.    if kb(n) == 1 P(n)=Pk; Q(n) = Qk; end   % Swing bus P
  69.      if kb(n) == 2  Q(n)=Qk;
  70.          Qgc = Q(n)*basemva + Qd(n) - Qsh(n);
  71.          if Qmax(n) ~= 0
  72.            if iter <= 20                 % Between the 1th & 6th iterations
  73.               if iter >= 10              % the Mvar of generator buses are
  74.                 if Qgc  < Qmin(n),       % tested. If not within limits Vm(n)
  75.                 Vm(n) = Vm(n) + 0.005;   % is changed in steps of 0.05 pu to
  76.                 elseif Qgc > Qmax(n),    % bring the generator Mvar within
  77.                 Vm(n) = Vm(n) - 0.005;end % the specified limits.
  78.               else, end
  79.            else,end
  80.          else,end
  81.      end
  82.    if kb(n) ~= 1
  83.    id = id+1;
  84.      DP(id) = P(n)-Pk;
  85.      DPV(id) = (P(n)-Pk)/Vm(n);
  86.    end
  87.    if kb(n) == 0
  88.    iv=iv+1;
  89.      DQ(iv) = Q(n)-Qk;
  90.      DQV(iv) = (Q(n)-Qk)/Vm(n);
  91.    end
  92. end
  93. Dd=-B1DPV';
  94. DV=-B2DQV';
  95. id=0;iv=0;
  96.   for n=1:nbus
  97.     if kb(n) ~= 1
  98.     id = id+1;
  99.     delta(n) = delta(n)+Dd(id); end
  100.     if kb(n) == 0
  101.     iv = iv+1;
  102.     Vm(n)=Vm(n)+DV(iv); end
  103.   end
  104.     maxerror=max(max(abs(DP)),max(abs(DQ)));
  105.    if iter == maxiter & maxerror > accuracy
  106.    fprintf('nWARNING: Iterative solution did not converged after ')
  107.    fprintf('%g', iter), fprintf(' iterations.nn')
  108.    fprintf('Press Enter to terminate the iterations and print the results n')
  109.    converge = 0; pause, else, end
  110.    
  111. end
  112. if converge ~= 1
  113.    tech= ('                      ITERATIVE SOLUTION DID NOT CONVERGE'); else, 
  114.    tech=('                   Power Flow Solution by Fast Decoupled Method');
  115. end   
  116. k=0;
  117. V = Vm.*cos(delta)+j*Vm.*sin(delta);
  118. deltad=180/pi*delta;
  119. clear A; clear DC; clear DX
  120. i=sqrt(-1);
  121. for n = 1:nbus
  122.      if kb(n) == 1
  123.      S(n)=P(n)+j*Q(n);
  124.      Pg(n) = P(n)*basemva + Pd(n);
  125.      Qg(n) = Q(n)*basemva + Qd(n) - Qsh(n);
  126.      k=k+1;
  127.      Pgg(k)=Pg(n);
  128.      elseif  kb(n) ==2
  129.      S(n)=P(n)+j*Q(n);
  130.      Qg(n) = Q(n)*basemva + Qd(n) - Qsh(n);
  131.      k=k+1;
  132.      Pgg(k)=Pg(n);
  133.      end
  134. yload(n) = (Pd(n)- j*Qd(n)+j*Qsh(n))/(basemva*Vm(n)^2);
  135. end
  136. busdata(:,3)=Vm'; busdata(:,4)=deltad';
  137. Pgt = sum(Pg);  Qgt = sum(Qg); Pdt = sum(Pd); Qdt = sum(Qd); Qsht = sum(Qsh);
  138. clear Pk Qk  DP DQ J11 J33 B1 B1inv B2 B2inv DPV  DQV Dd delta ib id ii iv jb jj