Source Code for Module pypower.gausspf

  1  # Copyright (C) 1996-2011 Power System Engineering Research Center (PSERC) 
  2  # Copyright (C) 2011 Richard Lincoln 
  3  # 
  4  # PYPOWER is free software: you can redistribute it and/or modify 
  5  # it under the terms of the GNU General Public License as published 
  6  # by the Free Software Foundation, either version 3 of the License, 
  7  # or (at your option) any later version. 
  8  # 
  9  # PYPOWER is distributed in the hope that it will be useful, 
 10  # but WITHOUT ANY WARRANTY; without even the implied warranty of 
 11  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
 12  # GNU General Public License for more details. 
 13  # 
 14  # You should have received a copy of the GNU General Public License 
 15  # along with PYPOWER. If not, see <http://www.gnu.org/licenses/>. 
 16   
 17  """Solves the power flow using a Gauss-Seidel method. 
 18  """ 
 19   
 20  import sys 
 21   
 22  from numpy import linalg, conj, r_, Inf, asscalar 
 23   
 24  from ppoption import ppoption 
 25   
 26   
 27 -def gausspf(Ybus, Sbus, V0, ref, pv, pq, ppopt=None): 
 28      """Solves the power flow using a Gauss-Seidel method. 
 29   
 30      Solves for bus voltages given the full system admittance matrix (for 
 31      all buses), the complex bus power injection vector (for all buses), 
 32      the initial vector of complex bus voltages, and column vectors with 
 33      the lists of bus indices for the swing bus, PV buses, and PQ buses, 
 34      respectively. The bus voltage vector contains the set point for 
 35      generator (including ref bus) buses, and the reference angle of the 
 36      swing bus, as well as an initial guess for remaining magnitudes and 
 37      angles. C{ppopt} is a PYPOWER options vector which can be used to 
 38      set the termination tolerance, maximum number of iterations, and 
 39      output options (see C{ppoption} for details). Uses default options 
 40      if this parameter is not given. Returns the final complex voltages, 
 41      a flag which indicates whether it converged or not, and the number 
 42      of iterations performed. 
 43   
 44      @see: L{runpf} 
 45   
 46      @author: Ray Zimmerman (PSERC Cornell) 
 47      @author: Alberto Borghetti (University of Bologna, Italy) 
 48      @author: Richard Lincoln 
 49      """ 
 50      ## default arguments 
 51      if ppopt is None: 
 52          ppopt = ppoption() 
 53   
 54      ## options 
 55      tol     = ppopt['PF_TOL'] 
 56      max_it  = ppopt['PF_MAX_IT_GS'] 
 57      verbose = ppopt['VERBOSE'] 
 58   
 59      ## initialize 
 60      converged = 0 
 61      i = 0 
 62      V = V0.copy() 
 63      #Va = angle(V) 
 64      Vm = abs(V) 
 65   
 66      ## set up indexing for updating V 
 67      npv = len(pv) 
 68      npq = len(pq) 
 69      pvpq = r_[pv, pq] 
 70   
 71      ## evaluate F(x0) 
 72      mis = V * conj(Ybus * V) - Sbus 
 73      F = r_[  mis[pvpq].real, 
 74               mis[pq].imag   ] 
 75   
 76      ## check tolerance 
 77      normF = linalg.norm(F, Inf) 
 78      if verbose > 0: 
 79          sys.stdout.write('(Gauss-Seidel)\n') 
 80      if verbose > 1: 
 81          sys.stdout.write('\n it    max P & Q mismatch (p.u.)') 
 82          sys.stdout.write('\n----  ---------------------------') 
 83          sys.stdout.write('\n%3d        %10.3e' % (i, normF)) 
 84      if normF < tol: 
 85          converged = 1 
 86          if verbose > 1: 
 87              sys.stdout.write('\nConverged!\n') 
 88   
 89      ## do Gauss-Seidel iterations 
 90      while (not converged and i < max_it): 
 91          ## update iteration counter 
 92          i = i + 1 
 93   
 94          ## update voltage 
 95          ## at PQ buses 
 96          for k in pq[range(npq)]: 
 97              tmp = (conj(Sbus[k] / V[k]) - Ybus[k, :] * V) / Ybus[k, k] 
 98              V[k] = V[k] + asscalar(tmp) 
 99   
100          ## at PV buses 
101          if npv: 
102              for k in pv[range(npv)]: 
103                  tmp = (V[k] * conj(Ybus[k,:] * V)).imag 
104                  Sbus[k] = Sbus[k].real + 1j * asscalar(tmp) 
105                  tmp = (conj(Sbus[k] / V[k]) - Ybus[k, :] * V) / Ybus[k, k] 
106                  V[k] = V[k] + asscalar(tmp) 
107  #               V[k] = Vm[k] * V[k] / abs(V[k]) 
108              V[pv] = Vm[pv] * V[pv] / abs(V[pv]) 
109   
110          ## evalute F(x) 
111          mis = V * conj(Ybus * V) - Sbus 
112          F = r_[  mis[pv].real, 
113                   mis[pq].real, 
114                   mis[pq].imag  ] 
115   
116          ## check for convergence 
117          normF = linalg.norm(F, Inf) 
118          if verbose > 1: 
119              sys.stdout.write('\n%3d        %10.3e' % (i, normF)) 
120          if normF < tol: 
121              converged = 1 
122              if verbose: 
123                  sys.stdout.write('\nGauss-Seidel power flow converged in ' 
124                                   '%d iterations.\n' % i) 
125   
126      if verbose: 
127          if not converged: 
128              sys.stdout.write('Gauss-Seidel power did not converge in %d ' 
129                               'iterations.' % i) 
130   
131      return V, converged, i 
132