Source Code for Module pypower.newtonpf

  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 full Newton's method. 
 18  """ 
 19   
 20  import sys 
 21   
 22  from numpy import array, angle, exp, linalg, conj, r_, Inf 
 23   
 24  from scipy.sparse import hstack, vstack 
 25  from scipy.sparse.linalg import spsolve 
 26   
 27  from dSbus_dV import dSbus_dV 
 28  from ppoption import ppoption 
 29   
 30   
 31 -def newtonpf(Ybus, Sbus, V0, ref, pv, pq, ppopt=None): 
 32      """Solves the power flow using a full Newton's method. 
 33   
 34      Solves for bus voltages given the full system admittance matrix (for 
 35      all buses), the complex bus power injection vector (for all buses), 
 36      the initial vector of complex bus voltages, and column vectors with 
 37      the lists of bus indices for the swing bus, PV buses, and PQ buses, 
 38      respectively. The bus voltage vector contains the set point for 
 39      generator (including ref bus) buses, and the reference angle of the 
 40      swing bus, as well as an initial guess for remaining magnitudes and 
 41      angles. C{ppopt} is a PYPOWER options vector which can be used to 
 42      set the termination tolerance, maximum number of iterations, and 
 43      output options (see L{ppoption} for details). Uses default options if 
 44      this parameter is not given. Returns the final complex voltages, a 
 45      flag which indicates whether it converged or not, and the number of 
 46      iterations performed. 
 47   
 48      @see: L{runpf} 
 49   
 50      @author: Ray Zimmerman (PSERC Cornell) 
 51      @author: Richard Lincoln 
 52      """ 
 53      ## default arguments 
 54      if ppopt is None: 
 55          ppopt = ppoption() 
 56   
 57      ## options 
 58      tol     = ppopt['PF_TOL'] 
 59      max_it  = ppopt['PF_MAX_IT'] 
 60      verbose = ppopt['VERBOSE'] 
 61   
 62      ## initialize 
 63      converged = 0 
 64      i = 0 
 65      V = V0 
 66      Va = angle(V) 
 67      Vm = abs(V) 
 68   
 69      ## set up indexing for updating V 
 70      pvpq = r_[pv, pq] 
 71      npv = len(pv) 
 72      npq = len(pq) 
 73      j1 = 0;         j2 = npv           ## j1:j2 - V angle of pv buses 
 74      j3 = j2;        j4 = j2 + npq      ## j3:j4 - V angle of pq buses 
 75      j5 = j4;        j6 = j4 + npq      ## j5:j6 - V mag of pq buses 
 76   
 77      ## evaluate F(x0) 
 78      mis = V * conj(Ybus * V) - Sbus 
 79      F = r_[  mis[pv].real, 
 80               mis[pq].real, 
 81               mis[pq].imag  ] 
 82   
 83      ## check tolerance 
 84      normF = linalg.norm(F, Inf) 
 85      if verbose > 0: 
 86          sys.stdout.write('(Newton)\n') 
 87      if verbose > 1: 
 88          sys.stdout.write('\n it    max P & Q mismatch (p.u.)') 
 89          sys.stdout.write('\n----  ---------------------------') 
 90          sys.stdout.write('\n%3d        %10.3e' % (i, normF)) 
 91      if normF < tol: 
 92          converged = 1 
 93          if verbose > 1: 
 94              sys.stdout.write('\nConverged!\n') 
 95   
 96      ## do Newton iterations 
 97      while (not converged and i < max_it): 
 98          ## update iteration counter 
 99          i = i + 1 
100   
101          ## evaluate Jacobian 
102          dS_dVm, dS_dVa = dSbus_dV(Ybus, V) 
103   
104          J11 = dS_dVa[array([pvpq]).T, pvpq].real 
105          J12 = dS_dVm[array([pvpq]).T, pq].real 
106          J21 = dS_dVa[array([pq]).T, pvpq].imag 
107          J22 = dS_dVm[array([pq]).T, pq].imag 
108   
109          J = vstack([ 
110                  hstack([J11, J12]), 
111                  hstack([J21, J22]) 
112              ], format="csr") 
113   
114          ## compute update step 
115          dx = -1 * spsolve(J, F) 
116   
117          ## update voltage 
118          if npv: 
119              Va[pv] = Va[pv] + dx[j1:j2] 
120          if npq: 
121              Va[pq] = Va[pq] + dx[j3:j4] 
122              Vm[pq] = Vm[pq] + dx[j5:j6] 
123          V = Vm * exp(1j * Va) 
124          Vm = abs(V)            ## update Vm and Va again in case 
125          Va = angle(V)          ## we wrapped around with a negative Vm 
126   
127          ## evalute F(x) 
128          mis = V * conj(Ybus * V) - Sbus 
129          F = r_[  mis[pv].real, 
130                   mis[pq].real, 
131                   mis[pq].imag  ] 
132   
133          ## check for convergence 
134          normF = linalg.norm(F, Inf) 
135          if verbose > 1: 
136              sys.stdout.write('\n%3d        %10.3e' % (i, normF)) 
137          if normF < tol: 
138              converged = 1 
139              if verbose: 
140                  sys.stdout.write("\nNewton's method power flow converged in " 
141                                   "%d iterations.\n" % i) 
142   
143      if verbose: 
144          if not converged: 
145              sys.stdout.write("\nNewton's method power did not converge in %d " 
146                               "iterations.\n" % i) 
147   
148      return V, converged, i 
149