pypower.runpf

53 """Runs a power flow. 54 55 Runs a power flow [full AC Newton's method by default] and optionally 56 returns the solved values in the data matrices, a flag which is C{True} if 57 the algorithm was successful in finding a solution, and the elapsed 58 time in seconds. All input arguments are optional. If C{casename} is 59 provided it specifies the name of the input data file or dict 60 containing the power flow data. The default value is 'case9'. 61 62 If the ppopt is provided it overrides the default PYPOWER options 63 vector and can be used to specify the solution algorithm and output 64 options among other things. If the 3rd argument is given the pretty 65 printed output will be appended to the file whose name is given in 66 C{fname}. If C{solvedcase} is specified the solved case will be written 67 to a case file in PYPOWER format with the specified name. If C{solvedcase} 68 ends with '.mat' it saves the case as a MAT-file otherwise it saves it 69 as a Python-file. 70 71 If the C{ENFORCE_Q_LIMS} options is set to C{True} [default is false] then 72 if any generator reactive power limit is violated after running the AC 73 power flow, the corresponding bus is converted to a PQ bus, with Qg at 74 the limit, and the case is re-run. The voltage magnitude at the bus 75 will deviate from the specified value in order to satisfy the reactive 76 power limit. If the reference bus is converted to PQ, the first 77 remaining PV bus will be used as the slack bus for the next iteration. 78 This may result in the real power output at this generator being 79 slightly off from the specified values. 80 81 Enforcing of generator Q limits inspired by contributions from Mu Lin, 82 Lincoln University, New Zealand (1/14/05). 83 84 @author: Ray Zimmerman (PSERC Cornell) 85 @author: Richard Lincoln 86 """ 87 ## default arguments 88 if casedata is None: 89 casedata = join(dirname(__file__), 'case9') 90 ppopt = ppoption(ppopt) 91 92 ## options 93 verbose = ppopt["VERBOSE"] 94 qlim = ppopt["ENFORCE_Q_LIMS"] ## enforce Q limits on gens? 95 dc = ppopt["PF_DC"] ## use DC formulation? 96 97 ## read data 98 ppc = loadcase(casedata) 99 100 ## add zero columns to branch for flows if needed 101 if ppc["branch"].shape[1] < QT: 102 ppc["branch"] = c_[ppc["branch"], 103 zeros((ppc["branch"].shape[0], 104 QT - ppc["branch"].shape[1] + 1))] 105 106 ## convert to internal indexing 107 ppc = ext2int(ppc) 108 baseMVA, bus, gen, branch = \ 109 ppc["baseMVA"], ppc["bus"], ppc["gen"], ppc["branch"] 110 111 ## get bus index lists of each type of bus 112 ref, pv, pq = bustypes(bus, gen) 113 114 ## generator info 115 on = find(gen[:, GEN_STATUS] > 0) ## which generators are on? 116 gbus = gen[on, GEN_BUS].astype(int) ## what buses are they at? 117 118 ##----- run the power flow ----- 119 t0 = time() 120 if verbose > 0: 121 v = ppver('all') 122 stdout.write('PYPOWER Version %s, %s' % (v["Version"], v["Date"])) 123 124 if dc: # DC formulation 125 if verbose: 126 stdout.write(' -- DC Power Flow\n') 127 128 ## initial state 129 Va0 = bus[:, VA] * (pi / 180) 130 131 ## build B matrices and phase shift injections 132 B, Bf, Pbusinj, Pfinj = makeBdc(baseMVA, bus, branch) 133 134 ## compute complex bus power injections [generation - load] 135 ## adjusted for phase shifters and real shunts 136 Pbus = makeSbus(baseMVA, bus, gen).real - Pbusinj - bus[:, GS] / baseMVA 137 138 ## "run" the power flow 139 Va = dcpf(B, Pbus, Va0, ref, pv, pq) 140 141 ## update data matrices with solution 142 branch[:, [QF, QT]] = zeros((branch.shape[0], 2)) 143 branch[:, PF] = (Bf * Va + Pfinj) * baseMVA 144 branch[:, PT] = -branch[:, PF] 145 bus[:, VM] = ones(bus.shape[0]) 146 bus[:, VA] = Va * (180 / pi) 147 ## update Pg for swing generator [note: other gens at ref bus are accounted for in Pbus] 148 ## Pg = Pinj + Pload + Gs 149 ## newPg = oldPg + newPinj - oldPinj 150 refgen = find(gbus == ref) ## which is[are] the reference gen[s]? 151 gen[on[refgen[0]], PG] = gen[on[refgen[0]], PG] + (B[ref, :] * Va - Pbus[ref]) * baseMVA 152 153 success = 1 154 else: ## AC formulation 155 if verbose > 0: 156 stdout.write(' -- AC Power Flow ') ## solver name and \n added later 157 158 ## initial state 159 # V0 = ones(bus.shape[0]) ## flat start 160 V0 = bus[:, VM] * exp(1j * pi/180 * bus[:, VA]) 161 V0[gbus] = gen[on, VG] / abs(V0[gbus]) * V0[gbus] 162 163 if qlim: 164 ref0 = ref ## save index and angle of 165 Varef0 = bus[ref0, VA] ## original reference bus 166 limited = [] ## list of indices of gens @ Q lims 167 fixedQg = zeros(gen.shape[0]) ## Qg of gens at Q limits 168 169 repeat = True 170 while repeat: 171 ## build admittance matrices 172 Ybus, Yf, Yt = makeYbus(baseMVA, bus, branch) 173 174 ## compute complex bus power injections [generation - load] 175 Sbus = makeSbus(baseMVA, bus, gen) 176 177 ## run the power flow 178 alg = ppopt["PF_ALG"] 179 if alg == 1: 180 V, success, _ = newtonpf(Ybus, Sbus, V0, ref, pv, pq, ppopt) 181 elif alg == 2 or alg == 3: 182 Bp, Bpp = makeB(baseMVA, bus, branch, alg) 183 V, success, _ = fdpf(Ybus, Sbus, V0, Bp, Bpp, ref, pv, pq, ppopt) 184 elif alg == 4: 185 V, success, _ = gausspf(Ybus, Sbus, V0, ref, pv, pq, ppopt) 186 else: 187 stderr.write('Only Newton''s method, fast-decoupled, and ' 188 'Gauss-Seidel power flow algorithms currently ' 189 'implemented.\n') 190 191 ## update data matrices with solution 192 bus, gen, branch = pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, pv, pq) 193 194 if qlim: ## enforce generator Q limits 195 ## find gens with violated Q constraints 196 mx = find( gen[:, GEN_STATUS] > 0 & gen[:, QG] > gen[:, QMAX] ) 197 mn = find( gen[:, GEN_STATUS] > 0 & gen[:, QG] < gen[:, QMIN] ) 198 199 if len(mx) > 0 or len(mn) > 0: ## we have some Q limit violations 200 if len(pv): 201 if verbose: 202 if len(mx) > 0: 203 print 'Gen %d [only one left] exceeds upper Q limit : INFEASIBLE PROBLEM\n' % mx 204 else: 205 print 'Gen %d [only one left] exceeds lower Q limit : INFEASIBLE PROBLEM\n' % mn 206 207 success = 0 208 break 209 210 ## one at a time? 211 if qlim == 2: ## fix largest violation, ignore the rest 212 k = argmax(r_[gen[mx, QG] - gen[mx, QMAX], 213 gen[mn, QMIN] - gen[mn, QG]]) 214 if k > len(mx): 215 mn = mn[k - len(mx)] 216 mx = [] 217 else: 218 mx = mx[k] 219 mn = [] 220 221 if verbose and len(mx) > 0: 222 print 'Gen %d at upper Q limit, converting to PQ bus\n' % mx 223 224 if verbose and len(mn) > 0: 225 print 'Gen %d at lower Q limit, converting to PQ bus\n' % mn 226 227 ## save corresponding limit values 228 fixedQg[mx] = gen[mx, QMAX] 229 fixedQg[mn] = gen[mn, QMIN] 230 mx = r_[mx, mn] 231 232 ## convert to PQ bus 233 gen[mx, QG] = fixedQg[mx] ## set Qg to binding limit 234 gen[mx, GEN_STATUS] = 0 ## temporarily turn off gen, 235 for i in range(mx): ## [one at a time, since they may be at same bus] 236 bi = gen[mx[i], GEN_BUS] ## adjust load accordingly, 237 bus[bi, [PD, QD]] = (bus[bi, [PD, QD]] - gen[mx[i], [PG, QG]]) 238 239 bus[gen[mx, GEN_BUS], BUS_TYPE] = PQ ## & set bus type to PQ 240 241 ## update bus index lists of each type of bus 242 ref_temp = ref 243 ref, pv, pq = bustypes(bus, gen) 244 if verbose and ref != ref_temp: 245 print 'Bus %d is new slack bus\n' % ref 246 247 limited = r_[limited, mx] 248 else: 249 repeat = 0 ## no more generator Q limits violated 250 else: 251 repeat = 0 ## don't enforce generator Q limits, once is enough 252 253 if qlim and len(limited) > 0: 254 ## restore injections from limited gens [those at Q limits] 255 gen[limited, QG] = fixedQg[limited] ## restore Qg value, 256 for i in range(limited): ## [one at a time, since they may be at same bus] 257 bi = gen[limited[i], GEN_BUS] ## re-adjust load, 258 bus[bi, [PD, QD]] = bus[bi, [PD, QD]] + gen[limited[i], [PG, QG]] 259 260 gen[limited, GEN_STATUS] = 1 ## and turn gen back on 261 if ref != ref0: 262 ## adjust voltage angles to make original ref bus correct 263 bus[:, VA] = bus[:, VA] - bus[ref0, VA] + Varef0 264 265 ppc["et"] = time() - t0 266 ppc["success"] = success 267 268 ##----- output results ----- 269 ## convert back to original bus numbering & print results 270 ppc["bus"], ppc["gen"], ppc["branch"] = bus, gen, branch 271 results = int2ext(ppc) 272 273 ## zero out result fields of out-of-service gens & branches 274 if len(results["order"]["gen"]["status"]["off"]) > 0: 275 results["gen"][ix_(results["order"]["gen"]["status"]["off"], [PG, QG])] = 0 276 277 if len(results["order"]["branch"]["status"]["off"]) > 0: 278 results["branch"][ix_(results["order"]["branch"]["status"]["off"], [PF, QF, PT, QT])] = 0 279 280 if fname: 281 fd = None 282 try: 283 fd = open(fname, "wb") 284 except Exception, detail: 285 stderr.write("Error opening %s: %s.\n" % (fname, detail)) 286 finally: 287 if fd is not None: 288 printpf(results, fd, ppopt) 289 fd.close() 290 291 printpf(results, stdout, ppopt) 292 293 ## save solved case 294 if solvedcase: 295 savecase(solvedcase, results) 296 297 return results, success

Source Code for Module pypower.runpf