Source Code for Module pypower.pfsoln

  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  """Updates bus, gen, branch data structures to match power flow soln. 
 18  """ 
 19   
 20  from numpy import asarray, angle, pi, conj, zeros, ones, finfo, c_, ix_ 
 21  from numpy import flatnonzero as find 
 22   
 23  from scipy.sparse import csr_matrix 
 24   
 25  from idx_bus import VM, VA, PD, QD 
 26  from idx_gen import GEN_BUS, GEN_STATUS, PG, QG, QMIN, QMAX 
 27  from idx_brch import F_BUS, T_BUS, BR_STATUS, PF, PT, QF, QT 
 28   
 29  EPS = finfo(float).eps 
 30   
 31   
 32 -def pfsoln(baseMVA, bus0, gen0, branch0, Ybus, Yf, Yt, V, ref, pv, pq): 
 33      """Updates bus, gen, branch data structures to match power flow soln. 
 34   
 35      @author: Ray Zimmerman (PSERC Cornell) 
 36      @author: Richard Lincoln 
 37      """ 
 38      ## initialize return values 
 39      bus     = bus0 
 40      gen     = gen0 
 41      branch  = branch0 
 42   
 43      ##----- update bus voltages ----- 
 44      bus[:, VM] = abs(V) 
 45      bus[:, VA] = angle(V) * 180 / pi 
 46   
 47      ##----- update Qg for all gens and Pg for swing bus ----- 
 48      ## generator info 
 49      on = find(gen[:, GEN_STATUS] > 0) ## which generators are on? 
 50      gbus = gen[on, GEN_BUS].astype(int)  ## what buses are they at? 
 51      refgen = find(gbus == ref)        ## which is(are) the reference gen(s)? 
 52   
 53      ## compute total injected bus powers 
 54      Sg = V[gbus] * conj(Ybus[gbus, :] * V) 
 55   
 56      ## update Qg for all generators 
 57      gen[:, QG] = zeros(gen.shape[0])              ## zero out all Qg 
 58      gen[on, QG] = Sg.imag * baseMVA + bus[gbus, QD]    ## inj Q + local Qd 
 59      ## ... at this point any buses with more than one generator will have 
 60      ## the total Q dispatch for the bus assigned to each generator. This 
 61      ## must be split between them. We do it first equally, then in proportion 
 62      ## to the reactive range of the generator. 
 63   
 64      if len(on) > 1: 
 65          ## build connection matrix, element i, j is 1 if gen on(i) at bus j is ON 
 66          nb = bus.shape[0] 
 67          ngon = on.shape[0] 
 68          Cg = csr_matrix((ones(ngon), (range(ngon), gbus)), (ngon, nb)) 
 69   
 70          ## divide Qg by number of generators at the bus to distribute equally 
 71          ngg = Cg * Cg.sum(0).T    ## ngon x 1, number of gens at this gen's bus 
 72          ngg = asarray(ngg).flatten()  # 1D array 
 73          gen[on, QG] = gen[on, QG] / ngg 
 74   
 75   
 76          ## divide proportionally 
 77          Cmin = csr_matrix((gen[on, QMIN], (range(ngon), gbus)), (ngon, nb)) 
 78          Cmax = csr_matrix((gen[on, QMAX], (range(ngon), gbus)), (ngon, nb)) 
 79          Qg_tot = Cg.T * gen[on, QG]## nb x 1 vector of total Qg at each bus 
 80          Qg_min = Cmin.sum(0).T       ## nb x 1 vector of min total Qg at each bus 
 81          Qg_max = Cmax.sum(0).T       ## nb x 1 vector of max total Qg at each bus 
 82          Qg_min = asarray(Qg_min).flatten()  # 1D array 
 83          Qg_max = asarray(Qg_max).flatten()  # 1D array 
 84          ## gens at buses with Qg range = 0 
 85          ig = find(Cg * Qg_min == Cg * Qg_max) 
 86          Qg_save = gen[on[ig], QG] 
 87          gen[on, QG] = gen[on, QMIN] + \ 
 88              (Cg * ((Qg_tot - Qg_min) / (Qg_max - Qg_min + EPS))) * \ 
 89                  (gen[on, QMAX] - gen[on, QMIN])    ##    ^ avoid div by 0 
 90          gen[on[ig], QG] = Qg_save  ## (terms are mult by 0 anyway) 
 91   
 92      ## update Pg for swing bus 
 93      ## inj P + local Pd 
 94      gen[on[refgen[0]], PG] = Sg[refgen[0]].real * baseMVA + bus[ref, PD] 
 95      if len(refgen) > 1:       ## more than one generator at the ref bus 
 96          ## subtract off what is generated by other gens at this bus 
 97          gen[on[refgen[0]], PG] = \ 
 98              gen[on[refgen[0]], PG] - sum(gen[on[ refgen[1:len(refgen)] ], PG]) 
 99   
100      ##----- update/compute branch power flows ----- 
101      out = find(branch[:, BR_STATUS] == 0)        ## out-of-service branches 
102      br =  find(branch[:, BR_STATUS]).astype(int) ## in-service branches 
103   
104      ## complex power at "from" bus 
105      Sf = V[ branch[br, F_BUS].astype(int) ] * conj(Yf[br, :] * V) * baseMVA 
106      ## complex power injected at "to" bus 
107      St = V[ branch[br, T_BUS].astype(int) ] * conj(Yt[br, :] * V) * baseMVA 
108      branch[ ix_(br, [PF, QF, PT, QT]) ] = c_[Sf.real, Sf.imag, St.real, St.imag] 
109      branch[ ix_(out, [PF, QF, PT, QT]) ] = zeros((len(out), 4)) 
110   
111      return bus, gen, branch 
112