Commit ⋅ e2nIEE/pandapower

    import logging

logger = logging.getLogger(__name__)
# import time

import numpy as np
from scipy.sparse.linalg import factorized

from pandapower.auxiliary import _clean_up, _add_ppc_options, _add_sc_options, _add_auxiliary_elements
from pandapower.pd2ppc import _pd2ppc
from pandapower.pd2ppc_zero import _pd2ppc_zero
from pandapower.results import _copy_results_ppci_to_ppc
from pandapower.shortcircuit.currents import _calc_ikss, _calc_ikss_1ph, _calc_ip, _calc_ith, _calc_branch_currents, \
    _calc_single_bus_sc
from pandapower.shortcircuit.currents import _calc_ikss, \
    _calc_ikss_1ph, _calc_ip, _calc_ith, _calc_branch_currents, \
    _calc_single_bus_sc, _calc_single_bus_sc_no_y_inv
from pandapower.shortcircuit.impedance import _calc_zbus, _calc_ybus, _calc_rx
from pandapower.shortcircuit.kappa import _add_kappa_to_ppc
from pandapower.shortcircuit.results import _extract_results, _extract_single_results

@@ -26,7 +29,8 @@


def calc_sc(net, fault="3ph", case='max', lv_tol_percent=10, topology="auto", ip=False,
            ith=False, tk_s=1., kappa_method="C", r_fault_ohm=0., x_fault_ohm=0.,
            branch_results=False, check_connectivity=True, return_all_currents=False):
            branch_results=False, check_connectivity=True, return_all_currents=False,
            bus=None, inverse_y=True):
    """
    Calculates minimal or maximal symmetrical short-circuit currents.
    The calculation is based on the method of the equivalent voltage source

@@ -85,6 +89,10 @@

        **return_all_currents** (bool, False) applies only if branch_results=True, if True short-circuit currents for
        each (branch, bus) tuple is returned otherwise only the max/min is returned

        **bus** (int, list, np.array, None) defines if short-circuit calculations should only be calculated for defined bus

        **inverse_y** (bool, True) defines if complete inverse should be used instead of LU factorization, factorization version is in experiment which should be faster and memory efficienter


    OUTPUT:


@@ -112,6 +120,12 @@

        logger.warning("Branch results are in beta mode and might not always be reliable, "
                       "especially for transformers")

    # Convert bus to numpy array for better performance
    if isinstance(bus, int):
        bus = np.array([bus])
    elif isinstance(bus, list):
        bus = np.array(bus)

    kappa = ith or ip
    net["_options"] = {}
    _add_ppc_options(net, calculate_voltage_angles=False, trafo_model="pi",

@@ -121,17 +135,19 @@

    _add_sc_options(net, fault=fault, case=case, lv_tol_percent=lv_tol_percent, tk_s=tk_s,
                    topology=topology, r_fault_ohm=r_fault_ohm, kappa_method=kappa_method,
                    x_fault_ohm=x_fault_ohm, kappa=kappa, ip=ip, ith=ith,
                    branch_results=branch_results, return_all_currents=return_all_currents)
                    branch_results=branch_results, return_all_currents=return_all_currents,
                    inverse_y=inverse_y)
    init_results(net, "sc")
    if fault == "3ph":
        _calc_sc(net)
    if fault == "2ph":
        _calc_sc(net)
    if fault == "1ph":
        _calc_sc_1ph(net)
    if fault in ("2ph", "3ph"):
        _calc_sc(net, bus)
    elif fault == "1ph":
        _calc_sc_1ph(net, bus)
    else:
        raise ValueError("Invalid fault %s" % fault)


def calc_single_sc(net, bus, fault="3ph", case='max', lv_tol_percent=10, check_connectivity=True):
def calc_single_sc(net, bus, fault="3ph", case='max', lv_tol_percent=10,
                   check_connectivity=True, inverse_y=True):
    """
    Calculates minimal or maximal symmetrical short-circuit currents.
    The calculation is based on the method of the equivalent voltage source

@@ -190,9 +206,10 @@

    _add_sc_options(net, fault=fault, case=case, lv_tol_percent=lv_tol_percent, tk_s=1.,
                    topology="auto", r_fault_ohm=0., kappa_method="C",
                    x_fault_ohm=0., kappa=False, ip=False, ith=False,
                    branch_results=True, return_all_currents=False)
                    branch_results=True, return_all_currents=False,
                    inverse_y=inverse_y)
    init_results(net, "sc")
    if fault == "3ph" or fault == "2ph":
    if fault in ("2ph", "3ph"):
        _calc_sc_single(net, bus)
    elif fault == "1ph":
        raise NotImplementedError("1ph short-circuits are not yet implemented")

@@ -204,48 +221,64 @@

    _add_auxiliary_elements(net)
    ppc, ppci = _pd2ppc(net)
    _calc_ybus(ppci)
    try:
        _calc_zbus(ppci)
    except Exception as e:
        _clean_up(net, res=False)
        raise (e)
    _calc_rx(net, ppci)
    _calc_ikss(net, ppci)
    _calc_single_bus_sc(net, ppci, bus)

    if net["_options"]["inverse_y"]:
        _calc_zbus(net, ppci)
        _calc_rx(net, ppci, bus=None)
        _calc_ikss(net, ppci, bus=None)
        _calc_single_bus_sc(net, ppci, bus)
    else:
        # Factorization Ybus once
        ppci["internal"]["ybus_fact"] = factorized(ppci["internal"]["Ybus"])

        _calc_rx(net, ppci, bus)
        _calc_ikss(net, ppci, bus)
        _calc_single_bus_sc_no_y_inv(net, ppci, bus)

        # Delete factorization object
        ppci["internal"].pop("ybus_fact")

    ppc = _copy_results_ppci_to_ppc(ppci, ppc, "sc")
    _extract_single_results(net, ppc)
    _clean_up(net)


def _calc_sc(net):
def _calc_sc(net, bus):
    _add_auxiliary_elements(net)
    ppc, ppci = _pd2ppc(net)
    _calc_ybus(ppci)
    try:
        _calc_zbus(ppci)
    except Exception as e:
        _clean_up(net, res=False)
        raise (e)
    _calc_rx(net, ppci)
    _add_kappa_to_ppc(net, ppci)
    _calc_ikss(net, ppci)

    if net["_options"]["inverse_y"]:
        _calc_zbus(net, ppci)
    else:
        # Factorization Ybus once
        ppci["internal"]["ybus_fact"] = factorized(ppci["internal"]["Ybus"])

    _calc_rx(net, ppci, bus)

    # kappa required inverse of Zbus, which is optimized
    if net["_options"]["kappa"]:
        _add_kappa_to_ppc(net, ppci)
    _calc_ikss(net, ppci, bus)

    if net["_options"]["ip"]:
        _calc_ip(net, ppci)
    if net["_options"]["ith"]:
        _calc_ith(net, ppci)

    if net._options["branch_results"]:
        _calc_branch_currents(net, ppci)
        _calc_branch_currents(net, ppci, bus)

    ppc = _copy_results_ppci_to_ppc(ppci, ppc, "sc")

    if net["_options"]["return_all_currents"]:
        _extract_results(net, ppc, ppc_0=None)
    else:
        _extract_results(net, ppc, ppc_0=None)
    _extract_results(net, ppc, ppc_0=None, bus=bus)
    _clean_up(net)

    if "ybus_fact" in ppci["internal"]:
        # Delete factorization object
        ppci["internal"].pop("ybus_fact")

def _calc_sc_1ph(net):

def _calc_sc_1ph(net, bus):
    """
    calculation method for single phase to ground short-circuit currents
    """

@@ -253,26 +286,28 @@

    # pos. seq bus impedance
    ppc, ppci = _pd2ppc(net)
    _calc_ybus(ppci)
    try:
        _calc_zbus(ppci)
    except Exception as e:
        _clean_up(net, res=False)
        raise (e)
    _calc_rx(net, ppci)
    _add_kappa_to_ppc(net, ppci)

    # zero seq bus impedance
    ppc_0, ppci_0 = _pd2ppc_zero(net)
    _calc_ybus(ppci_0)
    try:
        _calc_zbus(ppci_0)
    except Exception as e:
        _clean_up(net, res=False)
        raise (e)
    _calc_rx(net, ppci_0)
    _calc_ikss_1ph(net, ppci, ppci_0)

    if net["_options"]["inverse_y"]:
        _calc_zbus(net, ppci)
        _calc_zbus(net, ppci_0)
    else:
        # Factorization Ybus once
        ppci["internal"]["ybus_fact"] = factorized(ppci["internal"]["Ybus"])
        ppci_0["internal"]["ybus_fact"] = factorized(ppci_0["internal"]["Ybus"])

    _calc_rx(net, ppci, bus=bus)
    _add_kappa_to_ppc(net, ppci)

    _calc_rx(net, ppci_0, bus=bus)
    _calc_ikss_1ph(net, ppci, ppci_0, bus=bus)

    if net._options["branch_results"]:
        _calc_branch_currents(net, ppci)
        _calc_branch_currents(net, ppci, bus=bus)
    ppc_0 = _copy_results_ppci_to_ppc(ppci_0, ppc_0, "sc")
    ppc = _copy_results_ppci_to_ppc(ppci, ppc, "sc")
    _extract_results(net, ppc, ppc_0)
    _extract_results(net, ppc, ppc_0, bus=bus)
    _clean_up(net)

@@ -12,33 +12,51 @@

from pandapower.pypower.idx_gen import GEN_BUS, MBASE
from pandapower.shortcircuit.idx_brch import IKSS_F, IKSS_T, IP_F, IP_T, ITH_F, ITH_T
from pandapower.shortcircuit.idx_bus import C_MIN, C_MAX, KAPPA, R_EQUIV, IKSS1, IP, ITH, X_EQUIV, IKSS2, IKCV, M
from pandapower.shortcircuit.impedance import _calc_zbus_diag

from pandapower.pypower.pfsoln import pfsoln as pfsoln_pypower
from pandapower.pf.ppci_variables import _get_pf_variables_from_ppci

def _calc_ikss(net, ppc, bus=None):
    # Vectorized for multiple bus
    if bus is None:
        # Slice(None) is equal to : select
        bus_idx = slice(None)
    else:
        bus_idx = net._pd2ppc_lookups["bus"][bus] #bus where the short-circuit is calculated (j)

def _calc_ikss(net, ppc):
    fault = net._options["fault"]
    case = net._options["case"]
    c = ppc["bus"][:, C_MIN] if case == "min" else ppc["bus"][:, C_MAX]
    c = ppc["bus"][bus_idx, C_MIN] if case == "min" else ppc["bus"][bus_idx, C_MAX]
    ppc["internal"]["baseI"] = ppc["bus"][:, BASE_KV] * np.sqrt(3) / ppc["baseMVA"]
    z_equiv = abs(ppc["bus"][:, R_EQUIV] + ppc["bus"][:, X_EQUIV] * 1j)

    z_equiv = abs(ppc["bus"][bus_idx, R_EQUIV] + ppc["bus"][bus_idx, X_EQUIV] * 1j)
    if fault == "3ph":
        ppc["bus"][:, IKSS1] = c / z_equiv / ppc["bus"][:, BASE_KV] / np.sqrt(3) * ppc["baseMVA"]
        ppc["bus"][bus_idx, IKSS1] = c / z_equiv / ppc["bus"][bus_idx, BASE_KV] / np.sqrt(3) * ppc["baseMVA"]
    elif fault == "2ph":
        ppc["bus"][:, IKSS1] = c / z_equiv / ppc["bus"][:, BASE_KV] / 2 * ppc["baseMVA"]
        ppc["bus"][bus_idx, IKSS1] = c / z_equiv / ppc["bus"][bus_idx, BASE_KV] / 2 * ppc["baseMVA"]
    _current_source_current(net, ppc)


