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

其他行业

开发平台:

Matlab

  1. %   Power flow solution by Gauss-Seidel method
  2. %   Copyright (c) 1998 by  H. Saadat
  3. 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. DV(n)=0;
  19. end
  20. num = 0; AcurBus = 0; converge = 1;
  21. Vc = zeros(nbus,1)+j*zeros(nbus,1); Sc = zeros(nbus,1)+j*zeros(nbus,1);
  22. while exist('accel')~=1
  23.    accel = 1.3;
  24. end
  25. while exist('accuracy')~=1
  26.    accuracy = 0.001;
  27. end
  28. while exist('basemva')~=1
  29.    basemva= 100;
  30. end
  31. while exist('maxiter')~=1
  32.    maxiter = 100;
  33. end
  34. iter=0;
  35. maxerror=10;
  36. while maxerror >= accuracy & iter <= maxiter
  37. iter=iter+1;
  38.   for n = 1:nbus;
  39.   YV = 0+j*0;
  40.     for L = 1:nbr;
  41.             if nl(L) == n, k=nr(L);
  42.             YV = YV + Ybus(n,k)*V(k);
  43.             elseif nr(L) == n, k=nl(L);
  44.             YV = YV + Ybus(n,k)*V(k);
  45.             end
  46.     end
  47.        Sc = conj(V(n))*(Ybus(n,n)*V(n) + YV) ;
  48.        Sc = conj(Sc);
  49.        DP(n) = P(n) - real(Sc);
  50.        DQ(n) = Q(n) - imag(Sc);
  51.          if kb(n) == 1
  52.          S(n) =Sc; P(n) = real(Sc); Q(n) = imag(Sc); DP(n) =0; DQ(n)=0;
  53.          Vc(n) = V(n);
  54.          elseif kb(n) == 2
  55.          Q(n) = imag(Sc); S(n) = P(n) + j*Q(n);
  56.            if Qmax(n) ~= 0
  57.              Qgc = Q(n)*basemva + Qd(n) - Qsh(n);
  58.              if abs(DQ(n)) <= .005 & iter >= 10 % After 10 iterations
  59.                if DV(n) <= 0.045                % the Mvar of generator buses are
  60.                   if Qgc < Qmin(n),             % tested. If not within limits Vm(n)
  61.                   Vm(n) = Vm(n) + 0.005;        % is changed in steps of 0.005 pu
  62.                   DV(n) = DV(n)+.005;           % up to .05  pu in order to bring
  63.                   elseif Qgc > Qmax(n),         % the generator Mvar within the
  64.                   Vm(n) = Vm(n) - 0.005;        % specified limits.
  65.                   DV(n)=DV(n)+.005; end
  66.                else, end
  67.              else,end
  68.            else,end
  69.          end
  70.        if kb(n) ~= 1
  71.        Vc(n) = (conj(S(n))/conj(V(n)) - YV )/ Ybus(n,n);
  72.        else, end
  73.           if kb(n) == 0
  74.           V(n) = V(n) + accel*(Vc(n)-V(n));
  75.           elseif kb(n) == 2
  76.           VcI = imag(Vc(n));
  77.           VcR = sqrt(Vm(n)^2 - VcI^2);
  78.           Vc(n) = VcR + j*VcI;
  79.            V(n) = V(n) + accel*(Vc(n) -V(n));
  80.           end
  81.    end
  82.   maxerror=max( max(abs(real(DP))), max(abs(imag(DQ))) );
  83.    if iter == maxiter & maxerror > accuracy
  84.    fprintf('nWARNING: Iterative solution did not converged after ')
  85.    fprintf('%g', iter), fprintf(' iterations.nn')
  86.    fprintf('Press Enter to terminate the iterations and print the results n')
  87.    converge = 0; pause, else, end
  88.    
  89. end
  90. if converge ~= 1
  91.    tech= ('                      ITERATIVE SOLUTION DID NOT CONVERGE'); else, 
  92.    tech=('                   Power Flow Solution by Gauss-Seidel Method');
  93. end   
  94. k=0;
  95. for n = 1:nbus
  96.   Vm(n) = abs(V(n)); deltad(n) = angle(V(n))*180/pi;
  97.      if kb(n) == 1
  98.      S(n)=P(n)+j*Q(n);
  99.      Pg(n) = P(n)*basemva + Pd(n);
  100.      Qg(n) = Q(n)*basemva + Qd(n) - Qsh(n);
  101.      k=k+1;
  102.      Pgg(k)=Pg(n);
  103.      elseif  kb(n) ==2
  104.      k=k+1;
  105.      Pgg(k)=Pg(n);
  106.      S(n)=P(n)+j*Q(n);
  107.      Qg(n) = Q(n)*basemva + Qd(n) - Qsh(n);
  108.      end
  109. yload(n) = (Pd(n)- j*Qd(n)+j*Qsh(n))/(basemva*Vm(n)^2);
  110. end
  111. Pgt = sum(Pg);  Qgt = sum(Qg); Pdt = sum(Pd); Qdt = sum(Qd); Qsht = sum(Qsh);
  112. busdata(:,3)=Vm'; busdata(:,4)=deltad';
  113. clear  AcurBus  DP  DQ  DV  L Sc Vc VcI VcR YV converge delta