pylon/dc_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 a solver for DC power flow.
  19 """
  20
  21 #------------------------------------------------------------------------------
  22 #  Imports:
  23 #------------------------------------------------------------------------------
  24
  25 import time
  26 import logging
  27 import math
  28
  29 from numpy import array, linalg, pi, r_, ix_
  30
  31 from scipy.sparse.linalg import spsolve
  32
  33 from pylon.case import REFERENCE, PV, PQ
  34
  35 #------------------------------------------------------------------------------
  36 #  Logging:
  37 #------------------------------------------------------------------------------
  38
  39 logger = logging.getLogger(__name__)
  40
  41 #------------------------------------------------------------------------------
  42 #  "DCPF" class:
  43 #------------------------------------------------------------------------------
  44
  45 class DCPF(object):
  46     """ Solves DC power flow.
  47
  48     Based on dcpf.m from MATPOWER by Ray Zimmerman, developed at PSERC
  49     Cornell. See U{http://www.pserc.cornell.edu/matpower/} for more info.
  50     """
  51
  52     #--------------------------------------------------------------------------
  53     #  "object" interface:
  54     #--------------------------------------------------------------------------
  55
  56     def __init__(self, case):
  57         """ Initialises a DCPF instance.
  58         """
  59         #: Solved case.
  60         self.case = case
  61
  62         #: Vector of voltage phase angles.
  63         self.v_angle = None
  64
  65
  66     def solve(self):
  67         """ Solves a DC power flow.
  68         """
  69         case = self.case
  70         logger.info("Starting DC power flow [%s]." % case.name)
  71         t0 = time.time()
  72         # Update bus indexes.
  73         self.case.index_buses()
  74
  75         # Find the index of the refence bus.
  76         ref_idx = self._get_reference_index(case)
  77         if ref_idx < 0:
  78             return False
  79
  80         # Build the susceptance matrices.
  81         B, Bsrc, p_businj, p_srcinj = case.Bdc
  82         # Get the vector of initial voltage angles.
  83         v_angle_guess = self._get_v_angle_guess(case)
  84         # Calculate the new voltage phase angles.
  85         v_angle, p_ref = self._get_v_angle(case, B, v_angle_guess, p_businj,
  86                                            ref_idx)
  87         logger.debug("Bus voltage phase angles: \n%s" % v_angle)
  88         self.v_angle = v_angle
  89
  90         # Push the results to the case.
  91         self._update_model(case, B, Bsrc, v_angle, p_srcinj, p_ref, ref_idx)
  92
  93         logger.info("DC power flow completed in %.3fs." % (time.time() - t0))
  94
  95         return True
  96
  97     #--------------------------------------------------------------------------
  98     #  Reference bus index:
  99     #--------------------------------------------------------------------------
 100
 101     def _get_reference_index(self, case):
 102         """ Returns the index of the reference bus.
 103         """
 104         refs = [bus._i for bus in case.connected_buses if bus.type == REFERENCE]
 105         if len(refs) == 1:
 106             return refs [0]
 107         else:
 108             logger.error("Single swing bus required for DCPF.")
 109             return -1
 110
 111     #--------------------------------------------------------------------------
 112     #  Build voltage phase angle guess vector:
 113     #--------------------------------------------------------------------------
 114
 115     def _get_v_angle_guess(self, case):
 116         """ Make the vector of voltage phase guesses.
 117         """
 118         v_angle = array([bus.v_angle * (pi / 180.0)
 119                          for bus in case.connected_buses])
 120         return v_angle
 121
 122     #--------------------------------------------------------------------------
 123     #  Calculate voltage angles:
 124     #--------------------------------------------------------------------------
 125
 126     def _get_v_angle(self, case, B, v_angle_guess, p_businj, iref):
 127         """ Calculates the voltage phase angles.
 128         """
 129         buses = case.connected_buses
 130
 131         pv_idxs = [bus._i for bus in buses if bus.type == PV]
 132         pq_idxs = [bus._i for bus in buses if bus.type == PQ]
 133         pvpq_idxs = pv_idxs + pq_idxs
 134         pvpq_rows = [[i] for i in pvpq_idxs]
 135
 136         # Get the susceptance matrix with the column and row corresponding to
 137         # the reference bus removed.
 138         Bpvpq = B[pvpq_rows, pvpq_idxs]
 139
 140         Bref = B[pvpq_rows, [iref]]
 141
 142         # Bus active power injections (generation - load) adjusted for phase
 143         # shifters and real shunts.
 144         p_surplus = array([case.s_surplus(v).real for v in buses])
 145         g_shunt = array([bus.g_shunt for bus in buses])
 146         Pbus = (p_surplus - p_businj - g_shunt) / case.base_mva
 147
 148         Pbus.shape = len(Pbus), 1
 149
 150         A = Bpvpq
 151         b = Pbus[pvpq_idxs] - Bref * v_angle_guess[iref]
 152
 153 #        x, res, rank, s = linalg.lstsq(A.todense(), b)
 154         x = spsolve(A, b)
 155
 156         # Insert the reference voltage angle of the slack bus.
 157         v_angle = r_[x[:iref], v_angle_guess[iref], x[iref:]]
 158
 159         return v_angle, Pbus[iref]
 160
 161     #--------------------------------------------------------------------------
 162     #  Update model with solution:
 163     #--------------------------------------------------------------------------
 164
 165     def _update_model(self, case, B, Bsrc, v_angle, p_srcinj, p_ref, ref_idx):
 166         """ Updates the case with values computed from the voltage phase
 167             angle solution.
 168         """
 169         iref = ref_idx
 170         base_mva = case.base_mva
 171         buses = case.connected_buses
 172         branches = case.online_branches
 173
 174         p_from = (Bsrc * v_angle + p_srcinj) * base_mva
 175         p_to = -p_from
 176
 177         for i, branch in enumerate(branches):
 178             branch.p_from = p_from[i]
 179             branch.p_to = p_to[i]
 180             branch.q_from = 0.0
 181             branch.q_to = 0.0
 182
 183         for j, bus in enumerate(buses):
 184             bus.v_angle = v_angle[j] * (180 / pi)
 185             bus.v_magnitude = 1.0
 186
 187         # Update Pg for swing generator.
 188         g_ref = [g for g in case.generators if g.bus == buses[iref]][0]
 189         # Pg = Pinj + Pload + Gs
 190         # newPg = oldPg + newPinj - oldPinj
 191         p_inj = (B[iref, :] * v_angle - p_ref) * base_mva
 192         g_ref.p += p_inj[0]
 193
 194 # EOF -------------------------------------------------------------------------