def _calc_ikss_1ph(net, ppc, ppc_0):
    fault = net._options["fault"]
def _calc_ikss_1ph(net, ppc, ppc_0, bus=None):
    # Vectorized for multiple bus
    if bus is None:
        # Slice(None) is equal to : select
        bus_idx = slice(None)
    else:
        bus_idx = net._pd2ppc_lookups["bus"][bus] #bus where the short-circuit is calculated (j)

    case = net._options["case"]
    c = ppc["bus"][:, C_MIN] if case == "min" else ppc["bus"][:, C_MAX]
    c = ppc["bus"][bus_idx, C_MIN] if case == "min" else ppc["bus"][bus_idx, C_MAX]
    ppc["internal"]["baseI"] = ppc["bus"][:, BASE_KV] * np.sqrt(3) / ppc["baseMVA"]
    ppc_0["internal"]["baseI"] = ppc_0["bus"][:, BASE_KV] * np.sqrt(3) / ppc_0["baseMVA"]
    z_equiv = abs((ppc["bus"][:, R_EQUIV] + ppc["bus"][:, X_EQUIV] * 1j) * 2 +
                  (ppc_0["bus"][:, R_EQUIV] + ppc_0["bus"][:, X_EQUIV] * 1j))
    if fault == "1ph":
        ppc_0["bus"][:, IKSS1] = c / z_equiv / ppc_0["bus"][:, BASE_KV] * np.sqrt(3) * ppc_0["baseMVA"]
        ppc["bus"][:, IKSS1] = c / z_equiv / ppc_0["bus"][:, BASE_KV] * np.sqrt(3) * ppc_0["baseMVA"]
        _current_source_current(net, ppc)

    z_equiv = abs((ppc["bus"][bus_idx, R_EQUIV] + ppc["bus"][bus_idx, X_EQUIV] * 1j) * 2 +
                  (ppc_0["bus"][bus_idx, R_EQUIV] + ppc_0["bus"][bus_idx, X_EQUIV] * 1j))

    ppc_0["bus"][bus_idx, IKSS1] = c / z_equiv / ppc_0["bus"][bus_idx, BASE_KV] * np.sqrt(3) * ppc_0["baseMVA"]
    ppc["bus"][bus_idx, IKSS1] = c / z_equiv / ppc["bus"][bus_idx, BASE_KV] * np.sqrt(3) * ppc["baseMVA"]
    _current_source_current(net, ppc)


def _current_source_current(net, ppc):

@@ -56,13 +74,18 @@

    baseI = ppc["internal"]["baseI"]
    sgen_buses = sgen.bus.values
    sgen_buses_ppc = bus_lookup[sgen_buses]
    Zbus = ppc["internal"]["Zbus"]
    if not "k" in sgen:
        raise ValueError("Nominal to short-circuit current has to specified in net.sgen.k")
    i_sgen_pu = sgen.sn_mva.values / net.sn_mva * sgen.k.values
    buses, ikcv_pu, _ = _sum_by_group(sgen_buses_ppc, i_sgen_pu, i_sgen_pu)
    ppc["bus"][buses, IKCV] = ikcv_pu
    ppc["bus"][:, IKSS2] = abs(1 / np.diag(Zbus) * np.dot(Zbus, ppc["bus"][:, IKCV] * -1j) / baseI)
    if net["_options"]["inverse_y"]:
        Zbus = ppc["internal"]["Zbus"]
        ppc["bus"][:, IKSS2] = abs(1 / np.diag(Zbus) * np.dot(Zbus, ppc["bus"][:, IKCV] * -1j) / baseI)
    else:
        ybus_fact = ppc["internal"]["ybus_fact"]
        diagZ = _calc_zbus_diag(net, ppc)
        ppc["bus"][:, IKSS2] = abs(ybus_fact(ppc["bus"][:, IKCV] * -1j) / diagZ / baseI)
    ppc["bus"][buses, IKCV] /= baseI[buses]



@@ -83,84 +106,6 @@

    ppc["bus"][:, ITH] = ith


def _calc_branch_currents(net, ppc):
    case = net._options["case"]
    Zbus = ppc["internal"]["Zbus"]
    Yf = ppc["internal"]["Yf"]
    Yt = ppc["internal"]["Yt"]
    baseI = ppc["internal"]["baseI"]
    n = ppc["bus"].shape[0]
    fb = np.real(ppc["branch"][:, 0]).astype(int)
    tb = np.real(ppc["branch"][:, 1]).astype(int)
    minmax = np.nanmin if case == "min" else np.nanmax

    # calculate voltage source branch current
    V_ikss = (ppc["bus"][:, IKSS1] * baseI) * Zbus
    ikss1_all_f = np.conj(Yf.dot(V_ikss))
    ikss1_all_t = np.conj(Yt.dot(V_ikss))
    ikss1_all_f[abs(ikss1_all_f) < 1e-10] = 0.
    ikss1_all_t[abs(ikss1_all_t) < 1e-10] = 0.

    # add current source branch current if there is one
    current_sources = any(ppc["bus"][:, IKCV]) > 0
    if current_sources:
        current = np.tile(-ppc["bus"][:, IKCV], (n, 1))
        np.fill_diagonal(current, current.diagonal() + ppc["bus"][:, IKSS2])
        V = np.dot((current * baseI), Zbus).T
        fb = np.real(ppc["branch"][:, 0]).astype(int)
        tb = np.real(ppc["branch"][:, 1]).astype(int)
        ikss2_all_f = np.conj(Yf.dot(V))
        ikss2_all_t = np.conj(Yt.dot(V))
        ikss_all_f = abs(ikss1_all_f + ikss2_all_f)
        ikss_all_t = abs(ikss1_all_t + ikss2_all_t)
    else:
        ikss_all_f = abs(ikss1_all_f)
        ikss_all_t = abs(ikss1_all_t)

    if net._options["return_all_currents"]:
        ppc["internal"]["branch_ikss_f"] = ikss_all_f / baseI[fb, None]
        ppc["internal"]["branch_ikss_t"] = ikss_all_t / baseI[tb, None]
    else:
        ikss_all_f[ikss_all_f < 1e-10] = np.nan
        ikss_all_t[ikss_all_t < 1e-10] = np.nan
        ppc["branch"][:, IKSS_F] = minmax(ikss_all_f, axis=1) / baseI[fb]
        ppc["branch"][:, IKSS_T] = minmax(ikss_all_t, axis=1) / baseI[tb]

    if net._options["ip"]:
        kappa = ppc["bus"][:, KAPPA]
        if current_sources:
            ip_all_f = np.sqrt(2) * (ikss1_all_f * kappa + ikss2_all_f)
            ip_all_t = np.sqrt(2) * (ikss1_all_t * kappa + ikss2_all_t)
        else:
            ip_all_f = np.sqrt(2) * ikss1_all_f * kappa
            ip_all_t = np.sqrt(2) * ikss1_all_t * kappa

        if net._options["return_all_currents"]:
            ppc["internal"]["branch_ip_f"] = abs(ip_all_f) / baseI[fb, None]
            ppc["internal"]["branch_ip_t"] = abs(ip_all_t) / baseI[tb, None]
        else:
            ip_all_f[abs(ip_all_f) < 1e-10] = np.nan
            ip_all_t[abs(ip_all_t) < 1e-10] = np.nan
            ppc["branch"][:, IP_F] = minmax(abs(ip_all_f), axis=1) / baseI[fb]
            ppc["branch"][:, IP_T] = minmax(abs(ip_all_t), axis=1) / baseI[tb]

    if net._options["ith"]:
        n = 1
        m = ppc["bus"][:, M]
        ith_all_f = ikss_all_f * np.sqrt(m + n)
        ith_all_t = ikss_all_t * np.sqrt(m + n)

        if net._options["return_all_currents"]:
            ppc["internal"]["branch_ith_f"] = ith_all_f / baseI[fb, None]
            ppc["internal"]["branch_ith_t"] = ith_all_t / baseI[tb, None]
        else:
            ppc["branch"][:, ITH_F] = minmax(ith_all_f, axis=1) / baseI[fb]
            ppc["branch"][:, ITH_T] = minmax(ith_all_t, axis=1) / baseI[fb]





def _calc_ib_generator(net, ppci):
    Zbus = ppci["internal"]["Zbus"]
    baseI = ppci["internal"]["baseI"]

@@ -230,7 +175,7 @@



def _calc_single_bus_sc(net, ppc, bus):
    case = net._options["case"]
    # case = net._options["case"]
    bus_idx = net._pd2ppc_lookups["bus"][bus]
    n = ppc["bus"].shape[0]
    Zbus = ppc["internal"]["Zbus"]

@@ -239,7 +184,7 @@

    baseI = ppc["internal"]["baseI"]
    #    fb = np.real(ppc["branch"][:, 0]).astype(int)
    #    tb = np.real(ppc["branch"][:, 1]).astype(int)
    c = ppc["bus"][:, C_MIN] if case == "min" else ppc["bus"][:, C_MAX]
    # c = ppc["bus"][:, C_MIN] if case == "min" else ppc["bus"][:, C_MAX]

    # calculate voltage source branch current
    V_ikss = (ppc["bus"][:, IKSS1] * baseI) * Zbus

@@ -257,15 +202,167 @@

    #    ppc["branch"][:, IKSS_T] = abs(ikss_all_t[:, bus_idx] / baseI[tb])
    calc_branch_results(net, ppc, V)

def _calc_single_bus_sc_no_y_inv(net, ppc, bus):
    # Vectorized for multiple bus
    if bus is None:
        # Slice(None) is equal to : select
        bus_idx = slice(None)
    else:
        bus_idx = net._pd2ppc_lookups["bus"][bus] #bus where the short-circuit is calculated (j)

from pandapower.pypower.pfsoln import pfsoln as pfsoln_pypower
from pandapower.pf.ppci_variables import _get_pf_variables_from_ppci
    ybus = ppc["internal"]["Ybus"]
    ybus_fact = ppc["internal"]["ybus_fact"]
    # case = net._options["case"]
    baseI = ppc["internal"]["baseI"]
    # vqj = ppc["bus"][:, C_MIN] if case == "min" else ppc["bus"][:, C_MAX] #this is the source voltage in per unit (VQj)

    # Solve Ikss from voltage source
    n_bus = ybus.shape[0]

    # ybus_sub_mask = (np.arange(ybus.shape[0]) != bus_idx)
    # V_ikss = np.zeros(n_bus, dtype=np.complex)
    # V_ikss[bus_idx] = vqj[bus_idx]

    # # Solve Ax = b
    # b = np.zeros(n_bus-1, dtype=np.complex) -\
    #     (ybus[:, ~ybus_sub_mask].toarray())[ybus_sub_mask].ravel() * V_ikss[bus_idx]
    # ybus_sub = ybus[ybus_sub_mask, :][:, ybus_sub_mask]
    # x = spsolve(ybus_sub, b)

    # V_ikss[ybus_sub_mask] = x
    # I_ikss = np.zeros(n_bus, dtype=np.complex)
    # I_ikss[bus_idx] = np.dot(ybus[bus_idx, :].toarray(), V_ikss)
    # V = V_ikss

    # Version 2
    I_ikss = np.zeros(n_bus, dtype=np.complex)
    I_ikss[bus_idx] = ppc["bus"][bus_idx, IKSS1]
    V_ikss = ybus_fact(I_ikss * baseI)
    V = V_ikss

    #TODO include current sources
    current_sources = any(ppc["bus"][:, IKCV]) > 0
    if current_sources:
        current = -ppc["bus"][:, IKCV]
        current[bus_idx] += ppc["bus"][bus_idx, IKSS2]
        V_source = ybus_fact(current)
        V += V_source

    calc_branch_results(net, ppc, V)


