Based on makeYbus.m from MATPOWER by Ray Zimmerman, developed at
        PSERC Cornell. See U{http://www.pserc.cornell.edu/matpower/} for more
        information.

        @rtype: tuple
        @return: A triple consisting of the bus admittance matrix (i.e. for all
        buses) and the matrices Yf and Yt which, when multiplied by a complex
        voltage vector, yield the vector currents injected into each line from
        the "from" and "to" buses respectively of each line.

def getYbus(self, buses=None, branches=None):
        """ Based on makeYbus.m from MATPOWER by Ray Zimmerman, developed at
        PSERC Cornell. See U{http://www.pserc.cornell.edu/matpower/} for more
        information.

        @rtype: tuple
        @return: A triple consisting of the bus admittance matrix (i.e. for all
        buses) and the matrices Yf and Yt which, when multiplied by a complex
        voltage vector, yield the vector currents injected into each line from
        the "from" and "to" buses respectively of each line.
        """
        buses = self.buses if buses is None else buses
        branches = self.branches if branches is None else branches

        nb = len(buses)
        nl = len(branches)
        ib = array(range(nb), dtype=int32)
        il = array(range(nl), dtype=int32)

        online = array([e.online for e in branches])

        # Series admittance.
        r = array([e.r for e in branches])
        x = array([e.x for e in branches])
        Ys = online / (r + 1j * x)

        # Line charging susceptance.
        b = array([e.b for e in branches])
        Bc = online * b

        #  Transformer tap ratios.
        tap = ones(nl) # Default tap ratio = 1.0.
        # Indices of branches with non-zero tap ratio.
        i_trx = array([i for i, e in enumerate(branches) if e.ratio != 0.0],
                      dtype=int32)
        # Transformer off nominal turns ratio ( = 0 for lines ) (taps at
        # "from" bus, impedance at 'to' bus, i.e. ratio = Vf / Vt)"
        ratio = array([e.ratio for e in branches])

        # Set non-zero tap ratios.
        if len(i_trx) > 0:
            tap[i_trx] = ratio[i_trx]

        # Phase shifters.
        shift = array([e.phase_shift * pi / 180.0 for e in branches])

        tap = tap * exp(1j * shift)

        # Branch admittance matrix elements.
        Ytt = Ys + 1j * Bc / 2.0
        Yff = Ytt / (tap * conj(tap))
        Yft = -Ys / conj(tap)
        Ytf = -Ys / tap

        # Shunt admittance.
        g_shunt = array([v.g_shunt for v in buses])
        b_shunt = array([v.b_shunt for v in buses])
        Ysh = (g_shunt + 1j * b_shunt) / self.base_mva

        # Connection matrices.
        f = [e.from_bus._i for e in branches]
        t = [e.to_bus._i for e in branches]

        Cf = csc_matrix((ones(nl), (il, f)), shape=(nl, nb))
        Ct = csc_matrix((ones(nl), (il, t)), shape=(nl, nb))

        # Build bus admittance matrix
        i = r_[il, il]
        j = r_[f, t]
        Yf = csc_matrix((r_[Yff, Yft], (i, j)), (nl, nb))
        Yt = csc_matrix((r_[Ytf, Ytt], (i, j)), (nl, nb))

        # Branch admittances plus shunt admittances.
        Ysh_diag = csc_matrix((Ysh, (ib, ib)), shape=(nb, nb))
        Ybus = Cf.T * Yf + Ct.T * Yt + Ysh_diag
        assert Ybus.shape == (nb, nb)

        return Ybus, Yf, Yt