pylon/ac_pf.py

   1 #------------------------------------------------------------------------------
   2 # Copyright (C) 1996-2010 Power System Engineering Research Center (PSERC)
   3 # Copyright (C) 2007-2010 Richard Lincoln
   4 #
   5 # Licensed under the Apache License, Version 2.0 (the "License");
   6 # you may not use this file except in compliance with the License.
   7 # You may obtain a copy of the License at
   8 #
   9 #     http://www.apache.org/licenses/LICENSE-2.0
  10 #
  11 # Unless required by applicable law or agreed to in writing, software
  12 # distributed under the License is distributed on an "AS IS" BASIS,
  13 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14 # See the License for the specific language governing permissions and
  15 # limitations under the License.
  16 #------------------------------------------------------------------------------
  17
  18 """ Defines solvers for AC power flow.
  19
  20 Based on runpf.m from MATPOWER by Ray Zimmerman, developed at PSERC Cornell.
  21 See U{http://www.pserc.cornell.edu/matpower/} for more information.
  22 """
  23
  24 #------------------------------------------------------------------------------
  25 #  Imports:
  26 #------------------------------------------------------------------------------
  27
  28 import logging
  29 from time import time
  30
  31 from numpy import array, angle, pi, exp, linalg, multiply, conj, r_, Inf
  32
  33 from scipy.sparse import hstack, vstack
  34 from scipy.sparse.linalg import spsolve, splu
  35
  36 from pylon.case import PQ, PV, REFERENCE
  37
  38 #------------------------------------------------------------------------------
  39 #  Logging:
  40 #------------------------------------------------------------------------------
  41
  42 logger = logging.getLogger(__name__)
  43
  44 #------------------------------------------------------------------------------
  45 #  Constants:
  46 #------------------------------------------------------------------------------
  47
  48 BX = "BX"
  49 XB = "XB"
  50
  51 #------------------------------------------------------------------------------
  52 #  Exceptions:
  53 #------------------------------------------------------------------------------
  54
  55 class SlackBusError(Exception):
  56     """ No single slack bus error. """
  57
  58 #------------------------------------------------------------------------------
  59 #  "_ACPF" class:
  60 #------------------------------------------------------------------------------
  61
  62 class _ACPF(object):
  63     """ Defines a base class for AC power flow solvers.
  64
  65     Based on runpf.m from MATPOWER by Ray Zimmerman, developed at PSERC
  66     Cornell. See U{http://www.pserc.cornell.edu/matpower/} for more info.
  67     """
  68
  69     #--------------------------------------------------------------------------
  70     #  "object" interface:
  71     #--------------------------------------------------------------------------
  72
  73     def __init__(self, case, qlimit=False, tolerance=1e-08, iter_max=10,
  74                  verbose=True):
  75         #: Solved case.
  76         self.case = case
  77
  78         #: Enforce Q limits on generators.
  79         self.qlimit = qlimit
  80
  81         #: Convergence tolerance.
  82         self.tolerance = tolerance
  83
  84         #: Maximum number of iterations.
  85         self.iter_max = iter_max
  86
  87         #: Print progress information.
  88         self.verbose = verbose
  89
  90     #--------------------------------------------------------------------------
  91     #  "_ACPF" interface:
  92     #--------------------------------------------------------------------------
  93
  94     def solve(self):
  95         """ Runs a power flow
  96
  97         @rtype: dict
  98         @return: Solution dictionary with the following keys:
  99                    - C{V} - final complex voltages
 100                    - C{converged} - boolean value indicating if the solver
 101                      converged or not
 102                    - C{iterations} - the number of iterations performed
 103         """
 104         # Zero result attributes.
 105         self.case.reset()
 106
 107         # Retrieve the contents of the case.
 108         b, l, g, _, _, _, _ = self._unpack_case(self.case)
 109
 110         # Update bus indexes.
 111         self.case.index_buses(b)
 112
 113         # Index buses accoding to type.
 114 #        try:
 115 #            _, pq, pv, pvpq = self._index_buses(b)
 116 #        except SlackBusError:
 117 #            logger.error("Swing bus required for DCPF.")
 118 #            return {"converged": False}
 119
 120         refs, pq, pv, pvpq = self._index_buses(b)
 121         if len(refs) != 1:
 122             logger.error("Swing bus required for DCPF.")
 123             return {"converged": False}
 124
 125         # Start the clock.
 126         t0 = time()
 127
 128         # Build the vector of initial complex bus voltages.
 129         V0 = self._initial_voltage(b, g)
 130
 131         # Save index and angle of original reference bus.
 132 #        if self.qlimit:
 133 #            ref0 = ref
 134 #            Varef0 = b[ref0].Va
 135 #            # List of buses at Q limits.
 136 #            limits = []
 137 #            # Qg of generators at Q limits.
 138 #            fixedQg = matrix(0.0, (g.size[0], 1))
 139
 140         repeat = True
 141         while repeat:
 142             # Build admittance matrices.
 143             Ybus, Yf, Yt = self.case.getYbus(b, l)
 144
 145             # Compute complex bus power injections (generation - load).
 146             Sbus = self.case.getSbus(b)
 147
 148             # Run the power flow.
 149             V, converged, i = self._run_power_flow(Ybus, Sbus, V0, pv, pq, pvpq)
 150
 151             # Update case with solution.
 152             self.case.pf_solution(Ybus, Yf, Yt, V)
 153
 154             # Enforce generator Q limits.
 155             if self.qlimit:
 156                 raise NotImplementedError
 157             else:
 158                 repeat = False
 159
 160         elapsed = time() - t0
 161
 162         if converged and self.verbose:
 163             logger.info("AC power flow converged in %.3fs" % elapsed)
 164
 165         return {"converged": converged, "elapsed": elapsed, "iterations": i,
 166                 "V":V}
 167
 168
 169     def _unpack_case(self, case):
 170         """ Returns the contents of the case to be used in the OPF.
 171         """
 172         base_mva = case.base_mva
 173         b = case.connected_buses
 174         l = case.online_branches
 175         g = case.online_generators
 176         nb = len(b)
 177         nl = len(l)
 178         ng = len(g)
 179
 180         return b, l, g, nb, nl, ng, base_mva
 181
 182
 183     def _index_buses(self, buses):
 184         """ Set up indexing for updating v.
 185         """
 186         refs = [bus._i for bus in buses if bus.type == REFERENCE]
 187 #        if len(refs) != 1:
 188 #            raise SlackBusError
 189         pv = [bus._i for bus in buses if bus.type == PV]
 190         pq = [bus._i for bus in buses if bus.type == PQ]
 191         pvpq = pv + pq
 192
 193         return refs, pq, pv, pvpq
 194
 195
 196     def _initial_voltage(self, buses, generators):
 197         """ Returns the initial vector of complex bus voltages.
 198
 199         The bus voltage vector contains the set point for generator
 200         (including ref bus) buses, and the reference angle of the swing
 201         bus, as well as an initial guess for remaining magnitudes and
 202         angles.
 203         """
 204         Vm = array([bus.v_magnitude for bus in buses])
 205
 206         # Initial bus voltage angles in radians.
 207         Va = array([bus.v_angle * (pi / 180.0) for bus in buses])
 208
 209         V = Vm * exp(1j * Va)
 210
 211         # Get generator set points.
 212         for g in generators:
 213             i = g.bus._i
 214             V[i] = g.v_magnitude / abs(V[i]) * V[i]
 215
 216         return V
 217
 218
 219     def _run_power_flow(self, Ybus, Sbus, V0):
 220         """ Override this method in subclasses.
 221         """
 222         raise NotImplementedError
 223
 224 #------------------------------------------------------------------------------
 225 #  "NewtonPF" class:
 226 #------------------------------------------------------------------------------
 227
 228 class NewtonPF(_ACPF):
 229     """ Solves the power flow using full Newton's method.
 230
 231     Based on newtonpf.m from MATPOWER by Ray Zimmerman, developed at PSERC
 232     Cornell. See U{http://www.pserc.cornell.edu/matpower/} for more info.
 233     """
 234
 235     def _run_power_flow(self, Ybus, Sbus, V, pv, pq, pvpq, **kw_args):
 236         """ Solves the power flow using a full Newton's method.
 237         """
 238         Va = angle(V)
 239         Vm = abs(V)
 240
 241         # Initial evaluation of F(x0)...
 242         F = self._evaluate_function(Ybus, V, Sbus, pv, pq)
 243         # ...and convergency check.
 244         converged = self._check_convergence(F)
 245
 246         # Perform Newton iterations.
 247         i = 0
 248         while (not converged) and (i < self.iter_max):
 249             V, Vm, Va = self._one_iteration(F, Ybus, V, Vm, Va, pv, pq, pvpq)
 250             F = self._evaluate_function(Ybus, V, Sbus, pv, pq)
 251             converged = self._check_convergence(F)
 252             i += 1
 253
 254         if converged:
 255             if self.verbose:
 256                 logger.info("Newton's method power flow converged in %d "
 257                             "iterations." % i)
 258         else:
 259             logger.error("Newton's method power flow did not converge in %d "
 260                          "iterations." % i)
 261
 262         return V, converged, i
 263
 264
 265     def _one_iteration(self, F, Ybus, V, Vm, Va, pv, pq, pvpq):
 266         """ Performs one Newton iteration.
 267         """
 268         J = self._build_jacobian(Ybus, V, pv, pq, pvpq)
 269
 270         # Update step.
 271         dx = -1 * spsolve(J, F)
 272 #        dx = -1 * linalg.lstsq(J.todense(), F)[0]
 273
 274         # Update voltage vector.
 275         npv = len(pv)
 276         npq = len(pq)
 277         if npv > 0:
 278             Va[pv] = Va[pv] + dx[range(npv)]
 279         if npq > 0:
 280             Va[pq] = Va[pq] + dx[range(npv, npv + npq)]
 281             Vm[pq] = Vm[pq] + dx[range(npv + npq, npv + npq + npq)]
 282
 283         V = Vm * exp(1j * Va)
 284         Vm = abs(V) # Avoid wrapped round negative Vm.
 285         Va = angle(V)
 286
 287         return V, Vm, Va
 288
 289     #--------------------------------------------------------------------------
 290     #  Evaluate F(x):
 291     #--------------------------------------------------------------------------
 292
 293     def _evaluate_function(self, Ybus, V, Sbus, pv, pq):
 294         """ Evaluates F(x).
 295         """
 296         mis = multiply(V, conj(Ybus * V)) - Sbus
 297
 298         F = r_[mis[pv].real, mis[pq].real, mis[pq].imag]
 299
 300         return F
 301
 302     #--------------------------------------------------------------------------
 303     #  Check convergence:
 304     #--------------------------------------------------------------------------
 305
 306     def _check_convergence(self, F):
 307         """ Checks if the solution has converged to within the specified
 308             tolerance.
 309         """
 310         normF = linalg.norm(F, Inf)
 311
 312         if normF < self.tolerance:
 313             converged = True
 314         else:
 315             converged = False
 316             if self.verbose:
 317                 logger.info("Difference: %.3f" % (normF - self.tolerance))
 318
 319         return converged
 320
 321     #--------------------------------------------------------------------------
 322     #  Evaluate Jacobian:
 323     #--------------------------------------------------------------------------
 324
 325     def _build_jacobian(self, Ybus, V, pv, pq, pvpq):
 326         """ Returns the Jacobian matrix.
 327         """
 328         pq_col = [[i] for i in pq]
 329         pvpq_col = [[i] for i in pvpq]
 330
 331         dS_dVm, dS_dVa = self.case.dSbus_dV(Ybus, V)
 332
 333         J11 = dS_dVa[pvpq_col, pvpq].real
 334
 335         J12 = dS_dVm[pvpq_col, pq].real
 336         J21 = dS_dVa[pq_col, pvpq].imag
 337         J22 = dS_dVm[pq_col, pq].imag
 338
 339         J = vstack([
 340             hstack([J11, J12]),
 341             hstack([J21, J22])
 342         ], format="csr")
 343
 344         return J
 345
 346 #------------------------------------------------------------------------------
 347 #  "FastDecoupledPF" class:
 348 #------------------------------------------------------------------------------
 349
 350 class FastDecoupledPF(_ACPF):
 351     """ Solves the power flow using fast decoupled method.
 352
 353     Based on fdpf.m from MATPOWER by Ray Zimmerman, developed at PSERC
 354     Cornell. See U{http://www.pserc.cornell.edu/matpower/} for more info.
 355     """
 356     #--------------------------------------------------------------------------
 357     #  "object" interface:
 358     #--------------------------------------------------------------------------
 359
 360     def __init__(self, case, qlimit=False, tolerance=1e-08, iter_max=20,
 361                  verbose=True, method=XB):
 362         """ Initialises a new ACPF instance.
 363         """
 364         super(FastDecoupledPF, self).__init__(case, qlimit, tolerance,
 365                                               iter_max, verbose)
 366         #: Use XB or BX method?
 367         self.method = method
 368
 369
 370     def _run_power_flow(self, Ybus, Sbus, V, pv, pq, pvpq):
 371         """ Solves the power flow using a full Newton's method.
 372         """
 373         i = 0
 374         Va = angle(V)
 375         Vm = abs(V)
 376
 377         # FIXME: Do not repeat build for each Q limit loop.
 378         Bp, Bpp = self.case.makeB(method=self.method)
 379
 380         # Evaluate initial mismatch.
 381         P, Q = self._evaluate_mismatch(Ybus, V, Sbus, pq, pvpq)
 382
 383         if self.verbose:
 384             logger.info("iteration     max mismatch (p.u.)  \n")
 385             logger.info("type   #        P            Q     \n")
 386             logger.info("---- ----  -----------  -----------\n")
 387
 388         # Check tolerance.
 389         converged = self._check_convergence(P, Q, i, "P")
 390
 391         if converged and self.verbose:
 392             logger.info("Converged!")
 393
 394         # Reduce B matrices.
 395         pq_col = [[k] for k in pq]
 396         pvpq_col = [[k] for k in pvpq]
 397         Bp = Bp[pvpq_col, pvpq].tocsc() # splu requires a CSC matrix
 398         Bpp = Bpp[pq_col, pq].tocsc()
 399
 400         # Factor B matrices.
 401         Bp_solver = splu(Bp)
 402         Bpp_solver = splu(Bpp)
 403 #        L = decomp.lu(Bp.todense())
 404 #        LU, P = decomp.lu_factor(Bp.todense())
 405
 406         # Perform Newton iterations.
 407         while (not converged) and (i < self.iter_max):
 408             i += 1
 409             # Perform P iteration, update Va.
 410             V, Vm, Va = self._p_iteration(P, Bp_solver, Vm, Va, pvpq)
 411
 412             # Evalute mismatch.
 413             P, Q = self._evaluate_mismatch(Ybus, V, Sbus, pq, pvpq)
 414             # Check tolerance.
 415             converged = self._check_convergence(P, Q, i, "P")
 416
 417             if self.verbose and converged:
 418                 logger.info("Fast-decoupled power flow converged in %d "
 419                     "P-iterations and %d Q-iterations." % (i, i - 1))
 420                 break
 421
 422             # Perform Q iteration, update Vm.
 423             V, Vm, Va = self._q_iteration(Q, Bpp_solver, Vm, Va, pq)
 424
 425             # Evalute mismatch.
 426             P, Q = self._evaluate_mismatch(Ybus, V, Sbus, pq, pvpq)
 427             # Check tolerance.
 428             converged = self._check_convergence(P, Q, i, "Q")
 429
 430             if self.verbose and converged:
 431                 logger.info("Fast-decoupled power flow converged in %d "
 432                     "P-iterations and %d Q-iterations." % (i, i))
 433                 break
 434
 435         if self.verbose and not converged:
 436             logger.error("FDPF did not converge in %d iterations." % i)
 437
 438         return V, converged, i
 439
 440     #--------------------------------------------------------------------------
 441     #  Evaluate mismatch:
 442     #--------------------------------------------------------------------------
 443
 444     def _evaluate_mismatch(self, Ybus, V, Sbus, pq, pvpq):
 445         """ Evaluates the mismatch.
 446         """
 447         mis = (multiply(V, conj(Ybus * V)) - Sbus) / abs(V)
 448
 449         P = mis[pvpq].real
 450         Q = mis[pq].imag
 451
 452         return P, Q
 453
 454     #--------------------------------------------------------------------------
 455     #  Check convergence:
 456     #--------------------------------------------------------------------------
 457
 458     def _check_convergence(self, P, Q, i, type):
 459         """ Checks if the solution has converged to within the specified
 460         tolerance.
 461         """
 462         normP = linalg.norm(P, Inf)
 463         normQ = linalg.norm(Q, Inf)
 464
 465         if self.verbose:
 466             logger.info("  %s  %3d   %10.3e   %10.3e" % (type,i, normP, normQ))
 467
 468         if (normP < self.tolerance) and (normQ < self.tolerance):
 469             converged = True
 470         else:
 471             converged = False
 472
 473         return converged
 474
 475     #--------------------------------------------------------------------------
 476     #  P iterations:
 477     #--------------------------------------------------------------------------
 478
 479     def _p_iteration(self, P, Bp_solver, Vm, Va, pvpq):
 480         """ Performs a P iteration, updates Va.
 481         """
 482         dVa = -Bp_solver.solve(P)
 483
 484         # Update voltage.
 485         Va[pvpq] = Va[pvpq] + dVa
 486         V = Vm * exp(1j * Va)
 487
 488         return V, Vm, Va
 489
 490     #--------------------------------------------------------------------------
 491     #  Q iterations:
 492     #--------------------------------------------------------------------------
 493
 494     def _q_iteration(self, Q, Bpp_solver, Vm, Va, pq):
 495         """ Performs a Q iteration, updates Vm.
 496         """
 497         dVm = -Bpp_solver.solve(Q)
 498
 499         # Update voltage.
 500         Vm[pq] = Vm[pq] + dVm
 501         V = Vm * exp(1j * Va)
 502
 503         return V, Vm, Va
 504
 505 # EOF -------------------------------------------------------------------------