def calc_branch_results(net, ppci, V):
    Ybus = ppci["internal"]["Ybus"]
    Yf = ppci["internal"]["Yf"]
    Yt = ppci["internal"]["Yt"]
    baseMVA, bus, gen, branch, ref, _, pq, _, _, V0, ref_gens = _get_pf_variables_from_ppci(ppci)
    baseMVA, bus, gen, branch, ref, _, _, _, _, _, ref_gens = _get_pf_variables_from_ppci(ppci)
    bus, gen, branch = pfsoln_pypower(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, ref_gens)
    ppci["bus"], ppci["gen"], ppci["branch"] = bus, gen, branch


def _calc_branch_currents(net, ppc, bus):
    # Vectorized for multiple bus
    if bus is None:
        # Slice(None) is equal to select all
        bus = net.bus.index

    bus_idx = net._pd2ppc_lookups["bus"][bus]
    # Select only in service bus for sc calculation
    bus_idx = bus_idx[bus_idx < ppc['bus'].shape[0]]
    n_sc_bus = np.shape(bus_idx)[0]

    case = net._options["case"]

    Yf = ppc["internal"]["Yf"]
    Yt = ppc["internal"]["Yt"]
    baseI = ppc["internal"]["baseI"]
    n_bus = ppc["bus"].shape[0]
    fb = np.real(ppc["branch"][:, 0]).astype(int)
    tb = np.real(ppc["branch"][:, 1]).astype(int)
    minmax = np.nanmin if case == "min" else np.nanmax

    # calculate voltage source branch current
    if net["_options"]["inverse_y"]:
        Zbus = ppc["internal"]["Zbus"]
        V_ikss = (ppc["bus"][:, IKSS1] * baseI) * Zbus
        V_ikss = V_ikss[:, bus_idx]
    else:
        ybus_fact = ppc["internal"]["ybus_fact"]
        V_ikss = np.zeros((n_bus, n_sc_bus), dtype=np.complex)
        for ix, b in enumerate(bus_idx):
            ikss = np.zeros(n_bus, dtype=np.complex)
            ikss[b] = ppc["bus"][b, IKSS1] * baseI[b]
            V_ikss[:, ix] = ybus_fact(ikss)

    ikss1_all_f = np.conj(Yf.dot(V_ikss))
    ikss1_all_t = np.conj(Yt.dot(V_ikss))
    ikss1_all_f[abs(ikss1_all_f) < 1e-10] = 0.
    ikss1_all_t[abs(ikss1_all_t) < 1e-10] = 0.

    # add current source branch current if there is one
    current_sources = any(ppc["bus"][:, IKCV]) > 0
    if current_sources:
        current = np.tile(-ppc["bus"][:, IKCV], (n_sc_bus, 1))
        for ix, b in enumerate(bus_idx):
            current[ix, b] += ppc["bus"][b, IKSS2]

        # calculate voltage source branch current
        if net["_options"]["inverse_y"]:
            Zbus = ppc["internal"]["Zbus"]
            V = np.dot((current * baseI), Zbus).T
        else:
            ybus_fact = ppc["internal"]["ybus_fact"]
            V = np.zeros((n_bus, n_sc_bus), dtype=np.complex)
            for ix, b in enumerate(bus_idx):
                V[:, ix] = ybus_fact(current[ix, :] * baseI[b])

        fb = np.real(ppc["branch"][:, 0]).astype(int)
        tb = np.real(ppc["branch"][:, 1]).astype(int)
        ikss2_all_f = np.conj(Yf.dot(V))
        ikss2_all_t = np.conj(Yt.dot(V))

        ikss_all_f = abs(ikss1_all_f + ikss2_all_f)
        ikss_all_t = abs(ikss1_all_t + ikss2_all_t)
    else:
        ikss_all_f = abs(ikss1_all_f)
        ikss_all_t = abs(ikss1_all_t)

    if net._options["return_all_currents"]:
        ppc["internal"]["branch_ikss_f"] = ikss_all_f / baseI[fb, None]
        ppc["internal"]["branch_ikss_t"] = ikss_all_t / baseI[tb, None]
    else:
        ikss_all_f[abs(ikss_all_f) < 1e-10] = np.nan
        ikss_all_t[abs(ikss_all_t) < 1e-10] = np.nan
        ppc["branch"][:, IKSS_F] = np.nan_to_num(minmax(ikss_all_f, axis=1) / baseI[fb])
        ppc["branch"][:, IKSS_T] = np.nan_to_num(minmax(ikss_all_t, axis=1) / baseI[tb])

    if net._options["ip"]:
        kappa = ppc["bus"][:, KAPPA]
        if current_sources:
            ip_all_f = np.sqrt(2) * (ikss1_all_f * kappa[bus_idx] + ikss2_all_f)
            ip_all_t = np.sqrt(2) * (ikss1_all_t * kappa[bus_idx] + ikss2_all_t)
        else:
            ip_all_f = np.sqrt(2) * ikss1_all_f * kappa[bus_idx]
            ip_all_t = np.sqrt(2) * ikss1_all_t * kappa[bus_idx]

        if net._options["return_all_currents"]:
            ppc["internal"]["branch_ip_f"] = abs(ip_all_f) / baseI[fb, None]
            ppc["internal"]["branch_ip_t"] = abs(ip_all_t) / baseI[tb, None]
        else:
            ip_all_f[abs(ip_all_f) < 1e-10] = np.nan
            ip_all_t[abs(ip_all_t) < 1e-10] = np.nan
            ppc["branch"][:, IP_F] = np.nan_to_num(minmax(abs(ip_all_f), axis=1) / baseI[fb])
            ppc["branch"][:, IP_T] = np.nan_to_num(minmax(abs(ip_all_t), axis=1) / baseI[tb])

    if net._options["ith"]:
        n = 1
        m = ppc["bus"][bus_idx, M]
        ith_all_f = ikss_all_f * np.sqrt(m + n)
        ith_all_t = ikss_all_t * np.sqrt(m + n)

        if net._options["return_all_currents"]:
            ppc["internal"]["branch_ith_f"] = ith_all_f / baseI[fb, None]
            ppc["internal"]["branch_ith_t"] = ith_all_t / baseI[tb, None]
        else:
            ppc["branch"][:, ITH_F] = np.nan_to_num(minmax(ith_all_f, axis=1) / baseI[fb])
            ppc["branch"][:, ITH_T] = np.nan_to_num(minmax(ith_all_t, axis=1) / baseI[fb])

@@ -189,14 +189,23 @@


def _ppci_internal_to_ppc(result, ppc):
    for key, value in result["internal"].items():
        # Only for sc calculation
        # if branch current matrices have been stored they need to include out of service elements
        if key in ["branch_ikss_f", "branch_ikss_t", "branch_ip_f", "branch_ip_t", "branch_ith_f", "branch_ith_t"]:
            n_buses = np.shape(ppc['bus'])[0]
        if key in ["branch_ikss_f", "branch_ikss_t",
                   "branch_ip_f", "branch_ip_t",
                   "branch_ith_f", "branch_ith_t"]:

            # n_buses = np.shape(ppc['bus'])[0]
            n_branches = np.shape(ppc['branch'])[0]
            n_rows_result = np.shape(result['bus'])[0]
            update_matrix = np.empty((n_branches, n_buses)) * np.nan
            update_matrix[result["internal"]['branch_is'], :n_rows_result] = result["internal"][key]
            ppc['internal'][key] = np.copy(update_matrix)
            # n_rows_result = np.shape(result['bus'])[0]
            # update_matrix = np.empty((n_branches, n_buses)) * np.nan
            # update_matrix[result["internal"]['branch_is'], :n_rows_result] = result["internal"][key]

            # # To select only required buses and pad one column of nan value for oos bus
            update_matrix = np.empty((n_branches, value.shape[1]+1)) * 0.0
            update_matrix[result["internal"]['branch_is'],
                          :value.shape[1]] = result["internal"][key]
            ppc['internal'][key] = update_matrix
        else:
            ppc["internal"][key] = value


@@ -9,25 +9,32 @@


from pandapower.shortcircuit.idx_brch import IKSS_F, IKSS_T, IP_F, IP_T, ITH_F, ITH_T
from pandapower.shortcircuit.idx_bus import IKSS1, IP, ITH, IKSS2
from pandapower.pypower.idx_bus import VM, VA
from pandapower.pypower.idx_bus import VM, VA, BUS_TYPE
from pandapower.results_bus import _get_bus_idx, _set_buses_out_of_service
from pandapower.results import _get_aranged_lookup, _get_branch_results
from pandapower.shortcircuit.idx_bus import C_MIN, C_MAX

def _get_bus_ppc_idx_for_br_all_results(net, ppc, bus):
    bus_lookup = net._pd2ppc_lookups["bus"]
    if bus is None:
        bus = net.bus.index

    ppc_index = np.arange(np.shape(bus)[0])
    ppc_index[ppc["bus"][bus_lookup[ppc_index], BUS_TYPE] == 4] = -1
    return bus, ppc_index

def _extract_results(net, ppc, ppc_0):
    _get_bus_results(net, ppc, ppc_0)
def _extract_results(net, ppc, ppc_0, bus):
    _get_bus_results(net, ppc, ppc_0, bus)
    if net._options["branch_results"]:
        if net._options['return_all_currents']:
            _get_line_all_results(net, ppc)
            _get_trafo_all_results(net, ppc)
            _get_trafo3w_all_results(net, ppc)
            _get_line_all_results(net, ppc, bus)
            _get_trafo_all_results(net, ppc, bus)
            _get_trafo3w_all_results(net, ppc, bus)
        else:
            _get_line_results(net, ppc)
            _get_trafo_results(net, ppc)
            _get_trafo3w_results(net, ppc)


def _extract_single_results(net, ppc):
    for element in ["line", "trafo"]:
        net["res_%s_sc"%element] = pd.DataFrame(np.nan, index=net[element].index,

@@ -51,7 +58,10 @@

    net["res_bus_sc"]["va_degree"] = ppc["bus"][bus_idx, VA]


def _get_bus_results(net, ppc, ppc_0):
def _get_bus_results(net, ppc, ppc_0, bus):
    if bus is None:
        bus = slice(None)

    bus_lookup = net._pd2ppc_lookups["bus"]
    ppc_index = bus_lookup[net.bus.index]


@@ -64,6 +74,8 @@

    if net._options["ith"]:
        net.res_bus_sc["ith_ka"] = ppc["bus"][ppc_index, ITH]

    net.res_bus_sc = net.res_bus_sc.loc[bus, :]


def _get_line_results(net, ppc):
    branch_lookup = net._pd2ppc_lookups["branch"]

@@ -78,13 +90,13 @@

            net.res_line_sc["ith_ka"] = minmax(ppc["branch"][f:t, [ITH_F, ITH_T]].real, axis=1)


def _get_line_all_results(net, ppc):
def _get_line_all_results(net, ppc, bus):
    case = net._options["case"]
    bus_lookup = net._pd2ppc_lookups["bus"]
    ppc_index = bus_lookup[net.bus.index]

    bus, ppc_index = _get_bus_ppc_idx_for_br_all_results(net, ppc, bus)
    branch_lookup = net._pd2ppc_lookups["branch"]

    multindex = pd.MultiIndex.from_product([net.res_line_sc.index, net.bus.index], names=['line','bus'])
    multindex = pd.MultiIndex.from_product([net.res_line_sc.index, bus], names=['line','bus'])
    net.res_line_sc = net.res_line_sc.reindex(multindex)

    if "line" in branch_lookup:

@@ -100,7 +112,6 @@

            net.res_line_sc["ith_ka"] = minmax(ppc["internal"]["branch_ith_f"][f:t, ppc_index].real.reshape(-1, 1),
                                               ppc["internal"]["branch_ith_t"][f:t, ppc_index].real.reshape(-1, 1))


def _get_trafo_results(net, ppc):
    branch_lookup = net._pd2ppc_lookups["branch"]
    if "trafo" in branch_lookup:

@@ -109,12 +120,11 @@

        net.res_trafo_sc["ikss_lv_ka"] = ppc["branch"][f:t, IKSS_T].real


def _get_trafo_all_results(net, ppc):
    bus_lookup = net._pd2ppc_lookups["bus"]
    ppc_index = bus_lookup[net.bus.index]
def _get_trafo_all_results(net, ppc, bus):
    bus, ppc_index = _get_bus_ppc_idx_for_br_all_results(net, ppc, bus)
    branch_lookup = net._pd2ppc_lookups["branch"]

    multindex = pd.MultiIndex.from_product([net.res_trafo_sc.index, net.bus.index], names=['trafo', 'bus'])
    multindex = pd.MultiIndex.from_product([net.res_trafo_sc.index, bus], names=['trafo', 'bus'])
    net.res_trafo_sc = net.res_trafo_sc.reindex(multindex)

    if "trafo" in branch_lookup:

@@ -135,12 +145,11 @@

        net.res_trafo3w_sc["ikss_lv_ka"] = ppc["branch"][mv:lv, IKSS_T].real


def _get_trafo3w_all_results(net, ppc):
    bus_lookup = net._pd2ppc_lookups["bus"]
    ppc_index = bus_lookup[net.bus.index]
def _get_trafo3w_all_results(net, ppc, bus):
    bus, ppc_index = _get_bus_ppc_idx_for_br_all_results(net, ppc, bus)
    branch_lookup = net._pd2ppc_lookups["branch"]

    multindex = pd.MultiIndex.from_product([net.res_trafo3w_sc.index, net.bus.index], names=['trafo3w', 'bus'])
    multindex = pd.MultiIndex.from_product([net.res_trafo3w_sc.index, bus], names=['trafo3w', 'bus'])
    net.res_trafo3w_sc = net.res_trafo3w_sc.reindex(multindex)

    if "trafo3w" in branch_lookup:

@@ -7,6 +7,7 @@

import copy
import networkx as nx
import numpy as np
from scipy.sparse.linalg import factorized

from pandapower.pypower.idx_brch import F_BUS, T_BUS, BR_R, BR_X
from pandapower.pypower.idx_bus import BUS_I, GS, BS, BASE_KV

@@ -29,6 +30,9 @@

        raise ValueError("Unknown kappa method %s - specify B or C"%kappa_method)
    ppc["bus"][:, KAPPA] = kappa

def _kappa(rx):
    return 1.02 + .98 * np.exp(-3 * rx)

def _kappa_method_c(net, ppc):
    if net.f_hz == 50:
        fc = 20

@@ -48,14 +52,17 @@

    ppc_c["bus"][conductance, GS] = y_shunt.real[0]
    ppc_c["bus"][conductance, BS] = y_shunt.imag[0]
    _calc_ybus(ppc_c)
    _calc_zbus(ppc_c)
    _calc_rx(net, ppc_c)

    if net["_options"]["inverse_y"]:
        _calc_zbus(net, ppc_c)
    else:
        # Factorization Ybus once
        ppc_c["internal"]["ybus_fact"] = factorized(ppc_c["internal"]["Ybus"])

    _calc_rx(net, ppc_c, bus=None)
    rx_equiv_c = ppc_c["bus"][:, R_EQUIV] / ppc_c["bus"][:, X_EQUIV] * fc / net.f_hz
    return _kappa(rx_equiv_c)

def _kappa(rx):
    return 1.02 + .98 * np.exp(-3 * rx)

def _kappa_method_b(net, ppc):
    topology = net._options["topology"]
    kappa_max = np.full(ppc["bus"].shape[0], 2.)

@@ -13,23 +13,37 @@


from pandapower.shortcircuit.idx_bus import R_EQUIV, X_EQUIV
from pandapower.pypower.idx_bus import BASE_KV
from pandapower.auxiliary import _clean_up
try:
    from pandapower.pf.makeYbus_numba import makeYbus
except ImportError:
    from pandapower.pypower.makeYbus import makeYbus


def _calc_rx(net, ppc):
    Zbus = ppc["internal"]["Zbus"]
def _calc_rx(net, ppc, bus):
    # Vectorized for multiple bus
    if bus is None:
        # Slice(None) is equal to select all
        bus_idx = slice(None)
    else:
        bus_idx = net._pd2ppc_lookups["bus"][bus] #bus where the short-circuit is calculated (j)

    r_fault = net["_options"]["r_fault_ohm"]
    x_fault = net["_options"]["x_fault_ohm"]
    if r_fault > 0 or x_fault > 0:
        base_r = np.square(ppc["bus"][:, BASE_KV]) / ppc["baseMVA"]
        base_r = np.square(ppc["bus"][bus_idx, BASE_KV]) / ppc["baseMVA"]
        fault_impedance = (r_fault + x_fault * 1j) / base_r
        np.fill_diagonal(Zbus, Zbus.diagonal() + fault_impedance)
    z_equiv = np.diag(Zbus)
    ppc["bus"][:, R_EQUIV] = z_equiv.real
    ppc["bus"][:, X_EQUIV] = z_equiv.imag
    else:
        fault_impedance = 0 + 0j

    if net["_options"]["inverse_y"]:
        Zbus = ppc["internal"]["Zbus"]
        z_equiv = np.diag(Zbus)[bus_idx] + fault_impedance
    else:
        z_equiv = _calc_zbus_diag(net, ppc, bus) + fault_impedance
    ppc["bus"][bus_idx, R_EQUIV] = z_equiv.real
    ppc["bus"][bus_idx, X_EQUIV] = z_equiv.imag


def _calc_ybus(ppc):
    Ybus, Yf, Yt = makeYbus(ppc["baseMVA"], ppc["bus"],  ppc["branch"])

@@ -39,12 +53,60 @@

    ppc["internal"]["Yt"] = Yt
    ppc["internal"]["Ybus"] = Ybus

def _calc_zbus(ppc):
    Ybus = ppc["internal"]["Ybus"]
    sparsity = Ybus.nnz / Ybus.shape[0]**2
    if sparsity < 0.002:
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            ppc["internal"]["Zbus"] = inv_sparse(Ybus).toarray()

def _calc_zbus(net, ppc):
    try:
        Ybus = ppc["internal"]["Ybus"]
        sparsity = Ybus.nnz / Ybus.shape[0]**2
        if sparsity < 0.002:
            with warnings.catch_warnings():
                warnings.simplefilter("ignore")
                ppc["internal"]["Zbus"] = inv_sparse(Ybus).toarray()
        else:
            ppc["internal"]["Zbus"] = inv(Ybus.toarray())
    except Exception as e:
        _clean_up(net, res=False)
        raise (e)


def _calc_zbus_diag(net, ppc, bus=None):
    ybus_fact = ppc["internal"]["ybus_fact"]
    n_bus = ppc["bus"].shape[0]

    if bus is None:
        diagZ = np.zeros(n_bus, dtype=np.complex)
        for i in range(ppc["bus"].shape[0]):
            b = np.zeros(n_bus, dtype=np.complex)
            b[i] = 1 + 0j
            diagZ[i] = ybus_fact(b)[i]
        ppc["internal"]["diagZ"] = diagZ
        return diagZ
    else:
        ppc["internal"]["Zbus"] = inv(Ybus.toarray())
        if isinstance(bus, int):
            bus = np.array([bus])
        diagZ = np.zeros(np.shape(bus)[0], dtype=np.complex)
        for ix, b in enumerate(bus):
            bus_idx = net._pd2ppc_lookups["bus"][b] #bus where the short-circuit is calculated (j)
            b = np.zeros(n_bus, dtype=np.complex)
            b[bus_idx] = 1 + 0j
            diagZ[ix] = ybus_fact(b)[bus_idx]
        return diagZ

    # if bus is None:
    #     bus = net.bus.index

    # diagZ = np.zeros(np.shape(bus)[0], dtype=np.complex)
    # ix = 0

    # # Use windows size 32 to calculate Zbus
    # while ix < np.shape(bus)[0]:
    #     ix_end = min(ix+32, np.shape(bus)[0])
    #     bus_idx = net._pd2ppc_lookups["bus"][bus[ix: ix_end]]
    #     b = np.zeros((n_bus, (ix_end-ix)), dtype=np.complex)
    #     for this_ix, this_bus_ix in enumerate(bus_idx):
    #         b[this_bus_ix, this_ix] = 1 + 0j
    #     res = ybus_fact(b)
    #     for this_ix, this_bus_ix in enumerate(bus_idx):
    #         diagZ[ix] = res[this_bus_ix, this_ix]
    #     ix += 32
    # return diagZ

@@ -577,7 +577,8 @@



def _add_sc_options(net, fault, case, lv_tol_percent, tk_s, topology, r_fault_ohm,
                    x_fault_ohm, kappa, ip, ith, branch_results, kappa_method, return_all_currents):
                    x_fault_ohm, kappa, ip, ith, branch_results, kappa_method, return_all_currents,
                    inverse_y):
    """
    creates dictionary for pf, opf and short circuit calculations from input parameters.
    """

@@ -594,7 +595,8 @@

        "ith": ith,
        "branch_results": branch_results,
        "kappa_method": kappa_method,
        "return_all_currents": return_all_currents
        "return_all_currents": return_all_currents,
        "inverse_y": inverse_y
    }
    _add_options(net, options)


Files		•	•	Coverage
pandapower	15,695	13,655	2,040	87.00%
setup.py	14	0	14	0.00%
Project Totals (162 files)	15,709	13,655	2,054	86.92%

View CI Build Download

No yaml found.

Create your codecov.yml to customize your Codecov experience

Read the docs

Interact now

The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.

SVG

The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.

SVG

Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.

SVG

Learn more about the coverage sunburst graph

112	120		logger.warning("Branch results are in beta mode and might not always be reliable, "
113	121		"especially for transformers")
114	122
	123	+	# Convert bus to numpy array for better performance
	124	+	if isinstance(bus, int):
	125	+	bus = np.array([bus])
	126	+	elif isinstance(bus, list):
	127	+	bus = np.array(bus)
	128	+
115	129		kappa = ith or ip
116	130		net["_options"] = {}
117	131		_add_ppc_options(net, calculate_voltage_angles=False, trafo_model="pi",

121	135		_add_sc_options(net, fault=fault, case=case, lv_tol_percent=lv_tol_percent, tk_s=tk_s,
122	136		topology=topology, r_fault_ohm=r_fault_ohm, kappa_method=kappa_method,
123	137		x_fault_ohm=x_fault_ohm, kappa=kappa, ip=ip, ith=ith,
124		-	branch_results=branch_results, return_all_currents=return_all_currents)
	138	+	branch_results=branch_results, return_all_currents=return_all_currents,
	139	+	inverse_y=inverse_y)
125	140		init_results(net, "sc")
126		-	if fault == "3ph":
127		-	_calc_sc(net)
128		-	if fault == "2ph":
129		-	_calc_sc(net)
130		-	if fault == "1ph":
131		-	_calc_sc_1ph(net)
	141	+	if fault in ("2ph", "3ph"):
	142	+	_calc_sc(net, bus)
	143	+	elif fault == "1ph":
	144	+	_calc_sc_1ph(net, bus)
	145	+	else:
	146	+	raise ValueError("Invalid fault %s" % fault)
132	147
133	148
134		-	def calc_single_sc(net, bus, fault="3ph", case='max', lv_tol_percent=10, check_connectivity=True):
	149	+	def calc_single_sc(net, bus, fault="3ph", case='max', lv_tol_percent=10,
	150	+	check_connectivity=True, inverse_y=True):
135	151		"""
136	152		Calculates minimal or maximal symmetrical short-circuit currents.
137	153		The calculation is based on the method of the equivalent voltage source

190	206		_add_sc_options(net, fault=fault, case=case, lv_tol_percent=lv_tol_percent, tk_s=1.,
191	207		topology="auto", r_fault_ohm=0., kappa_method="C",
192	208		x_fault_ohm=0., kappa=False, ip=False, ith=False,
193		-	branch_results=True, return_all_currents=False)
	209	+	branch_results=True, return_all_currents=False,
	210	+	inverse_y=inverse_y)
194	211		init_results(net, "sc")
195		-	if fault == "3ph" or fault == "2ph":
	212	+	if fault in ("2ph", "3ph"):
196	213		_calc_sc_single(net, bus)
197	214		elif fault == "1ph":
198	215		raise NotImplementedError("1ph short-circuits are not yet implemented")

204	221		_add_auxiliary_elements(net)
205	222		ppc, ppci = _pd2ppc(net)
206	223		_calc_ybus(ppci)
207		-	try:
208		-	_calc_zbus(ppci)
209		-	except Exception as e:
210		-	_clean_up(net, res=False)
211		-	raise (e)
212		-	_calc_rx(net, ppci)
213		-	_calc_ikss(net, ppci)
214		-	_calc_single_bus_sc(net, ppci, bus)
	224	+
	225	+	if net["_options"]["inverse_y"]:
	226	+	_calc_zbus(net, ppci)
	227	+	_calc_rx(net, ppci, bus=None)
	228	+	_calc_ikss(net, ppci, bus=None)
	229	+	_calc_single_bus_sc(net, ppci, bus)
	230	+	else:
	231	+	# Factorization Ybus once
	232	+	ppci["internal"]["ybus_fact"] = factorized(ppci["internal"]["Ybus"])
	233	+
	234	+	_calc_rx(net, ppci, bus)
	235	+	_calc_ikss(net, ppci, bus)
	236	+	_calc_single_bus_sc_no_y_inv(net, ppci, bus)
	237	+
	238	+	# Delete factorization object
	239	+	ppci["internal"].pop("ybus_fact")
	240	+
215	241		ppc = _copy_results_ppci_to_ppc(ppci, ppc, "sc")
216	242		_extract_single_results(net, ppc)
217	243		_clean_up(net)
218	244
219	245
220		-	def _calc_sc(net):
	246	+	def _calc_sc(net, bus):
221	247		_add_auxiliary_elements(net)
222	248		ppc, ppci = _pd2ppc(net)
223	249		_calc_ybus(ppci)
224		-	try:
225		-	_calc_zbus(ppci)
226		-	except Exception as e:
227		-	_clean_up(net, res=False)
228		-	raise (e)
229		-	_calc_rx(net, ppci)
230		-	_add_kappa_to_ppc(net, ppci)
231		-	_calc_ikss(net, ppci)
	250	+
	251	+	if net["_options"]["inverse_y"]:
	252	+	_calc_zbus(net, ppci)
	253	+	else:
	254	+	# Factorization Ybus once
	255	+	ppci["internal"]["ybus_fact"] = factorized(ppci["internal"]["Ybus"])
	256	+
	257	+	_calc_rx(net, ppci, bus)
	258	+
	259	+	# kappa required inverse of Zbus, which is optimized
	260	+	if net["_options"]["kappa"]:
	261	+	_add_kappa_to_ppc(net, ppci)
	262	+	_calc_ikss(net, ppci, bus)
	263	+
232	264		if net["_options"]["ip"]:
233	265		_calc_ip(net, ppci)
234	266		if net["_options"]["ith"]:
235	267		_calc_ith(net, ppci)
	268	+
236	269		if net._options["branch_results"]:
237		-	_calc_branch_currents(net, ppci)
	270	+	_calc_branch_currents(net, ppci, bus)
238	271
239	272		ppc = _copy_results_ppci_to_ppc(ppci, ppc, "sc")
240		-
241		-	if net["_options"]["return_all_currents"]:
242		-	_extract_results(net, ppc, ppc_0=None)
243		-	else:
244		-	_extract_results(net, ppc, ppc_0=None)
	273	+	_extract_results(net, ppc, ppc_0=None, bus=bus)
245	274		_clean_up(net)
246	275
	276	+	if "ybus_fact" in ppci["internal"]:
	277	+	# Delete factorization object
	278	+	ppci["internal"].pop("ybus_fact")
247	279
248		-	def _calc_sc_1ph(net):
	280	+
	281	+	def _calc_sc_1ph(net, bus):
249	282		"""
250	283		calculation method for single phase to ground short-circuit currents
251	284		"""

253	286		# pos. seq bus impedance
254	287		ppc, ppci = _pd2ppc(net)
255	288		_calc_ybus(ppci)
256		-	try:
257		-	_calc_zbus(ppci)
258		-	except Exception as e:
259		-	_clean_up(net, res=False)
260		-	raise (e)
261		-	_calc_rx(net, ppci)
262		-	_add_kappa_to_ppc(net, ppci)
	289	+
263	290		# zero seq bus impedance
264	291		ppc_0, ppci_0 = _pd2ppc_zero(net)
265	292		_calc_ybus(ppci_0)
266		-	try:
267		-	_calc_zbus(ppci_0)
268		-	except Exception as e:
269		-	_clean_up(net, res=False)
270		-	raise (e)
271		-	_calc_rx(net, ppci_0)
272		-	_calc_ikss_1ph(net, ppci, ppci_0)
	293	+
	294	+	if net["_options"]["inverse_y"]:
	295	+	_calc_zbus(net, ppci)
	296	+	_calc_zbus(net, ppci_0)
	297	+	else:
	298	+	# Factorization Ybus once
	299	+	ppci["internal"]["ybus_fact"] = factorized(ppci["internal"]["Ybus"])
	300	+	ppci_0["internal"]["ybus_fact"] = factorized(ppci_0["internal"]["Ybus"])
	301	+
	302	+	_calc_rx(net, ppci, bus=bus)
	303	+	_add_kappa_to_ppc(net, ppci)
	304	+
	305	+	_calc_rx(net, ppci_0, bus=bus)
	306	+	_calc_ikss_1ph(net, ppci, ppci_0, bus=bus)
	307	+
273	308		if net._options["branch_results"]:
274		-	_calc_branch_currents(net, ppci)
	309	+	_calc_branch_currents(net, ppci, bus=bus)
275	310		ppc_0 = _copy_results_ppci_to_ppc(ppci_0, ppc_0, "sc")
276	311		ppc = _copy_results_ppci_to_ppc(ppci, ppc, "sc")
277		-	_extract_results(net, ppc, ppc_0)
	312	+	_extract_results(net, ppc, ppc_0, bus=bus)
278	313		_clean_up(net)

10	10		import logging
11	11
12	12		logger = logging.getLogger(__name__)
13		-	# import time
	13	+
	14	+	import numpy as np
	15	+	from scipy.sparse.linalg import factorized
14	16
15	17		from pandapower.auxiliary import _clean_up, _add_ppc_options, _add_sc_options, _add_auxiliary_elements
16	18		from pandapower.pd2ppc import _pd2ppc
17	19		from pandapower.pd2ppc_zero import _pd2ppc_zero
18	20		from pandapower.results import _copy_results_ppci_to_ppc
19		-	from pandapower.shortcircuit.currents import _calc_ikss, _calc_ikss_1ph, _calc_ip, _calc_ith, _calc_branch_currents, \
20		-	_calc_single_bus_sc
	21	+	from pandapower.shortcircuit.currents import _calc_ikss, \
	22	+	_calc_ikss_1ph, _calc_ip, _calc_ith, _calc_branch_currents, \
	23	+	_calc_single_bus_sc, _calc_single_bus_sc_no_y_inv
21	24		from pandapower.shortcircuit.impedance import _calc_zbus, _calc_ybus, _calc_rx
22	25		from pandapower.shortcircuit.kappa import _add_kappa_to_ppc
23	26		from pandapower.shortcircuit.results import _extract_results, _extract_single_results

26	29
27	30		def calc_sc(net, fault="3ph", case='max', lv_tol_percent=10, topology="auto", ip=False,
28	31		ith=False, tk_s=1., kappa_method="C", r_fault_ohm=0., x_fault_ohm=0.,
29		-	branch_results=False, check_connectivity=True, return_all_currents=False):
	32	+	branch_results=False, check_connectivity=True, return_all_currents=False,
	33	+	bus=None, inverse_y=True):
30	34		"""
31	35		Calculates minimal or maximal symmetrical short-circuit currents.
32	36		The calculation is based on the method of the equivalent voltage source

85	89		return_all_currents (bool, False) applies only if branch_results=True, if True short-circuit currents for
86	90		each (branch, bus) tuple is returned otherwise only the max/min is returned
87	91
	92	+	bus (int, list, np.array, None) defines if short-circuit calculations should only be calculated for defined bus
	93	+
	94	+	inverse_y (bool, True) defines if complete inverse should be used instead of LU factorization, factorization version is in experiment which should be faster and memory efficienter
	95	+
88	96
89	97		OUTPUT:
90	98

12	12		from pandapower.pypower.idx_gen import GEN_BUS, MBASE
13	13		from pandapower.shortcircuit.idx_brch import IKSS_F, IKSS_T, IP_F, IP_T, ITH_F, ITH_T
14	14		from pandapower.shortcircuit.idx_bus import C_MIN, C_MAX, KAPPA, R_EQUIV, IKSS1, IP, ITH, X_EQUIV, IKSS2, IKCV, M
	15	+	from pandapower.shortcircuit.impedance import _calc_zbus_diag
15	16
	17	+	from pandapower.pypower.pfsoln import pfsoln as pfsoln_pypower
	18	+	from pandapower.pf.ppci_variables import _get_pf_variables_from_ppci
	19	+
	20	+	def _calc_ikss(net, ppc, bus=None):
	21	+	# Vectorized for multiple bus
	22	+	if bus is None:
	23	+	# Slice(None) is equal to : select
	24	+	bus_idx = slice(None)
	25	+	else:
	26	+	bus_idx = net._pd2ppc_lookups["bus"][bus] #bus where the short-circuit is calculated (j)
16	27
17		-	def _calc_ikss(net, ppc):
18	28		fault = net._options["fault"]
19	29		case = net._options["case"]
20		-	c = ppc["bus"][:, C_MIN] if case == "min" else ppc["bus"][:, C_MAX]
	30	+	c = ppc["bus"][bus_idx, C_MIN] if case == "min" else ppc["bus"][bus_idx, C_MAX]
21	31		ppc["internal"]["baseI"] = ppc["bus"][:, BASE_KV] * np.sqrt(3) / ppc["baseMVA"]
22		-	z_equiv = abs(ppc["bus"][:, R_EQUIV] + ppc["bus"][:, X_EQUIV] * 1j)
	32	+
	33	+	z_equiv = abs(ppc["bus"][bus_idx, R_EQUIV] + ppc["bus"][bus_idx, X_EQUIV] * 1j)
23	34		if fault == "3ph":
24		-	ppc["bus"][:, IKSS1] = c / z_equiv / ppc["bus"][:, BASE_KV] / np.sqrt(3) * ppc["baseMVA"]
	35	+	ppc["bus"][bus_idx, IKSS1] = c / z_equiv / ppc["bus"][bus_idx, BASE_KV] / np.sqrt(3) * ppc["baseMVA"]
25	36		elif fault == "2ph":
26		-	ppc["bus"][:, IKSS1] = c / z_equiv / ppc["bus"][:, BASE_KV] / 2 * ppc["baseMVA"]
	37	+	ppc["bus"][bus_idx, IKSS1] = c / z_equiv / ppc["bus"][bus_idx, BASE_KV] / 2 * ppc["baseMVA"]
27	38		_current_source_current(net, ppc)
28	39
29	40
30		-	def _calc_ikss_1ph(net, ppc, ppc_0):
31		-	fault = net._options["fault"]
	41	+	def _calc_ikss_1ph(net, ppc, ppc_0, bus=None):
	42	+	# Vectorized for multiple bus
	43	+	if bus is None:
	44	+	# Slice(None) is equal to : select
	45	+	bus_idx = slice(None)
	46	+	else:
	47	+	bus_idx = net._pd2ppc_lookups["bus"][bus] #bus where the short-circuit is calculated (j)
	48	+
32	49		case = net._options["case"]
33		-	c = ppc["bus"][:, C_MIN] if case == "min" else ppc["bus"][:, C_MAX]
	50	+	c = ppc["bus"][bus_idx, C_MIN] if case == "min" else ppc["bus"][bus_idx, C_MAX]
34	51		ppc["internal"]["baseI"] = ppc["bus"][:, BASE_KV] * np.sqrt(3) / ppc["baseMVA"]
35	52		ppc_0["internal"]["baseI"] = ppc_0["bus"][:, BASE_KV] * np.sqrt(3) / ppc_0["baseMVA"]
36		-	z_equiv = abs((ppc["bus"][:, R_EQUIV] + ppc["bus"][:, X_EQUIV] * 1j) * 2 +
37		-	(ppc_0["bus"][:, R_EQUIV] + ppc_0["bus"][:, X_EQUIV] * 1j))
38		-	if fault == "1ph":
39		-	ppc_0["bus"][:, IKSS1] = c / z_equiv / ppc_0["bus"][:, BASE_KV] * np.sqrt(3) * ppc_0["baseMVA"]
40		-	ppc["bus"][:, IKSS1] = c / z_equiv / ppc_0["bus"][:, BASE_KV] * np.sqrt(3) * ppc_0["baseMVA"]
41		-	_current_source_current(net, ppc)
	53	+
	54	+	z_equiv = abs((ppc["bus"][bus_idx, R_EQUIV] + ppc["bus"][bus_idx, X_EQUIV] * 1j) * 2 +
	55	+	(ppc_0["bus"][bus_idx, R_EQUIV] + ppc_0["bus"][bus_idx, X_EQUIV] * 1j))
	56	+
	57	+	ppc_0["bus"][bus_idx, IKSS1] = c / z_equiv / ppc_0["bus"][bus_idx, BASE_KV] * np.sqrt(3) * ppc_0["baseMVA"]
	58	+	ppc["bus"][bus_idx, IKSS1] = c / z_equiv / ppc["bus"][bus_idx, BASE_KV] * np.sqrt(3) * ppc["baseMVA"]
	59	+	_current_source_current(net, ppc)
42	60
43	61
44	62		def _current_source_current(net, ppc):

56	74		baseI = ppc["internal"]["baseI"]
57	75		sgen_buses = sgen.bus.values
58	76		sgen_buses_ppc = bus_lookup[sgen_buses]
59		-	Zbus = ppc["internal"]["Zbus"]
60	77		if not "k" in sgen:
61	78		raise ValueError("Nominal to short-circuit current has to specified in net.sgen.k")
62	79		i_sgen_pu = sgen.sn_mva.values / net.sn_mva * sgen.k.values
63	80		buses, ikcv_pu, _ = _sum_by_group(sgen_buses_ppc, i_sgen_pu, i_sgen_pu)
64	81		ppc["bus"][buses, IKCV] = ikcv_pu
65		-	ppc["bus"][:, IKSS2] = abs(1 / np.diag(Zbus) * np.dot(Zbus, ppc["bus"][:, IKCV] * -1j) / baseI)
	82	+	if net["_options"]["inverse_y"]:
	83	+	Zbus = ppc["internal"]["Zbus"]
	84	+	ppc["bus"][:, IKSS2] = abs(1 / np.diag(Zbus) * np.dot(Zbus, ppc["bus"][:, IKCV] * -1j) / baseI)
	85	+	else:
	86	+	ybus_fact = ppc["internal"]["ybus_fact"]
	87	+	diagZ = _calc_zbus_diag(net, ppc)
	88	+	ppc["bus"][:, IKSS2] = abs(ybus_fact(ppc["bus"][:, IKCV] * -1j) / diagZ / baseI)
66	89		ppc["bus"][buses, IKCV] /= baseI[buses]
67	90
68	91

83	106		ppc["bus"][:, ITH] = ith
84	107
85	108
86		-	def _calc_branch_currents(net, ppc):
87		-	case = net._options["case"]
88		-	Zbus = ppc["internal"]["Zbus"]
89		-	Yf = ppc["internal"]["Yf"]
90		-	Yt = ppc["internal"]["Yt"]
91		-	baseI = ppc["internal"]["baseI"]
92		-	n = ppc["bus"].shape[0]
93		-	fb = np.real(ppc["branch"][:, 0]).astype(int)
94		-	tb = np.real(ppc["branch"][:, 1]).astype(int)
95		-	minmax = np.nanmin if case == "min" else np.nanmax
96		-
97		-	# calculate voltage source branch current
98		-	V_ikss = (ppc["bus"][:, IKSS1] * baseI) * Zbus
99		-	ikss1_all_f = np.conj(Yf.dot(V_ikss))
100		-	ikss1_all_t = np.conj(Yt.dot(V_ikss))
101		-	ikss1_all_f[abs(ikss1_all_f) < 1e-10] = 0.
102		-	ikss1_all_t[abs(ikss1_all_t) < 1e-10] = 0.
103		-
104		-	# add current source branch current if there is one
105		-	current_sources = any(ppc["bus"][:, IKCV]) > 0
106		-	if current_sources:
107		-	current = np.tile(-ppc["bus"][:, IKCV], (n, 1))
108		-	np.fill_diagonal(current, current.diagonal() + ppc["bus"][:, IKSS2])
109		-	V = np.dot((current * baseI), Zbus).T
110		-	fb = np.real(ppc["branch"][:, 0]).astype(int)
111		-	tb = np.real(ppc["branch"][:, 1]).astype(int)
112		-	ikss2_all_f = np.conj(Yf.dot(V))
113		-	ikss2_all_t = np.conj(Yt.dot(V))
114		-	ikss_all_f = abs(ikss1_all_f + ikss2_all_f)
115		-	ikss_all_t = abs(ikss1_all_t + ikss2_all_t)
116		-	else:
117		-	ikss_all_f = abs(ikss1_all_f)
118		-	ikss_all_t = abs(ikss1_all_t)
119		-
120		-	if net._options["return_all_currents"]:
121		-	ppc["internal"]["branch_ikss_f"] = ikss_all_f / baseI[fb, None]
122		-	ppc["internal"]["branch_ikss_t"] = ikss_all_t / baseI[tb, None]
123		-	else:
124		-	ikss_all_f[ikss_all_f < 1e-10] = np.nan
125		-	ikss_all_t[ikss_all_t < 1e-10] = np.nan
126		-	ppc["branch"][:, IKSS_F] = minmax(ikss_all_f, axis=1) / baseI[fb]
127		-	ppc["branch"][:, IKSS_T] = minmax(ikss_all_t, axis=1) / baseI[tb]
128		-
129		-	if net._options["ip"]:
130		-	kappa = ppc["bus"][:, KAPPA]
131		-	if current_sources:
132		-	ip_all_f = np.sqrt(2) * (ikss1_all_f * kappa + ikss2_all_f)
133		-	ip_all_t = np.sqrt(2) * (ikss1_all_t * kappa + ikss2_all_t)
134		-	else:
135		-	ip_all_f = np.sqrt(2) * ikss1_all_f * kappa
136		-	ip_all_t = np.sqrt(2) * ikss1_all_t * kappa
137		-
138		-	if net._options["return_all_currents"]:
139		-	ppc["internal"]["branch_ip_f"] = abs(ip_all_f) / baseI[fb, None]
140		-	ppc["internal"]["branch_ip_t"] = abs(ip_all_t) / baseI[tb, None]
141		-	else:
142		-	ip_all_f[abs(ip_all_f) < 1e-10] = np.nan
143		-	ip_all_t[abs(ip_all_t) < 1e-10] = np.nan
144		-	ppc["branch"][:, IP_F] = minmax(abs(ip_all_f), axis=1) / baseI[fb]
145		-	ppc["branch"][:, IP_T] = minmax(abs(ip_all_t), axis=1) / baseI[tb]
146		-
147		-	if net._options["ith"]:
148		-	n = 1
149		-	m = ppc["bus"][:, M]
150		-	ith_all_f = ikss_all_f * np.sqrt(m + n)
151		-	ith_all_t = ikss_all_t * np.sqrt(m + n)
152		-
153		-	if net._options["return_all_currents"]:
154		-	ppc["internal"]["branch_ith_f"] = ith_all_f / baseI[fb, None]
155		-	ppc["internal"]["branch_ith_t"] = ith_all_t / baseI[tb, None]
156		-	else:
157		-	ppc["branch"][:, ITH_F] = minmax(ith_all_f, axis=1) / baseI[fb]
158		-	ppc["branch"][:, ITH_T] = minmax(ith_all_t, axis=1) / baseI[fb]
159		-
160		-
161		-
162		-
163		-
164	109		def _calc_ib_generator(net, ppci):
165	110		Zbus = ppci["internal"]["Zbus"]
166	111		baseI = ppci["internal"]["baseI"]

230	175
231	176
232	177		def _calc_single_bus_sc(net, ppc, bus):
233		-	case = net._options["case"]
	178	+	# case = net._options["case"]
234	179		bus_idx = net._pd2ppc_lookups["bus"][bus]
235	180		n = ppc["bus"].shape[0]
236	181		Zbus = ppc["internal"]["Zbus"]

239	184		baseI = ppc["internal"]["baseI"]
240	185		# fb = np.real(ppc["branch"][:, 0]).astype(int)
241	186		# tb = np.real(ppc["branch"][:, 1]).astype(int)
242		-	c = ppc["bus"][:, C_MIN] if case == "min" else ppc["bus"][:, C_MAX]
	187	+	# c = ppc["bus"][:, C_MIN] if case == "min" else ppc["bus"][:, C_MAX]
243	188
244	189		# calculate voltage source branch current
245	190		V_ikss = (ppc["bus"][:, IKSS1] * baseI) * Zbus

257	202		# ppc["branch"][:, IKSS_T] = abs(ikss_all_t[:, bus_idx] / baseI[tb])
258	203		calc_branch_results(net, ppc, V)
259	204
	205	+	def _calc_single_bus_sc_no_y_inv(net, ppc, bus):
	206	+	# Vectorized for multiple bus
	207	+	if bus is None:
	208	+	# Slice(None) is equal to : select
	209	+	bus_idx = slice(None)
	210	+	else:
	211	+	bus_idx = net._pd2ppc_lookups["bus"][bus] #bus where the short-circuit is calculated (j)
260	212
261		-	from pandapower.pypower.pfsoln import pfsoln as pfsoln_pypower
262		-	from pandapower.pf.ppci_variables import _get_pf_variables_from_ppci
	213	+	ybus = ppc["internal"]["Ybus"]
	214	+	ybus_fact = ppc["internal"]["ybus_fact"]
	215	+	# case = net._options["case"]
	216	+	baseI = ppc["internal"]["baseI"]
	217	+	# vqj = ppc["bus"][:, C_MIN] if case == "min" else ppc["bus"][:, C_MAX] #this is the source voltage in per unit (VQj)
	218	+
	219	+	# Solve Ikss from voltage source
	220	+	n_bus = ybus.shape[0]
	221	+
	222	+	# ybus_sub_mask = (np.arange(ybus.shape[0]) != bus_idx)
	223	+	# V_ikss = np.zeros(n_bus, dtype=np.complex)
	224	+	# V_ikss[bus_idx] = vqj[bus_idx]
	225	+
	226	+	# # Solve Ax = b
	227	+	# b = np.zeros(n_bus-1, dtype=np.complex) -\
	228	+	# (ybus[:, ~ybus_sub_mask].toarray())[ybus_sub_mask].ravel() * V_ikss[bus_idx]
	229	+	# ybus_sub = ybus[ybus_sub_mask, :][:, ybus_sub_mask]
	230	+	# x = spsolve(ybus_sub, b)
	231	+
	232	+	# V_ikss[ybus_sub_mask] = x
	233	+	# I_ikss = np.zeros(n_bus, dtype=np.complex)
	234	+	# I_ikss[bus_idx] = np.dot(ybus[bus_idx, :].toarray(), V_ikss)
	235	+	# V = V_ikss
	236	+
	237	+	# Version 2
	238	+	I_ikss = np.zeros(n_bus, dtype=np.complex)
	239	+	I_ikss[bus_idx] = ppc["bus"][bus_idx, IKSS1]
	240	+	V_ikss = ybus_fact(I_ikss * baseI)
	241	+	V = V_ikss
	242	+
	243	+	#TODO include current sources
	244	+	current_sources = any(ppc["bus"][:, IKCV]) > 0
	245	+	if current_sources:
	246	+	current = -ppc["bus"][:, IKCV]
	247	+	current[bus_idx] += ppc["bus"][bus_idx, IKSS2]
	248	+	V_source = ybus_fact(current)
	249	+	V += V_source
	250	+
	251	+	calc_branch_results(net, ppc, V)
263	252
264	253
265	254		def calc_branch_results(net, ppci, V):
266	255		Ybus = ppci["internal"]["Ybus"]
267	256		Yf = ppci["internal"]["Yf"]
268	257		Yt = ppci["internal"]["Yt"]
269		-	baseMVA, bus, gen, branch, ref, _, pq, _, _, V0, ref_gens = _get_pf_variables_from_ppci(ppci)
	258	+	baseMVA, bus, gen, branch, ref, _, _, _, _, _, ref_gens = _get_pf_variables_from_ppci(ppci)
270	259		bus, gen, branch = pfsoln_pypower(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, ref_gens)
271	260		ppci["bus"], ppci["gen"], ppci["branch"] = bus, gen, branch
	261	+
	262	+
	263	+	def _calc_branch_currents(net, ppc, bus):
	264	+	# Vectorized for multiple bus
	265	+	if bus is None:
	266	+	# Slice(None) is equal to select all
	267	+	bus = net.bus.index
	268	+
	269	+	bus_idx = net._pd2ppc_lookups["bus"][bus]
	270	+	# Select only in service bus for sc calculation
	271	+	bus_idx = bus_idx[bus_idx < ppc['bus'].shape[0]]
	272	+	n_sc_bus = np.shape(bus_idx)[0]
	273	+
	274	+	case = net._options["case"]
	275	+
	276	+	Yf = ppc["internal"]["Yf"]
	277	+	Yt = ppc["internal"]["Yt"]
	278	+	baseI = ppc["internal"]["baseI"]
	279	+	n_bus = ppc["bus"].shape[0]
	280	+	fb = np.real(ppc["branch"][:, 0]).astype(int)
	281	+	tb = np.real(ppc["branch"][:, 1]).astype(int)
	282	+	minmax = np.nanmin if case == "min" else np.nanmax
	283	+
	284	+	# calculate voltage source branch current
	285	+	if net["_options"]["inverse_y"]:
	286	+	Zbus = ppc["internal"]["Zbus"]
	287	+	V_ikss = (ppc["bus"][:, IKSS1] * baseI) * Zbus
	288	+	V_ikss = V_ikss[:, bus_idx]
	289	+	else:
	290	+	ybus_fact = ppc["internal"]["ybus_fact"]
	291	+	V_ikss = np.zeros((n_bus, n_sc_bus), dtype=np.complex)
	292	+	for ix, b in enumerate(bus_idx):
	293	+	ikss = np.zeros(n_bus, dtype=np.complex)
	294	+	ikss[b] = ppc["bus"][b, IKSS1] * baseI[b]
	295	+	V_ikss[:, ix] = ybus_fact(ikss)
	296	+
	297	+	ikss1_all_f = np.conj(Yf.dot(V_ikss))
	298	+	ikss1_all_t = np.conj(Yt.dot(V_ikss))
	299	+	ikss1_all_f[abs(ikss1_all_f) < 1e-10] = 0.
	300	+	ikss1_all_t[abs(ikss1_all_t) < 1e-10] = 0.
	301	+
	302	+	# add current source branch current if there is one
	303	+	current_sources = any(ppc["bus"][:, IKCV]) > 0
	304	+	if current_sources:
	305	+	current = np.tile(-ppc["bus"][:, IKCV], (n_sc_bus, 1))
	306	+	for ix, b in enumerate(bus_idx):
	307	+	current[ix, b] += ppc["bus"][b, IKSS2]
	308	+
	309	+	# calculate voltage source branch current
	310	+	if net["_options"]["inverse_y"]:
	311	+	Zbus = ppc["internal"]["Zbus"]
	312	+	V = np.dot((current * baseI), Zbus).T
	313	+	else:
	314	+	ybus_fact = ppc["internal"]["ybus_fact"]
	315	+	V = np.zeros((n_bus, n_sc_bus), dtype=np.complex)
	316	+	for ix, b in enumerate(bus_idx):
	317	+	V[:, ix] = ybus_fact(current[ix, :] * baseI[b])
	318	+
	319	+	fb = np.real(ppc["branch"][:, 0]).astype(int)
	320	+	tb = np.real(ppc["branch"][:, 1]).astype(int)
	321	+	ikss2_all_f = np.conj(Yf.dot(V))
	322	+	ikss2_all_t = np.conj(Yt.dot(V))
	323	+
	324	+	ikss_all_f = abs(ikss1_all_f + ikss2_all_f)
	325	+	ikss_all_t = abs(ikss1_all_t + ikss2_all_t)
	326	+	else:
	327	+	ikss_all_f = abs(ikss1_all_f)
	328	+	ikss_all_t = abs(ikss1_all_t)
	329	+
	330	+	if net._options["return_all_currents"]:
	331	+	ppc["internal"]["branch_ikss_f"] = ikss_all_f / baseI[fb, None]
	332	+	ppc["internal"]["branch_ikss_t"] = ikss_all_t / baseI[tb, None]
	333	+	else:
	334	+	ikss_all_f[abs(ikss_all_f) < 1e-10] = np.nan
	335	+	ikss_all_t[abs(ikss_all_t) < 1e-10] = np.nan
	336	+	ppc["branch"][:, IKSS_F] = np.nan_to_num(minmax(ikss_all_f, axis=1) / baseI[fb])
	337	+	ppc["branch"][:, IKSS_T] = np.nan_to_num(minmax(ikss_all_t, axis=1) / baseI[tb])
	338	+
	339	+	if net._options["ip"]:
	340	+	kappa = ppc["bus"][:, KAPPA]
	341	+	if current_sources:
	342	+	ip_all_f = np.sqrt(2) * (ikss1_all_f * kappa[bus_idx] + ikss2_all_f)
	343	+	ip_all_t = np.sqrt(2) * (ikss1_all_t * kappa[bus_idx] + ikss2_all_t)
	344	+	else:
	345	+	ip_all_f = np.sqrt(2) * ikss1_all_f * kappa[bus_idx]
	346	+	ip_all_t = np.sqrt(2) * ikss1_all_t * kappa[bus_idx]
	347	+
	348	+	if net._options["return_all_currents"]:
	349	+	ppc["internal"]["branch_ip_f"] = abs(ip_all_f) / baseI[fb, None]
	350	+	ppc["internal"]["branch_ip_t"] = abs(ip_all_t) / baseI[tb, None]
	351	+	else:
	352	+	ip_all_f[abs(ip_all_f) < 1e-10] = np.nan
	353	+	ip_all_t[abs(ip_all_t) < 1e-10] = np.nan
	354	+	ppc["branch"][:, IP_F] = np.nan_to_num(minmax(abs(ip_all_f), axis=1) / baseI[fb])
	355	+	ppc["branch"][:, IP_T] = np.nan_to_num(minmax(abs(ip_all_t), axis=1) / baseI[tb])
	356	+
	357	+	if net._options["ith"]:
	358	+	n = 1
	359	+	m = ppc["bus"][bus_idx, M]
	360	+	ith_all_f = ikss_all_f * np.sqrt(m + n)
	361	+	ith_all_t = ikss_all_t * np.sqrt(m + n)
	362	+
	363	+	if net._options["return_all_currents"]:
	364	+	ppc["internal"]["branch_ith_f"] = ith_all_f / baseI[fb, None]
	365	+	ppc["internal"]["branch_ith_t"] = ith_all_t / baseI[tb, None]
	366	+	else:
	367	+	ppc["branch"][:, ITH_F] = np.nan_to_num(minmax(ith_all_f, axis=1) / baseI[fb])
	368	+	ppc["branch"][:, ITH_T] = np.nan_to_num(minmax(ith_all_t, axis=1) / baseI[fb])

189	189
190	190		def _ppci_internal_to_ppc(result, ppc):
191	191		for key, value in result["internal"].items():
	192	+	# Only for sc calculation
192	193		# if branch current matrices have been stored they need to include out of service elements
193		-	if key in ["branch_ikss_f", "branch_ikss_t", "branch_ip_f", "branch_ip_t", "branch_ith_f", "branch_ith_t"]:
194		-	n_buses = np.shape(ppc['bus'])[0]
	194	+	if key in ["branch_ikss_f", "branch_ikss_t",
	195	+	"branch_ip_f", "branch_ip_t",
	196	+	"branch_ith_f", "branch_ith_t"]:
	197	+
	198	+	# n_buses = np.shape(ppc['bus'])[0]
195	199		n_branches = np.shape(ppc['branch'])[0]
196		-	n_rows_result = np.shape(result['bus'])[0]
197		-	update_matrix = np.empty((n_branches, n_buses)) * np.nan
198		-	update_matrix[result["internal"]['branch_is'], :n_rows_result] = result["internal"][key]
199		-	ppc['internal'][key] = np.copy(update_matrix)
	200	+	# n_rows_result = np.shape(result['bus'])[0]
	201	+	# update_matrix = np.empty((n_branches, n_buses)) * np.nan
	202	+	# update_matrix[result["internal"]['branch_is'], :n_rows_result] = result["internal"][key]
	203	+
	204	+	# # To select only required buses and pad one column of nan value for oos bus
	205	+	update_matrix = np.empty((n_branches, value.shape[1]+1)) * 0.0
	206	+	update_matrix[result["internal"]['branch_is'],
	207	+	:value.shape[1]] = result["internal"][key]
	208	+	ppc['internal'][key] = update_matrix
200	209		else:
201	210		ppc["internal"][key] = value
202	211

594	595		"ith": ith,
595	596		"branch_results": branch_results,
596	597		"kappa_method": kappa_method,
597		-	"return_all_currents": return_all_currents
	598	+	"return_all_currents": return_all_currents,
	599	+	"inverse_y": inverse_y
598	600		}
599	601		_add_options(net, options)
600	602

Navigation	Overlay
`t` Navigate files	`h` Toggle hits
`y` Change url to tip of branch	`m` Toggle misses
`b / v` Jump to prev/next hit line	`p` Toggle partial
`z / x` Jump to prev/next missed or partial line	`1..9` Toggle flags
`shift + o` Open current page in GitHub	`a` Toggle all on
`/ or ?` Show keyboard shortcuts dialog	`c` Toggle context lines or commits

9	9
10	10		from pandapower.shortcircuit.idx_brch import IKSS_F, IKSS_T, IP_F, IP_T, ITH_F, ITH_T
11	11		from pandapower.shortcircuit.idx_bus import IKSS1, IP, ITH, IKSS2
12		-	from pandapower.pypower.idx_bus import VM, VA
	12	+	from pandapower.pypower.idx_bus import VM, VA, BUS_TYPE
13	13		from pandapower.results_bus import _get_bus_idx, _set_buses_out_of_service
14	14		from pandapower.results import _get_aranged_lookup, _get_branch_results
15	15		from pandapower.shortcircuit.idx_bus import C_MIN, C_MAX
16	16
	17	+	def _get_bus_ppc_idx_for_br_all_results(net, ppc, bus):
	18	+	bus_lookup = net._pd2ppc_lookups["bus"]
	19	+	if bus is None:
	20	+	bus = net.bus.index
	21	+
	22	+	ppc_index = np.arange(np.shape(bus)[0])
	23	+	ppc_index[ppc["bus"][bus_lookup[ppc_index], BUS_TYPE] == 4] = -1
	24	+	return bus, ppc_index
17	25
18		-	def _extract_results(net, ppc, ppc_0):
19		-	_get_bus_results(net, ppc, ppc_0)
	26	+	def _extract_results(net, ppc, ppc_0, bus):
	27	+	_get_bus_results(net, ppc, ppc_0, bus)
20	28		if net._options["branch_results"]:
21	29		if net._options['return_all_currents']:
22		-	_get_line_all_results(net, ppc)
23		-	_get_trafo_all_results(net, ppc)
24		-	_get_trafo3w_all_results(net, ppc)
	30	+	_get_line_all_results(net, ppc, bus)
	31	+	_get_trafo_all_results(net, ppc, bus)
	32	+	_get_trafo3w_all_results(net, ppc, bus)
25	33		else:
26	34		_get_line_results(net, ppc)
27	35		_get_trafo_results(net, ppc)
28	36		_get_trafo3w_results(net, ppc)
29	37
30		-
31	38		def _extract_single_results(net, ppc):
32	39		for element in ["line", "trafo"]:
33	40		net["res_%s_sc"%element] = pd.DataFrame(np.nan, index=net[element].index,

51	58		net["res_bus_sc"]["va_degree"] = ppc["bus"][bus_idx, VA]
52	59
53	60
54		-	def _get_bus_results(net, ppc, ppc_0):
	61	+	def _get_bus_results(net, ppc, ppc_0, bus):
	62	+	if bus is None:
	63	+	bus = slice(None)
	64	+
55	65		bus_lookup = net._pd2ppc_lookups["bus"]
56	66		ppc_index = bus_lookup[net.bus.index]
57	67

64	74		if net._options["ith"]:
65	75		net.res_bus_sc["ith_ka"] = ppc["bus"][ppc_index, ITH]
66	76
	77	+	net.res_bus_sc = net.res_bus_sc.loc[bus, :]
	78	+
67	79
68	80		def _get_line_results(net, ppc):
69	81		branch_lookup = net._pd2ppc_lookups["branch"]

78	90		net.res_line_sc["ith_ka"] = minmax(ppc["branch"][f:t, [ITH_F, ITH_T]].real, axis=1)
79	91
80	92
81		-	def _get_line_all_results(net, ppc):
	93	+	def _get_line_all_results(net, ppc, bus):
82	94		case = net._options["case"]
83		-	bus_lookup = net._pd2ppc_lookups["bus"]
84		-	ppc_index = bus_lookup[net.bus.index]
	95	+
	96	+	bus, ppc_index = _get_bus_ppc_idx_for_br_all_results(net, ppc, bus)
85	97		branch_lookup = net._pd2ppc_lookups["branch"]
86	98
87		-	multindex = pd.MultiIndex.from_product([net.res_line_sc.index, net.bus.index], names=['line','bus'])
	99	+	multindex = pd.MultiIndex.from_product([net.res_line_sc.index, bus], names=['line','bus'])
88	100		net.res_line_sc = net.res_line_sc.reindex(multindex)
89	101
90	102		if "line" in branch_lookup:

100	112		net.res_line_sc["ith_ka"] = minmax(ppc["internal"]["branch_ith_f"][f:t, ppc_index].real.reshape(-1, 1),
101	113		ppc["internal"]["branch_ith_t"][f:t, ppc_index].real.reshape(-1, 1))
102	114
103		-
104	115		def _get_trafo_results(net, ppc):
105	116		branch_lookup = net._pd2ppc_lookups["branch"]
106	117		if "trafo" in branch_lookup:

109	120		net.res_trafo_sc["ikss_lv_ka"] = ppc["branch"][f:t, IKSS_T].real
110	121
111	122
112		-	def _get_trafo_all_results(net, ppc):
113		-	bus_lookup = net._pd2ppc_lookups["bus"]
114		-	ppc_index = bus_lookup[net.bus.index]
	123	+	def _get_trafo_all_results(net, ppc, bus):
	124	+	bus, ppc_index = _get_bus_ppc_idx_for_br_all_results(net, ppc, bus)
115	125		branch_lookup = net._pd2ppc_lookups["branch"]
116	126
117		-	multindex = pd.MultiIndex.from_product([net.res_trafo_sc.index, net.bus.index], names=['trafo', 'bus'])
	127	+	multindex = pd.MultiIndex.from_product([net.res_trafo_sc.index, bus], names=['trafo', 'bus'])
118	128		net.res_trafo_sc = net.res_trafo_sc.reindex(multindex)
119	129
120	130		if "trafo" in branch_lookup:

135	145		net.res_trafo3w_sc["ikss_lv_ka"] = ppc["branch"][mv:lv, IKSS_T].real
136	146
137	147
138		-	def _get_trafo3w_all_results(net, ppc):
139		-	bus_lookup = net._pd2ppc_lookups["bus"]
140		-	ppc_index = bus_lookup[net.bus.index]
	148	+	def _get_trafo3w_all_results(net, ppc, bus):
	149	+	bus, ppc_index = _get_bus_ppc_idx_for_br_all_results(net, ppc, bus)
141	150		branch_lookup = net._pd2ppc_lookups["branch"]
142	151
143		-	multindex = pd.MultiIndex.from_product([net.res_trafo3w_sc.index, net.bus.index], names=['trafo3w', 'bus'])
	152	+	multindex = pd.MultiIndex.from_product([net.res_trafo3w_sc.index, bus], names=['trafo3w', 'bus'])
144	153		net.res_trafo3w_sc = net.res_trafo3w_sc.reindex(multindex)
145	154
146	155		if "trafo3w" in branch_lookup:

7	7		import copy
8	8		import networkx as nx
9	9		import numpy as np
	10	+	from scipy.sparse.linalg import factorized
10	11
11	12		from pandapower.pypower.idx_brch import F_BUS, T_BUS, BR_R, BR_X
12	13		from pandapower.pypower.idx_bus import BUS_I, GS, BS, BASE_KV

29	30		raise ValueError("Unknown kappa method %s - specify B or C"%kappa_method)
30	31		ppc["bus"][:, KAPPA] = kappa
31	32
	33	+	def _kappa(rx):
	34	+	return 1.02 + .98 * np.exp(-3 * rx)
	35	+
32	36		def _kappa_method_c(net, ppc):
33	37		if net.f_hz == 50:
34	38		fc = 20

48	52		ppc_c["bus"][conductance, GS] = y_shunt.real[0]
49	53		ppc_c["bus"][conductance, BS] = y_shunt.imag[0]
50	54		_calc_ybus(ppc_c)
51		-	_calc_zbus(ppc_c)
52		-	_calc_rx(net, ppc_c)
	55	+
	56	+	if net["_options"]["inverse_y"]:
	57	+	_calc_zbus(net, ppc_c)
	58	+	else:
	59	+	# Factorization Ybus once
	60	+	ppc_c["internal"]["ybus_fact"] = factorized(ppc_c["internal"]["Ybus"])
	61	+
	62	+	_calc_rx(net, ppc_c, bus=None)
53	63		rx_equiv_c = ppc_c["bus"][:, R_EQUIV] / ppc_c["bus"][:, X_EQUIV] * fc / net.f_hz
54	64		return _kappa(rx_equiv_c)
55	65
56		-	def _kappa(rx):
57		-	return 1.02 + .98 * np.exp(-3 * rx)
58		-
59	66		def _kappa_method_b(net, ppc):
60	67		topology = net._options["topology"]
61	68		kappa_max = np.full(ppc["bus"].shape[0], 2.)

39	53		ppc["internal"]["Yt"] = Yt
40	54		ppc["internal"]["Ybus"] = Ybus
41	55
42		-	def _calc_zbus(ppc):
43		-	Ybus = ppc["internal"]["Ybus"]
44		-	sparsity = Ybus.nnz / Ybus.shape[0]**2
45		-	if sparsity < 0.002:
46		-	with warnings.catch_warnings():
47		-	warnings.simplefilter("ignore")
48		-	ppc["internal"]["Zbus"] = inv_sparse(Ybus).toarray()
	56	+
	57	+	def _calc_zbus(net, ppc):
	58	+	try:
	59	+	Ybus = ppc["internal"]["Ybus"]
	60	+	sparsity = Ybus.nnz / Ybus.shape[0]**2
	61	+	if sparsity < 0.002:
	62	+	with warnings.catch_warnings():
	63	+	warnings.simplefilter("ignore")
	64	+	ppc["internal"]["Zbus"] = inv_sparse(Ybus).toarray()
	65	+	else:
	66	+	ppc["internal"]["Zbus"] = inv(Ybus.toarray())
	67	+	except Exception as e:
	68	+	_clean_up(net, res=False)
	69	+	raise (e)
	70	+
	71	+
	72	+	def _calc_zbus_diag(net, ppc, bus=None):
	73	+	ybus_fact = ppc["internal"]["ybus_fact"]
	74	+	n_bus = ppc["bus"].shape[0]
	75	+
	76	+	if bus is None:
	77	+	diagZ = np.zeros(n_bus, dtype=np.complex)
	78	+	for i in range(ppc["bus"].shape[0]):
	79	+	b = np.zeros(n_bus, dtype=np.complex)
	80	+	b[i] = 1 + 0j
	81	+	diagZ[i] = ybus_fact(b)[i]
	82	+	ppc["internal"]["diagZ"] = diagZ
	83	+	return diagZ
49	84		else:
50		-	ppc["internal"]["Zbus"] = inv(Ybus.toarray())
	85	+	if isinstance(bus, int):
	86	+	bus = np.array([bus])
	87	+	diagZ = np.zeros(np.shape(bus)[0], dtype=np.complex)
	88	+	for ix, b in enumerate(bus):
	89	+	bus_idx = net._pd2ppc_lookups["bus"][b] #bus where the short-circuit is calculated (j)
	90	+	b = np.zeros(n_bus, dtype=np.complex)
	91	+	b[bus_idx] = 1 + 0j
	92	+	diagZ[ix] = ybus_fact(b)[bus_idx]
	93	+	return diagZ
	94	+
	95	+	# if bus is None:
	96	+	# bus = net.bus.index
	97	+
	98	+	# diagZ = np.zeros(np.shape(bus)[0], dtype=np.complex)
	99	+	# ix = 0
	100	+
	101	+	# # Use windows size 32 to calculate Zbus
	102	+	# while ix < np.shape(bus)[0]:
	103	+	# ix_end = min(ix+32, np.shape(bus)[0])
	104	+	# bus_idx = net._pd2ppc_lookups["bus"][bus[ix: ix_end]]
	105	+	# b = np.zeros((n_bus, (ix_end-ix)), dtype=np.complex)
	106	+	# for this_ix, this_bus_ix in enumerate(bus_idx):
	107	+	# b[this_bus_ix, this_ix] = 1 + 0j
	108	+	# res = ybus_fact(b)
	109	+	# for this_ix, this_bus_ix in enumerate(bus_idx):
	110	+	# diagZ[ix] = res[this_bus_ix, this_ix]
	111	+	# ix += 32
	112	+	# return diagZ

577	577
578	578
579	579		def _add_sc_options(net, fault, case, lv_tol_percent, tk_s, topology, r_fault_ohm,
580		-	x_fault_ohm, kappa, ip, ith, branch_results, kappa_method, return_all_currents):
	580	+	x_fault_ohm, kappa, ip, ith, branch_results, kappa_method, return_all_currents,
	581	+	inverse_y):
581	582		"""
582	583		creates dictionary for pf, opf and short circuit calculations from input parameters.
583	584		"""