Commit ⋅ e2nIEE/pandapower

#from pandapower.pd2ppc import _ppc2ppci, _init_ppc
from pandapower.pypower.idx_brch import BR_B, BR_R, BR_X, F_BUS, T_BUS, branch_cols, BR_STATUS, SHIFT, TAP
from pandapower.pypower.idx_bus import BASE_KV, BS, GS
from pandapower.build_branch import _calc_tap_from_dataframe, _transformer_correction_factor, _calc_nominal_ratio_from_dataframe
from pandapower.build_branch import _calc_tap_from_dataframe, _transformer_correction_factor, _calc_nominal_ratio_from_dataframe,\
     get_trafo_values, _trafo_df_from_trafo3w, _calc_branch_values_from_trafo_df
from pandapower.build_branch import _switch_branches, _branches_with_oos_buses, _initialize_branch_lookup, _end_temperature_correction_factor

def _pd2ppc_zero(net, sequence=0):

@@ -72,10 +73,14 @@

        branch_sc.fill(np.nan)
        ppc["branch"] = np.hstack((ppc["branch"], branch_sc ))
    ppc["branch"][:, :13] = np.array([0, 0, 0, 0, 0, 250, 250, 250, 1, 0, 1, -360, 360])

    _add_line_sc_impedance_zero(net, ppc)
    _add_trafo_sc_impedance_zero(net, ppc)
    if "trafo3w" in lookup:
        raise NotImplementedError("Three winding transformers are not implemented for unbalanced calculations")
    if mode == "sc":
        _add_trafo3w_sc_impedance_zero(net, ppc)
    else:
        if "trafo3w" in lookup:
            raise NotImplementedError("Three winding transformers are not implemented for unbalanced calculations")


def _add_trafo_sc_impedance_zero(net, ppc, trafo_df=None):

@@ -116,7 +121,7 @@

        # Just put pos seq parameter if zero seq parameter is zero
        if not "vkr0_percent" in trafos:
            raise ValueError("Real part of short circuit voltage Vk0(Real) needs to be specified for transformer \
                             modelling \n Try : net.trafo[\"vkr0_percent\"] = net.trafo[\"vkr_percent\"]" )        
                             modelling \n Try : net.trafo[\"vkr0_percent\"] = net.trafo[\"vkr_percent\"]" )
        vkr0_percent = trafos["vkr0_percent"].values.astype(float) if \
            trafos["vkr0_percent"].values.astype(float).all() != 0. else \
            trafos["vkr_percent"].values.astype(float)

@@ -129,15 +134,15 @@

            #        and mag0_percent = 10 ... 100  Zm0/ Zsc0
            # --pg 50 DigSilent Power Factory Transformer manual
            raise ValueError("Magnetizing impedance to vk0 ratio needs to be specified for transformer \
                             modelling  \n Try : net.trafo[\"mag0_percent\"] = 100" )           
                             modelling  \n Try : net.trafo[\"mag0_percent\"] = 100" )
        mag0_ratio = trafos.mag0_percent.values.astype(float)
        if not "mag0_rx" in trafos:
            raise ValueError("Magnetizing impedance R/X ratio needs to be specified for transformer \
                             modelling \n Try : net.trafo[\"mag0_rx\"] = 0 " )        
                             modelling \n Try : net.trafo[\"mag0_rx\"] = 0 " )
        mag0_rx = trafos["mag0_rx"].values.astype(float)
        if not "si0_hv_partial" in trafos:
            raise ValueError("Zero sequence short circuit impedance partition towards HV side needs to be specified for transformer \
                             modelling \n Try : net.trafo[\"si0_hv_partial\"] = 0.9 " )          
                             modelling \n Try : net.trafo[\"si0_hv_partial\"] = 0.9 " )
        si0_hv_partial = trafos.si0_hv_partial.values.astype(float)
        parallel = trafos.parallel.values.astype(float)
        in_service = trafos["in_service"].astype(int)

@@ -198,7 +203,7 @@

198	203		z1 = si0_hv_partial * z0_k
199	204		z2 = (1 - si0_hv_partial) * z0_k
200	205		z3 = z0_mag
201		-	z_temp = z1 * z2 + z2 * z3 + z1 * z3
	206	+	z_temp = z1 * z2 + z2 * z3 + z1 * z3
202	207		za = z_temp / z2
203	208		# za = z_temp / (z2+z3)
204	209		zb = z_temp / z1

@@ -208,7 +213,7 @@

        YAB = 1 / zc.astype(complex)
        YAN = 1 / za.astype(complex)
        YBN = 1 / zb.astype(complex)
        

#        YAB_AN = (zc + za) /(zc * za).astype(complex)  # Series conn YAB and YAN
#        YAB_BN = (zc + zb) / (zc * zb).astype(complex)  # Series conn YAB and YBN


@@ -243,7 +248,7 @@


            gs_all = np.hstack([gs_all, y.real * in_service])
            bs_all = np.hstack([bs_all, y.imag * in_service])
            


        elif vector_group == "YNyn":
            ppc["branch"][ppc_idx, BR_STATUS] = in_service

@@ -296,6 +301,7 @@

    ppc["bus"][buses, BS] += bs
    del net.trafo["_ppc_idx"]


def _add_ext_grid_sc_impedance_zero(net, ppc):
    mode = net["_options"]["mode"]


@@ -326,7 +332,7 @@


    z_grid = c / s_sc
    if mode == 'pf_3ph':
        z_grid = c / (s_sc/3)                       
        z_grid = c / (s_sc/3)
    x_grid = z_grid / np.sqrt(rx ** 2 + 1)
    r_grid = rx * x_grid
    eg["r"] = r_grid

@@ -348,7 +354,7 @@


def _add_line_sc_impedance_zero(net, ppc):
    branch_lookup = net["_pd2ppc_lookups"]["branch"]
    mode = net["_options"]["mode"]                            
    mode = net["_options"]["mode"]
    if not "line" in branch_lookup:
        return
    line = net["line"]

@@ -360,7 +366,7 @@

    tb = bus_lookup[line["to_bus"].values]
    baseR = np.square(ppc["bus"][fb, BASE_KV]) / net.sn_mva
    if mode == 'pf_3ph':
        baseR = np.square(ppc["bus"][fb, BASE_KV]) / (3*net.sn_mva)                     
        baseR = np.square(ppc["bus"][fb, BASE_KV]) / (3*net.sn_mva)
    f, t = branch_lookup["line"]
    # line zero sequence impedance
    ppc["branch"][f:t, F_BUS] = fb

@@ -373,3 +379,65 @@

373	379		ppc["branch"][f:t, BR_X] = line["x0_ohm_per_km"].values * length / baseR / parallel
374	380		ppc["branch"][f:t, BR_B] = (2 * net["f_hz"] * math.pi * line["c0_nf_per_km"].values * 1e-9 * baseR * length * parallel)
375	381		ppc["branch"][f:t, BR_STATUS] = line["in_service"].astype(int)
	382	+
	383	+
	384	+	def _add_trafo3w_sc_impedance_zero(net, ppc):
	385	+	branch_lookup = net["_pd2ppc_lookups"]["branch"]
	386	+	if not "trafo3w" in branch_lookup:
	387	+	return
	388	+	bus_lookup = net["_pd2ppc_lookups"]["bus"]
	389	+	branch = ppc["branch"]
	390	+	f, t = net["_pd2ppc_lookups"]["branch"]["trafo3w"]
	391	+	trafo_df = _trafo_df_from_trafo3w(net, seq=0)
	392	+	hv_bus = get_trafo_values(trafo_df, "hv_bus").astype(int)
	393	+	lv_bus = get_trafo_values(trafo_df, "lv_bus").astype(int)
	394	+	in_service = get_trafo_values(trafo_df, "in_service").astype(int)
	395	+	branch[f:t, F_BUS] = bus_lookup[hv_bus]
	396	+	branch[f:t, T_BUS] = bus_lookup[lv_bus]
	397	+
	398	+	r, x, _, ratio, shift = _calc_branch_values_from_trafo_df(net, ppc, trafo_df, seq=0)
	399	+
	400	+	n_t3 = net.trafo3w.shape[0]
	401	+	for t3_ix in np.arange(n_t3):
	402	+	t3 = net.trafo3w.iloc[t3_ix, :]
	403	+
	404	+	if t3.vector_group.lower() == "ynyd":
	405	+	# Correction for YnYD
	406	+	# z3->y3
	407	+	ys = 1 / ((x[t3_ix+n_t32] 1j + r[t3_ix+n_t32]) ratio[t3_ix+n_t32] * 2)
	408	+	aux_bus = bus_lookup[lv_bus[t3_ix]]
	409	+	ppc["bus"][aux_bus, BS] += ys.imag
	410	+	ppc["bus"][aux_bus, GS] += ys.real
	411	+
	412	+	# Set z2/z3 to almost 0 to avoid isolated bus
	413	+	x[[t3_ix+n_t3, t3_ix+n_t3*2]] = 1e20
	414	+	r[[t3_ix+n_t3, t3_ix+n_t3*2]] = 1e20
	415	+	elif t3.vector_group.lower() == "yynd":
	416	+	# z3->y3
	417	+	ys = 1 / ((x[t3_ix+n_t32] 1j + r[t3_ix+n_t32]) ratio[t3_ix+n_t32] * 2)
	418	+	aux_bus = bus_lookup[lv_bus[t3_ix]]
	419	+	ppc["bus"][aux_bus, BS] += ys.imag
	420	+	ppc["bus"][aux_bus, GS] += ys.real
	421	+
	422	+	# Set z1/z3 to almost 0 to avoid isolated bus
	423	+	x[[t3_ix, t3_ix+n_t3*2]] = 1e20
	424	+	r[[t3_ix, t3_ix+n_t3*2]] = 1e20
	425	+	elif t3.vector_group.lower() == "ynynd":
	426	+	# z3->y3
	427	+	ys = 1 / ((x[t3_ix+n_t32] 1j + r[t3_ix+n_t32]) ratio[t3_ix+n_t32] * 2)
	428	+	aux_bus = bus_lookup[lv_bus[t3_ix]]
	429	+	ppc["bus"][aux_bus, BS] += ys.imag
	430	+	ppc["bus"][aux_bus, GS] += ys.real
	431	+
	432	+	# Set z3 to almost 0 to avoid isolated bus
	433	+	x[t3_ix+n_t3*2] = 1e20
	434	+	r[t3_ix+n_t3*2] = 1e20
	435	+	else:
	436	+	raise UserWarning(f"{t3.vector_group} not supported yet for trafo3w!")
	437	+
	438	+	branch[f:t, BR_R] = r
	439	+	branch[f:t, BR_X] = x
	440	+	branch[f:t, BR_B] = 0
	441	+	branch[f:t, TAP] = ratio
	442	+	branch[f:t, SHIFT] = shift
	443	+	branch[f:t, BR_STATUS] = in_service

@@ -11,57 +11,72 @@

from pandapower.pypower.idx_bus import BASE_KV
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.idx_bus import C_MIN, C_MAX, KAPPA, R_EQUIV, IKSS1, IP, ITH,\
    X_EQUIV, IKSS2, IKCV, M, R_EQUIV_OHM, X_EQUIV_OHM, V_G, K_SG, SKSS
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, ppci, bus_idx):
    fault = net._options["fault"]
    case = net._options["case"]
    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"]
    c = ppci["bus"][bus_idx, C_MIN] if case == "min" else ppci["bus"][bus_idx, C_MAX]
    ppci["internal"]["baseI"] = ppci["bus"][:, BASE_KV] * np.sqrt(3) / ppci["baseMVA"]

    # Only for test, should correspondant to PF result
    baseZ = ppci["bus"][bus_idx, BASE_KV] ** 2 / ppci["baseMVA"]
    ppci["bus"][bus_idx, R_EQUIV_OHM] = baseZ * ppci["bus"][bus_idx, R_EQUIV]
    ppci["bus"][bus_idx, X_EQUIV_OHM] = baseZ * ppci["bus"][bus_idx, X_EQUIV]

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

    if fault == "3ph":
        ppci["bus"][bus_idx, SKSS] = np.sqrt(3) * ppci["bus"][bus_idx, IKSS1] * ppci["bus"][bus_idx, BASE_KV]
    elif fault == "2ph":
        ppci["bus"][bus_idx, SKSS] = ppci["bus"][bus_idx, IKSS1] * ppci["bus"][bus_idx, BASE_KV] / np.sqrt(3)

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)
    # Correct voltage of generator bus inside power station
    if np.any(~np.isnan(ppci["bus"][:, K_SG])):
        gen_bus_idx = bus_idx[~np.isnan(ppci["bus"][bus_idx, K_SG])]
        ppci["bus"][gen_bus_idx, IKSS1] *=\
            (ppci["bus"][gen_bus_idx, V_G] / ppci["bus"][gen_bus_idx, BASE_KV])
        ppci["bus"][gen_bus_idx, SKSS] *=\
            (ppci["bus"][gen_bus_idx, V_G] / ppci["bus"][gen_bus_idx, BASE_KV])

    _current_source_current(net, ppci)


def _calc_ikss_1ph(net, ppci, ppci_0, bus_idx):
    case = net._options["case"]
    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"]
    c = ppci["bus"][bus_idx, C_MIN] if case == "min" else ppci["bus"][bus_idx, C_MAX]
    ppci["internal"]["baseI"] = ppci["bus"][:, BASE_KV] * np.sqrt(3) / ppci["baseMVA"]
    ppci_0["internal"]["baseI"] = ppci_0["bus"][:, BASE_KV] * np.sqrt(3) / ppci_0["baseMVA"]

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

    # Only for test, should correspondant to PF result
    baseZ = ppci["bus"][bus_idx, BASE_KV] ** 2 / ppci["baseMVA"]
    ppci["bus"][bus_idx, R_EQUIV_OHM] = baseZ * ppci['bus'][bus_idx, R_EQUIV]
    ppci["bus"][bus_idx, X_EQUIV_OHM] = baseZ * ppci['bus'][bus_idx, X_EQUIV]
    ppci_0["bus"][bus_idx, R_EQUIV_OHM] = baseZ * ppci_0['bus'][bus_idx, R_EQUIV]
    ppci_0["bus"][bus_idx, X_EQUIV_OHM] = baseZ * ppci_0['bus'][bus_idx, X_EQUIV]

    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))
    ppci_0["bus"][bus_idx, IKSS1] = c / z_equiv / ppci_0["bus"][bus_idx, BASE_KV] * np.sqrt(3) * ppci_0["baseMVA"]
    ppci["bus"][bus_idx, IKSS1] = c / z_equiv / ppci["bus"][bus_idx, BASE_KV] * np.sqrt(3) * ppci["baseMVA"]

    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)
    _current_source_current(net, ppci)


def _current_source_current(net, ppc):
    ppc["bus"][:, IKCV] = 0
    ppc["bus"][:, IKSS2] = 0
def _current_source_current(net, ppci):
    ppci["bus"][:, IKCV] = 0
    ppci["bus"][:, IKSS2] = 0
    bus_lookup = net["_pd2ppc_lookups"]["bus"]
    if not False in net.sgen.current_source.values:
        sgen = net.sgen[net._is_elements["sgen"]]

@@ -71,227 +86,135 @@

        return
    if any(pd.isnull(sgen.sn_mva)):
        raise ValueError("sn_mva needs to be specified for all sgens in net.sgen.sn_mva")
    baseI = ppc["internal"]["baseI"]
    baseI = ppci["internal"]["baseI"]
    sgen_buses = sgen.bus.values
    sgen_buses_ppc = bus_lookup[sgen_buses]
    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
    ppci["bus"][buses, IKCV] = ikcv_pu
    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)
        Zbus = ppci["internal"]["Zbus"]
        ppci["bus"][:, IKSS2] = abs(1 / np.diag(Zbus) * np.dot(Zbus, ppci["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]
        ybus_fact = ppci["internal"]["ybus_fact"]
        diagZ = _calc_zbus_diag(net, ppci)
        ppci["bus"][:, IKSS2] = abs(ybus_fact(ppci["bus"][:, IKCV] * -1j) / diagZ / baseI)
    ppci["bus"][buses, IKCV] /= baseI[buses]


def _calc_ip(net, ppc):
    ip = np.sqrt(2) * (ppc["bus"][:, KAPPA] * ppc["bus"][:, IKSS1] + ppc["bus"][:, IKSS2])
    ppc["bus"][:, IP] = ip
def _calc_ip(net, ppci):
    ip = np.sqrt(2) * (ppci["bus"][:, KAPPA] * ppci["bus"][:, IKSS1] + ppci["bus"][:, IKSS2])
    ppci["bus"][:, IP] = ip


def _calc_ith(net, ppc):
def _calc_ith(net, ppci):
    tk_s = net["_options"]["tk_s"]
    kappa = ppc["bus"][:, KAPPA]
    kappa = ppci["bus"][:, KAPPA]
    f = 50
    n = 1
    m = (np.exp(4 * f * tk_s * np.log(kappa - 1)) - 1) / (2 * f * tk_s * np.log(kappa - 1))
    m[np.where(kappa > 1.99)] = 0
    ppc["bus"][:, M] = m
    ith = (ppc["bus"][:, IKSS1] + ppc["bus"][:, IKSS2]) * np.sqrt(m + n)
    ppc["bus"][:, ITH] = ith


def _calc_ib_generator(net, ppci):
    Zbus = ppci["internal"]["Zbus"]
    baseI = ppci["internal"]["baseI"]
    tk_s = net._options['tk_s']
    c = 1.1

    z_equiv = ppci["bus"][:, R_EQUIV] + ppci["bus"][:, X_EQUIV] * 1j
    I_ikss = c / z_equiv / ppci["bus"][:, BASE_KV] / np.sqrt(3) * ppci["baseMVA"]

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

    gen = net["gen"][net._is_elements["gen"]]
    gen_vn_kv = gen.vn_kv.values
    ppci["bus"][:, M] = m
    ith = (ppci["bus"][:, IKSS1] + ppci["bus"][:, IKSS2]) * np.sqrt(m + n)
    ppci["bus"][:, ITH] = ith

    gen_buses = ppci['gen'][:, GEN_BUS].astype(np.int64)
    gen_mbase = ppci['gen'][:, MBASE]
    gen_i_rg = gen_mbase / (np.sqrt(3) * gen_vn_kv)

    gen_buses_ppc, gen_sn_mva, I_rG = _sum_by_group(gen_buses, gen_mbase, gen_i_rg)
# TODO: Ib for generation close bus
# def _calc_ib_generator(net, ppci):
#     # Zbus = ppci["internal"]["Zbus"]
#     # baseI = ppci["internal"]["baseI"]
#     tk_s = net._options['tk_s']
#     c = 1.1

    # shunt admittance of generator buses and generator short circuit current
    # YS = ppci["bus"][gen_buses_ppc, GS] + ppci["bus"][gen_buses_ppc, BS] * 1j
    # I_kG = V_ikss.T[:, gen_buses_ppc] * YS / baseI[gen_buses_ppc]
#     z_equiv = ppci["bus"][:, R_EQUIV] + ppci["bus"][:, X_EQUIV] * 1j
#     I_ikss = c / z_equiv / ppci["bus"][:, BASE_KV] / np.sqrt(3) * ppci["baseMVA"]

    xdss_pu = gen.xdss_pu.values
    rdss_pu = gen.rdss_pu.values
    cosphi = gen.cos_phi.values
    X_dsss = xdss_pu * np.square(gen_vn_kv) / gen_mbase
    R_dsss = rdss_pu * np.square(gen_vn_kv) / gen_mbase
#     # calculate voltage source branch current
#     # I_ikss = ppci["bus"][:, IKSS1]
#     # V_ikss = (I_ikss * baseI) * Zbus

    K_G = ppci['bus'][gen_buses, BASE_KV] / gen_vn_kv * c / (1 + xdss_pu * np.sin(np.arccos(cosphi)))
    Z_G = (R_dsss + 1j * X_dsss)
#     gen = net["gen"][net._is_elements["gen"]]
#     gen_vn_kv = gen.vn_kv.values

    I_kG = c * ppci['bus'][gen_buses, BASE_KV] / np.sqrt(3) / (Z_G * K_G) * ppci["baseMVA"]
#     # Check difference ext_grid and gen
#     gen_buses = ppci['gen'][:, GEN_BUS].astype(np.int64)
#     gen_mbase = ppci['gen'][:, MBASE]
#     gen_i_rg = gen_mbase / (np.sqrt(3) * gen_vn_kv)

    dV_G = 1j * X_dsss * K_G * I_kG
    V_Is = c * ppci['bus'][gen_buses, BASE_KV] / np.sqrt(3)
#     gen_buses_ppc, gen_sn_mva, I_rG = _sum_by_group(gen_buses, gen_mbase, gen_i_rg)

    # I_kG_contribution = I_kG.sum(axis=1)
    # ratio_SG_ikss = I_kG_contribution / I_ikss
    # close_to_SG = ratio_SG_ikss > 5e-2
#     # shunt admittance of generator buses and generator short circuit current
#     # YS = ppci["bus"][gen_buses_ppc, GS] + ppci["bus"][gen_buses_ppc, BS] * 1j
#     # I_kG = V_ikss.T[:, gen_buses_ppc] * YS / baseI[gen_buses_ppc]

    close_to_SG = I_kG / I_rG > 2
#     xdss_pu = gen.xdss_pu.values
#     rdss_pu = gen.rdss_pu.values
#     cosphi = gen.cos_phi.values
#     X_dsss = xdss_pu * np.square(gen_vn_kv) / gen_mbase
#     R_dsss = rdss_pu * np.square(gen_vn_kv) / gen_mbase

    if tk_s == 2e-2:
        mu = 0.84 + 0.26 * np.exp(-0.26 * abs(I_kG) / I_rG)
    elif tk_s == 5e-2:
        mu = 0.71 + 0.51 * np.exp(-0.3 * abs(I_kG) / I_rG)
    elif tk_s == 10e-2:
        mu = 0.62 + 0.72 * np.exp(-0.32 * abs(I_kG) / I_rG)
    elif tk_s >= 25e-2:
        mu = 0.56 + 0.94 * np.exp(-0.38 * abs(I_kG) / I_rG)
    else:
        raise UserWarning('not implemented for other tk_s than 20ms, 50ms, 100ms and >=250ms')
#     K_G = ppci['bus'][gen_buses, BASE_KV] / gen_vn_kv * c / (1 + xdss_pu * np.sin(np.arccos(cosphi)))
#     Z_G = (R_dsss + 1j * X_dsss)

    mu = np.clip(mu, 0, 1)
#     I_kG = c * ppci['bus'][gen_buses, BASE_KV] / np.sqrt(3) / (Z_G * K_G) * ppci["baseMVA"]

    I_ikss_G = abs(I_ikss - np.sum((1 - mu) * I_kG, axis=1))
#     dV_G = 1j * X_dsss * K_G * I_kG
#     V_Is = c * ppci['bus'][gen_buses, BASE_KV] / np.sqrt(3)

    # I_ikss_G = I_ikss - np.sum(abs(V_ikss.T[:, gen_buses_ppc]) * (1-mu) * I_kG, axis=1)
#     # I_kG_contribution = I_kG.sum(axis=1)
#     # ratio_SG_ikss = I_kG_contribution / I_ikss
#     # close_to_SG = ratio_SG_ikss > 5e-2

    I_ikss_G = abs(I_ikss - np.sum(dV_G / V_Is * (1 - mu) * I_kG, axis=1))
#     close_to_SG = I_kG / I_rG > 2

    return I_ikss_G
#     if tk_s == 2e-2:
#         mu = 0.84 + 0.26 * np.exp(-0.26 * abs(I_kG) / I_rG)
#     elif tk_s == 5e-2:
#         mu = 0.71 + 0.51 * np.exp(-0.3 * abs(I_kG) / I_rG)
#     elif tk_s == 10e-2:
#         mu = 0.62 + 0.72 * np.exp(-0.32 * abs(I_kG) / I_rG)
#     elif tk_s >= 25e-2:
#         mu = 0.56 + 0.94 * np.exp(-0.38 * abs(I_kG) / I_rG)
#     else:
#         raise UserWarning('not implemented for other tk_s than 20ms, 50ms, 100ms and >=250ms')

#     mu = np.clip(mu, 0, 1)

def _calc_single_bus_sc(net, ppc, bus):
    # case = net._options["case"]
    bus_idx = net._pd2ppc_lookups["bus"][bus]
    n = ppc["bus"].shape[0]
    Zbus = ppc["internal"]["Zbus"]
    #    Yf = ppc["internal"]["Yf"]
    #    Yt = ppc["internal"]["Yf"]
    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]
#     I_ikss_G = abs(I_ikss - np.sum((1 - mu) * I_kG, axis=1))

    # calculate voltage source branch current
    V_ikss = (ppc["bus"][:, IKSS1] * baseI) * Zbus
    V = V_ikss[:, bus_idx]
    #    ikss_all_f = np.conj(Yf.dot(V_ikss))
    #    ikss_all_t = np.conj(Yt.dot(V_ikss))
    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_source = np.dot((current * baseI), Zbus).T
        V = V + V_source[:, bus_idx]
    # add current source branch current if there is one
    #    ppc["branch"][:, IKSS_F] = abs(ikss_all_f[:, bus_idx] / baseI[fb])
    #    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)

    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)
#     # I_ikss_G = I_ikss - np.sum(abs(V_ikss.T[:, gen_buses_ppc]) * (1-mu) * I_kG, axis=1)

#     I_ikss_G = abs(I_ikss - np.sum(dV_G / V_Is * (1 - mu) * I_kG, axis=1))

def calc_branch_results(net, ppci, V):
    Ybus = ppci["internal"]["Ybus"]
    Yf = ppci["internal"]["Yf"]
    Yt = ppci["internal"]["Yt"]
    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
#     return I_ikss_G


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]]
def _calc_branch_currents(net, ppci, bus_idx):
    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

    Yf = ppci["internal"]["Yf"]
    Yt = ppci["internal"]["Yt"]
    baseI = ppci["internal"]["baseI"]
    n_bus = ppci["bus"].shape[0]
    fb = np.real(ppci["branch"][:, 0]).astype(int)
    tb = np.real(ppci["branch"][:, 1]).astype(int)

    # calculate voltage source branch current
    if net["_options"]["inverse_y"]:
        Zbus = ppc["internal"]["Zbus"]
        V_ikss = (ppc["bus"][:, IKSS1] * baseI) * Zbus
        Zbus = ppci["internal"]["Zbus"]
        V_ikss = (ppci["bus"][:, IKSS1] * baseI) * Zbus
        V_ikss = V_ikss[:, bus_idx]
    else:
        ybus_fact = ppc["internal"]["ybus_fact"]
        ybus_fact = ppci["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]
            ikss[b] = ppci["bus"][b, IKSS1] * baseI[b]
            V_ikss[:, ix] = ybus_fact(ikss)

    ikss1_all_f = np.conj(Yf.dot(V_ikss))

@@ -300,24 +223,24 @@

    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
    current_sources = any(ppci["bus"][:, IKCV]) > 0
    if current_sources:
        current = np.tile(-ppc["bus"][:, IKCV], (n_sc_bus, 1))
        current = np.tile(-ppci["bus"][:, IKCV], (n_sc_bus, 1))
        for ix, b in enumerate(bus_idx):
            current[ix, b] += ppc["bus"][b, IKSS2]
            current[ix, b] += ppci["bus"][b, IKSS2]

        # calculate voltage source branch current
        if net["_options"]["inverse_y"]:
            Zbus = ppc["internal"]["Zbus"]
            Zbus = ppci["internal"]["Zbus"]
            V = np.dot((current * baseI), Zbus).T
        else:
            ybus_fact = ppc["internal"]["ybus_fact"]
            ybus_fact = ppci["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)
        fb = np.real(ppci["branch"][:, 0]).astype(int)
        tb = np.real(ppci["branch"][:, 1]).astype(int)
        ikss2_all_f = np.conj(Yf.dot(V))
        ikss2_all_t = np.conj(Yt.dot(V))


@@ -328,16 +251,16 @@

        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]
        ppci["internal"]["branch_ikss_f"] = ikss_all_f / baseI[fb, None]
        ppci["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])
        ppci["branch"][:, IKSS_F] = np.nan_to_num(minmax(ikss_all_f, axis=1) / baseI[fb])
        ppci["branch"][:, IKSS_T] = np.nan_to_num(minmax(ikss_all_t, axis=1) / baseI[tb])

    if net._options["ip"]:
        kappa = ppc["bus"][:, KAPPA]
        kappa = ppci["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)

@@ -346,23 +269,27 @@

            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]
            ppci["internal"]["branch_ip_f"] = abs(ip_all_f) / baseI[fb, None]
            ppci["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])
            ppci["branch"][:, IP_F] = np.nan_to_num(minmax(abs(ip_all_f), axis=1) / baseI[fb])
            ppci["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]
        m = ppci["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]
            ppci["internal"]["branch_ith_f"] = ith_all_f / baseI[fb, None]
            ppci["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])
            ppci["branch"][:, ITH_F] = np.nan_to_num(minmax(ith_all_f, axis=1) / baseI[fb])
            ppci["branch"][:, ITH_T] = np.nan_to_num(minmax(ith_all_t, axis=1) / baseI[fb])

    # Update bus index for branch results
    if net._options["return_all_currents"]:
        ppci["internal"]["br_res_ks_ppci_bus"] = bus_idx

@@ -3,8 +3,6 @@

# Copyright (c) 2016-2021 by University of Kassel and Fraunhofer Institute for Energy Economics
# and Energy System Technology (IEE), Kassel. All rights reserved.



from pandapower.pypower.idx_brch import branch_cols as start

IKSS_F    = start + 0

@@ -17,5 +15,7 @@

IK_T      = start + 7
IB_F      = start + 8
IB_T      = start + 9
K_T       = start + 10
K_ST       = start + 11

branch_cols_sc = 10

@@ -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.

@@ -506,12 +506,7 @@

        ppc["bus"][b, GS] = vp
        ppc["bus"][b, BS] = -vq


def _add_gen_impedances_ppc(net, ppc):
    _add_ext_grid_sc_impedance(net, ppc)
    _add_gen_sc_impedance(net, ppc)


# Short circuit relevant routines
def _add_ext_grid_sc_impedance(net, ppc):
    from pandapower.shortcircuit.idx_bus import C_MAX, C_MIN
    mode = net._options["mode"]

@@ -532,11 +527,11 @@

        c = 1.1
    if not "s_sc_%s_mva" % case in eg:
        raise ValueError("short circuit apparent power s_sc_%s_mva needs to be specified for "% case +
                         "external grid \n Try: net.ext_grid['s_sc_max_mva'] = 1000" )
                          "external grid \n Try: net.ext_grid['s_sc_max_mva'] = 1000" )
    s_sc = eg["s_sc_%s_mva" % case].values/ppc['baseMVA']
    if not "rx_%s" % case in eg:
        raise ValueError("short circuit R/X rate rx_%s needs to be specified for external grid \n Try: net.ext_grid['rx_max'] = 0.1" %
                         case)
                          case)
    rx = eg["rx_%s" % case].values

    z_grid = c / s_sc

@@ -553,44 +548,6 @@

    ppc["bus"][buses, BS] = bs * ppc['baseMVA']
    return gs * ppc['baseMVA'], bs * ppc['baseMVA']


def _add_gen_sc_impedance(net, ppc):
    from pandapower.shortcircuit.idx_bus import C_MAX
    gen = net["gen"][net._is_elements["gen"]]
    if len(gen) == 0:
        return
    gen_buses = gen.bus.values
    bus_lookup = net["_pd2ppc_lookups"]["bus"]
    gen_buses_ppc = bus_lookup[gen_buses]
    vn_gen = gen.vn_kv.values
    sn_gen = gen.sn_mva.values

    rdss_pu = gen.rdss_pu.values
    xdss_pu = gen.xdss_pu.values
    gens_without_r = np.isnan(rdss_pu)
    if gens_without_r.any():
        #  use the estimations from the IEC standard for generators without defined rdss_pu
        lv_gens = (vn_gen <= 1.) & gens_without_r
        hv_gens = (vn_gen > 1.) & gens_without_r
        large_hv_gens = (sn_gen >= 100) & hv_gens
        small_hv_gens = (sn_gen < 100) & hv_gens
        rdss_pu[lv_gens] = 0.15 * xdss_pu[lv_gens]
        rdss_pu[large_hv_gens] = 0.05 * xdss_pu[large_hv_gens]
        rdss_pu[small_hv_gens] = 0.07 * xdss_pu[small_hv_gens]

    vn_net = ppc["bus"][gen_buses_ppc, BASE_KV]
    cmax = ppc["bus"][gen_buses_ppc, C_MAX]
    phi_gen = np.arccos(gen.cos_phi.values)
    kg = vn_gen / vn_net * cmax / (1 + xdss_pu * np.sin(phi_gen))

    z_gen = (rdss_pu + xdss_pu * 1j) * kg / sn_gen
    y_gen = 1 / z_gen

    buses, gs, bs = _sum_by_group(gen_buses_ppc, y_gen.real, y_gen.imag)
    ppc["bus"][buses, GS] = gs
    ppc["bus"][buses, BS] = bs


def _add_motor_impedances_ppc(net, ppc):
    if net._options["case"] == "min":
        return

@@ -20,5 +20,14 @@

IB          = start + 10
ITH         = start + 11
IK          = start + 12
R_EQUIV_OHM = start + 13
X_EQUIV_OHM = start + 14
K_G         = start + 15
K_SG        = start + 16
V_G         = start + 17
PS_TRAFO_IX = start + 18
GS_P        = start + 19
BS_P        = start + 20
SKSS        = start + 21

bus_cols_sc = 13

@@ -0,0 +1,187 @@

1	+	# -- coding: utf-8 --
2	+
3	+	# Copyright (c) 2016-2021 by University of Kassel and Fraunhofer Institute for Energy Economics
4	+	# and Energy System Technology (IEE), Kassel. All rights reserved.
5	+
6	+	from copy import deepcopy
7	+	import numpy as np
8	+
9	+	from pandapower.pd2ppc import _pd2ppc, _ppc2ppci
10	+	from pandapower.auxiliary import _add_auxiliary_elements, _sum_by_group
11	+	from pandapower.build_branch import get_trafo_values, _transformer_correction_factor
12	+
13	+	from pandapower.pypower.idx_bus import BASE_KV, GS, BS
14	+	from pandapower.pypower.idx_brch import BR_X, BR_R
15	+	from pandapower.shortcircuit.idx_bus import C_MAX, K_G, K_SG, V_G,\
16	+	PS_TRAFO_IX, GS_P, BS_P
17	+	from pandapower.shortcircuit.idx_brch import K_T, K_ST
18	+
19	+
20	+	def _get_is_ppci_bus(net, bus):
21	+	is_bus = bus[np.in1d(bus, net._is_elements_final["bus_is_idx"])]
22	+	ppci_bus = np.unique(net._pd2ppc_lookups["bus"][is_bus])
23	+	return ppci_bus
24	+
25	+
26	+	def _init_ppc(net):
27	+	_check_sc_data_integrity(net)
28	+	_add_auxiliary_elements(net)
29	+	ppc, _ = _pd2ppc(net)
30	+
31	+	# Init the required columns to nan
32	+	ppc["bus"][:, [K_G, K_SG, V_G, PS_TRAFO_IX, GS_P, BS_P,]] = np.nan
33	+	ppc["branch"][:, [K_T, K_ST]] = np.nan
34	+
35	+	# Add parameter K into ppc
36	+	_add_kt(net, ppc)
37	+	_add_gen_sc_z_kg_ks(net, ppc)
38	+
39	+	ppci = _ppc2ppci(ppc, net)
40	+	return ppc, ppci
41	+
42	+
43	+	def _check_sc_data_integrity(net):
44	+	if not net.gen.empty:
45	+	for col in ("power_station_trafo", "pg_percent"):
46	+	if col not in net.gen.columns:
47	+	net.gen[col] = np.nan
48	+
49	+	if not net.trafo.empty:
50	+	for col in ("pt_percent",):
51	+	if col not in net.trafo.columns:
52	+	net.trafo[col] = np.nan
53	+
54	+
55	+	def _add_kt(net, ppc):
56	+	bus_lookup = net["_pd2ppc_lookups"]["bus"]
57	+	branch = ppc["branch"]
58	+	# "trafo/trafo3w" are already corrected in pd2ppc, write parameter kt for trafo
59	+	if not net.trafo.empty:
60	+	f, t = net["_pd2ppc_lookups"]["branch"]["trafo"]
61	+	trafo_df = net["trafo"]
62	+	cmax = ppc["bus"][bus_lookup[get_trafo_values(trafo_df, "lv_bus")], C_MAX]
63	+	kt = _transformer_correction_factor(trafo_df.vk_percent, trafo_df.vkr_percent,
64	+	trafo_df.sn_mva, cmax)
65	+	branch[f:t, K_T] = kt
66	+
67	+
68	+	def _add_gen_sc_z_kg_ks(net, ppc):
69	+	gen = net["gen"][net._is_elements["gen"]]
70	+	if len(gen) == 0:
71	+	return
72	+	gen_buses = gen.bus.values
73	+	bus_lookup = net["_pd2ppc_lookups"]["bus"]
74	+	gen_buses_ppc = bus_lookup[gen_buses]
75	+	vn_gen = gen.vn_kv.values
76	+	sn_gen = gen.sn_mva.values
77	+
78	+	rdss_ohm = gen.rdss_ohm.values
79	+	xdss_pu = gen.xdss_pu.values
80	+	pg_percent = gen.pg_percent.values
81	+	vn_net = ppc["bus"][gen_buses_ppc, BASE_KV]
82	+	sin_phi_gen = np.sin(np.arccos(gen.cos_phi.values))
83	+
84	+	gen_baseZ = (vn_net ** 2 / ppc["baseMVA"])
85	+	r_gen, x_gen = rdss_ohm, xdss_pu * vn_gen ** 2 / sn_gen
86	+	z_gen = (r_gen + x_gen * 1j)
87	+	z_gen_pu = z_gen / gen_baseZ
88	+	y_gen_pu = 1 / z_gen_pu
89	+
90	+	buses, gs, bs = _sum_by_group(gen_buses_ppc, y_gen_pu.real, y_gen_pu.imag)
91	+	ppc["bus"][buses, GS] = gs
92	+	ppc["bus"][buses, BS] = bs
93	+
94	+	# Calculate K_G
95	+	cmax = ppc["bus"][gen_buses_ppc, C_MAX]
96	+	kg = (1/(1+pg_percent/100)) * cmax / (1 + (xdss_pu * sin_phi_gen))
97	+	ppc["bus"][gen_buses_ppc, K_G] = kg
98	+	ppc["bus"][gen_buses_ppc, V_G] = vn_gen
99	+
100	+	# Calculate G,B (r,x) on generator for peak current calculation
101	+	r_gen_p = np.full_like(r_gen, fill_value=np.nan)
102	+	lv_gens, hv_gens = (vn_gen <= 1.), (vn_gen > 1.)
103	+	small_hv_gens = (sn_gen < 100) & hv_gens
104	+	large_hv_gens = (sn_gen >= 100) & hv_gens
105	+	if np.any(lv_gens):
106	+	r_gen_p[lv_gens] = 0.15 * x_gen[lv_gens]
107	+	if np.any(small_hv_gens):
108	+	r_gen_p[small_hv_gens] = 0.07 * x_gen[small_hv_gens]
109	+	if np.any(large_hv_gens):
110	+	r_gen_p[large_hv_gens] = 0.05 * x_gen[large_hv_gens]
111	+	z_gen_p = (r_gen_p + x_gen * 1j)
112	+	z_gen_p_pu = z_gen_p / gen_baseZ
113	+	y_gen_p_pu = 1 / z_gen_p_pu
114	+
115	+	_, ppc["bus"][buses, GS_P], ppc["bus"][buses, BS_P] =\
116	+	_sum_by_group(gen_buses_ppc, y_gen_p_pu.real, y_gen_p_pu.imag)
117	+
118	+	# Calculate K_S on power station configuration
119	+	if np.any(~np.isnan(net.gen.power_station_trafo.values)):
120	+	ps_gen_ix = net.gen.loc[~np.isnan(net.gen.power_station_trafo.values),:].index.values
121	+	ps_trafo_ix = net.gen.loc[ps_gen_ix, "power_station_trafo"].values.astype(int)
122	+	ps_trafo_with_tap_mask = (~np.isnan(net.trafo.loc[ps_trafo_ix, "pt_percent"]))
123	+	ps_gen_buses_ppc = bus_lookup[net.gen.loc[ps_gen_ix, "bus"]]
124	+	f, t = net["_pd2ppc_lookups"]["branch"]["trafo"]
125	+
126	+	ps_cmax = ppc["bus"][ps_gen_buses_ppc, C_MAX]
127	+	v_trafo_hv, v_trafo_lv =\
128	+	net.trafo.loc[ps_trafo_ix, "vn_hv_kv"].values, net.trafo.loc[ps_trafo_ix, "vn_lv_kv"].values
129	+	v_q = net.bus.loc[net.trafo.loc[ps_trafo_ix, "hv_bus"], "vn_kv"].values
130	+	v_g = vn_gen[ps_gen_ix]
131	+	x_g = xdss_pu[ps_gen_ix]
132	+	x_t = net.trafo.loc[ps_trafo_ix, "vk_percent"].values / 100
133	+
134	+	if np.any(ps_trafo_with_tap_mask):
135	+	ks = (v_q2/v_g2) * (v_trafo_lv2/v_trafo_hv2) *\
136	+	ps_cmax / (1 + np.abs(x_g - x_t) * sin_phi_gen[ps_gen_ix])
137	+	ppc["bus"][ps_gen_buses_ppc[ps_trafo_with_tap_mask], K_SG] = ks[ps_trafo_with_tap_mask]
138	+	ppc["branch"][np.arange(f, t)[ps_trafo_ix][ps_trafo_with_tap_mask], K_ST] = ks[ps_trafo_with_tap_mask]
139	+
140	+	if np.any(~ps_trafo_with_tap_mask):
141	+	kso = (v_q/v_g/(1 + pg_percent[ps_gen_ix] / 100)) * (v_trafo_lv/v_trafo_hv) *\
142	+	ps_cmax / (1 + x_g * sin_phi_gen[ps_gen_ix])
143	+	ppc["bus"][ps_gen_buses_ppc[~ps_trafo_with_tap_mask], K_SG] = kso[~ps_trafo_with_tap_mask]
144	+	ppc["branch"][np.arange(f, t)[ps_trafo_ix][~ps_trafo_with_tap_mask], K_ST] = kso[~ps_trafo_with_tap_mask]
145	+
146	+	ppc["bus"][ps_gen_buses_ppc, PS_TRAFO_IX] = np.arange(f, t)[ps_trafo_ix]
147	+
148	+
149	+	def _create_k_updated_ppci(net, ppci_orig, ppci_bus):
150	+	ppci = deepcopy(ppci_orig)
151	+
152	+	non_ps_gen_bus = ppci_bus[np.isnan(ppci["bus"][ppci_bus, K_SG])]
153	+	ps_gen_bus = ppci_bus[~np.isnan(ppci["bus"][ppci_bus, K_SG])]
154	+
155	+	ps_gen_bus_mask = ~np.isnan(ppci["bus"][:, K_SG])
156	+	ps_trafo_mask = ~np.isnan(ppci["branch"][:, K_ST])
157	+	if np.any(ps_gen_bus_mask):
158	+	ppci["bus"][np.ix_(ps_gen_bus_mask, [GS, BS, GS_P, BS_P])] /=\
159	+	ppci["bus"][np.ix_(ps_gen_bus_mask, [K_SG])]
160	+	ppci["branch"][np.ix_(ps_trafo_mask, [BR_X, BR_R])] *=\
161	+	ppci["branch"][np.ix_(ps_trafo_mask, [K_ST])] / ppci["branch"][np.ix_(ps_trafo_mask, [K_T])]
162	+
163	+	gen_bus_mask = np.isnan(ppci["bus"][:, K_SG]) & (~np.isnan(ppci["bus"][:, K_G]))
164	+	if np.any(gen_bus_mask):
165	+	ppci["bus"][np.ix_(gen_bus_mask, [GS, BS, GS_P, BS_P])] /=\
166	+	ppci["bus"][np.ix_(gen_bus_mask, [K_G])]
167	+
168	+	bus_ppci = {}
169	+	if ps_gen_bus.size > 0:
170	+	for bus in ps_gen_bus:
171	+	ppci_gen = deepcopy(ppci)
172	+	if not np.isfinite(ppci_gen["bus"][bus, K_SG]):
173	+	raise UserWarning("Parameter error of K SG")
174	+	# Correct ps gen bus
175	+	ppci_gen["bus"][bus, [GS, BS, GS_P, BS_P]] /=\
176	+	(ppci_gen["bus"][bus, K_G] / ppci_gen["bus"][bus, K_SG])
177	+
178	+	# Correct ps transfomer
179	+	trafo_ix = ppci_gen["bus"][bus, PS_TRAFO_IX].astype(int)
180	+	# TODO: Check this!
181	+	# ppci_gen["branch"][trafo_ix, [BR_X, BR_R]] *= \
182	+	# (ppci_gen["branch"][trafo_ix, K_T] / ppci_gen["branch"][trafo_ix, K_ST])
183	+	ppci_gen["branch"][trafo_ix, [BR_X, BR_R]] /= \
184	+	ppci_gen["branch"][trafo_ix, K_ST]
185	+
186	+	bus_ppci.update({bus: ppci_gen})
187	+	return non_ps_gen_bus, ppci, bus_ppci

@@ -18,19 +18,21 @@

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, _calc_single_bus_sc_no_y_inv

from pandapower.shortcircuit.currents import _calc_ikss,\
    _calc_ikss_1ph, _calc_ip, _calc_ith, _calc_branch_currents
from pandapower.shortcircuit.impedance import _calc_zbus, _calc_ybus, _calc_rx
from pandapower.shortcircuit.ppc_conversion import _init_ppc, _create_k_updated_ppci, _get_is_ppci_bus
from pandapower.shortcircuit.kappa import _add_kappa_to_ppc
from pandapower.shortcircuit.results import _extract_results, _extract_single_results
from pandapower.shortcircuit.results import _extract_results, _copy_result_to_ppci_orig
from pandapower.results import init_results


def calc_sc(net, fault="3ph", case='max', lv_tol_percent=10, topology="auto", ip=False,
def calc_sc(net, bus=None,
            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,
            bus=None, inverse_y=True):
            inverse_y=True):
    """
    Calculates minimal or maximal symmetrical short-circuit currents.
    The calculation is based on the method of the equivalent voltage source

@@ -46,6 +48,8 @@

    INPUT:
        **net** (pandapowerNet) pandapower Network

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

        ***fault** (str, 3ph) type of fault

            - "3ph" for three-phase

@@ -68,7 +72,7 @@


        **ip** (bool, False) if True, calculate aperiodic short-circuit current

        **Ith** (bool, False) if True, calculate equivalent thermical short-circuit current Ith
        **ith** (bool, False) if True, calculate equivalent thermical short-circuit current Ith

        **topology** (str, "auto") define option for meshing (only relevant for ip and ith)


@@ -89,8 +93,6 @@

        **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



@@ -106,12 +108,13 @@

            "Only 3ph, 2ph and 1ph short-circuit currents implemented")

    if len(net.gen) and (ip or ith):
        logger.warning("aperiodic and thermal short-circuit currents are only implemented for "
        logger.warning("aperiodic, thermal short-circuit currents are only implemented for "
                       "faults far from generators!")

    if case not in ['max', 'min']:
        raise ValueError('case can only be "min" or "max" for minimal or maximal short "\
                                "circuit current')

    if topology not in ["meshed", "radial", "auto"]:
        raise ValueError(
            'specify network structure as "meshed", "radial" or "auto"')

@@ -120,24 +123,27 @@

        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):
    # Convert bus to numpy array
    if bus is None:
        bus = net.bus.index.values
    elif isinstance(bus, int):
        bus = np.array([bus])
    elif isinstance(bus, list):
    elif not isinstance(bus, np.ndarray):
        bus = np.array(bus)

    kappa = ith or ip
    net["_options"] = {}
    _add_ppc_options(net, calculate_voltage_angles=False, trafo_model="pi",
                     check_connectivity=check_connectivity, mode="sc", switch_rx_ratio=2,
                     init_vm_pu="flat", init_va_degree="flat", enforce_q_lims=False,
                     recycle=None)
                      check_connectivity=check_connectivity, mode="sc", switch_rx_ratio=2,
                      init_vm_pu="flat", init_va_degree="flat", enforce_q_lims=False,
                      recycle=None)
    _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,
                    inverse_y=inverse_y)
    init_results(net, "sc")

    if fault in ("2ph", "3ph"):
        _calc_sc(net, bus)
    elif fault == "1ph":

@@ -146,128 +152,51 @@

        raise ValueError("Invalid fault %s" % fault)


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
    according to DIN/IEC EN 60909.
    The initial short-circuit alternating current *ikss* is the basis of the short-circuit
    calculation and is therefore always calculated.
    Other short-circuit currents can be calculated from *ikss* with the conversion factors defined
    in DIN/IEC EN 60909.

    The output is stored in the net.res_bus_sc table as a short_circuit current
    for each bus.

    INPUT:
        **net** (pandapowerNet) pandapower Network

        ***fault** (str, 3ph) type of fault

            - "3ph" for three-phase

            - "2ph" for two-phase short-circuits

        **case** (str, "max")

            - "max" for maximal current calculation

            - "min" for minimal current calculation

        **lv_tol_percent** (int, 10) voltage tolerance in low voltage grids

            - 6 for 6% voltage tolerance

            - 10 for 10% voltage olerance

        **r_fault_ohm** (float, 0) fault resistance in Ohm

        **x_fault_ohm** (float, 0) fault reactance in Ohm

    OUTPUT:

    EXAMPLE:
        calc_sc(net)
def _calc_current(net, ppci_orig, bus):
    # Select required ppci bus
    ppci_bus = _get_is_ppci_bus(net, bus)

        print(net.res_bus_sc)
    """
    if fault not in ["3ph", "2ph"]:
        raise NotImplementedError("Only 3ph and 2ph short-circuit currents implemented")

    if case not in ['max', 'min']:
        raise ValueError('case can only be "min" or "max" for minimal or maximal short "\
                                "circuit current')
    net["_options"] = {}
    _add_ppc_options(net, calculate_voltage_angles=False, trafo_model="pi",
                     check_connectivity=check_connectivity, mode="sc", switch_rx_ratio=2,
                     init_vm_pu="flat", init_va_degree="flat", enforce_q_lims=False,
                     recycle=None)
    _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,
                    inverse_y=inverse_y)
    init_results(net, "sc")
    if fault in ("2ph", "3ph"):
        _calc_sc_single(net, bus)
    elif fault == "1ph":
        raise NotImplementedError("1ph short-circuits are not yet implemented")
    else:
        raise ValueError("Invalid fault %s" % fault)
    # update ppci
    non_ps_gen_ppci_bus, non_ps_gen_ppci, ps_gen_bus_ppci_dict =\
        _create_k_updated_ppci(net, ppci_orig, ppci_bus=ppci_bus)

    # For each ps_gen_bus one unique ppci is required
    ps_gen_ppci_bus = list(ps_gen_bus_ppci_dict.keys())

def _calc_sc_single(net, bus):
    _add_auxiliary_elements(net)
    ppc, ppci = _pd2ppc(net)
    _calc_ybus(ppci)
    for calc_bus in ps_gen_ppci_bus+[non_ps_gen_ppci_bus]:
        if isinstance(calc_bus, np.ndarray):
            # Use ppci for general bus
            this_ppci, this_ppci_bus = non_ps_gen_ppci, calc_bus
        else:
            # Use specific ps_gen_bus ppci
            this_ppci, this_ppci_bus = ps_gen_bus_ppci_dict[calc_bus], np.array([calc_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_ybus(this_ppci)
        if net["_options"]["inverse_y"]:
            _calc_zbus(net, this_ppci)
        else:
            # Factorization Ybus once
            this_ppci["internal"]["ybus_fact"] = factorized(this_ppci["internal"]["Ybus"])

        _calc_rx(net, ppci, bus)
        _calc_ikss(net, ppci, bus)
        _calc_single_bus_sc_no_y_inv(net, ppci, bus)
        _calc_rx(net, this_ppci, this_ppci_bus)
        _calc_ikss(net, this_ppci, this_ppci_bus)
        _add_kappa_to_ppc(net, this_ppci)
        if net["_options"]["ip"]:
            _calc_ip(net, this_ppci)
        if net["_options"]["ith"]:
            _calc_ith(net, this_ppci)

        # Delete factorization object
        ppci["internal"].pop("ybus_fact")
        if net._options["branch_results"]:
            _calc_branch_currents(net, this_ppci, this_ppci_bus)

    ppc = _copy_results_ppci_to_ppc(ppci, ppc, "sc")
    _extract_single_results(net, ppc)
    _clean_up(net)
        _copy_result_to_ppci_orig(ppci_orig, this_ppci, this_ppci_bus,
                                  calc_options=net._options)


def _calc_sc(net, bus):
    _add_auxiliary_elements(net)
    ppc, ppci = _pd2ppc(net)
    _calc_ybus(ppci)
    ppc, ppci = _init_ppc(net)

    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, bus)
    _calc_current(net, ppci, bus)

    ppc = _copy_results_ppci_to_ppc(ppci, ppc, "sc")
    _extract_results(net, ppc, ppc_0=None, bus=bus)

@@ -299,14 +228,16 @@

        ppci["internal"]["ybus_fact"] = factorized(ppci["internal"]["Ybus"])
        ppci_0["internal"]["ybus_fact"] = factorized(ppci_0["internal"]["Ybus"])

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

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

    if net._options["branch_results"]:
        _calc_branch_currents(net, ppci, bus=bus)
        _calc_branch_currents(net, ppci, ppci_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, bus=bus)

@@ -1, +1, @@

from pandapower.shortcircuit.calc_sc import calc_sc, calc_single_sc

@@ -3,26 +3,62 @@

# Copyright (c) 2016-2021 by University of Kassel and Fraunhofer Institute for Energy Economics
# and Energy System Technology (IEE), Kassel. All rights reserved.


import numpy as np
import pandas as pd

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, 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
from pandapower.shortcircuit.idx_bus import IKSS1, IP, ITH, IKSS2, R_EQUIV_OHM, X_EQUIV_OHM, SKSS
from pandapower.pypower.idx_bus import BUS_TYPE

BRANCH_RESULTS_KEYS = ("branch_ikss_f", "branch_ikss_t",
                       "branch_ip_f", "branch_ip_t",
                       "branch_ith_f", "branch_ith_t")


def _copy_result_to_ppci_orig(ppci_orig, ppci, ppci_bus, calc_options):
    if ppci_orig is ppci:
        return

    ppci_orig["bus"][ppci_bus, :] = ppci["bus"][ppci_bus, :]
    if calc_options["branch_results"]:
        if calc_options["return_all_currents"]:
            ppci_orig["internal"]["br_res_ks_ppci_bus"] =\
                ppci_bus if "br_res_ks_ppci_bus" not in ppci_orig["internal"]\
                else np.r_[ppci_orig["internal"]["br_res_ks_ppci_bus"], ppci_bus]

            for res_key in BRANCH_RESULTS_KEYS:
                # Skip not required data points
                if res_key not in ppci["internal"]:
                    continue

                if res_key not in ppci_orig["internal"]:
                    ppci_orig["internal"][res_key] = ppci["internal"][res_key]
                else:
                    ppci_orig["internal"][res_key] = np.c_[ppci_orig["internal"][res_key],
                                                           ppci["internal"][res_key]]
        else:
            case = calc_options["case"]
            branch_results_cols = [IKSS_F, IKSS_T, IP_F, IP_T, ITH_F, ITH_T]
            if case == "max":
                ppci_orig["branch"][:, branch_results_cols] =\
                    np.maximum(np.nan_to_num(ppci["branch"][:, branch_results_cols]),
                               np.nan_to_num(ppci_orig["branch"][:, branch_results_cols]))
            else:
                ppci_orig["branch"][:, branch_results_cols] =\
                    np.minimum(np.nan_to_num(ppci["branch"][:, branch_results_cols], nan=1e10),
                               np.nan_to_num(ppci_orig["branch"][:, branch_results_cols], nan=1e10))


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 = bus_lookup[bus]
    ppc_index[ppc["bus"][bus_lookup[ppc_index], BUS_TYPE] == 4] = -1
    return bus, ppc_index


def _extract_results(net, ppc, ppc_0, bus):
    _get_bus_results(net, ppc, ppc_0, bus)
    if net._options["branch_results"]:

@@ -35,45 +71,27 @@

            _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,
                                                columns=net["_empty_res_%s"%element].columns,
                                                dtype='float')
    _get_single_bus_results(net, ppc)
    net["_options"]["ac"] = True
    net["_options"]["trafo_loading"] = "current"
    bus_lookup_aranged = _get_aranged_lookup(net)
    bus_pq = np.zeros(shape=(len(net["bus"].index), 2), dtype=np.float)
    _get_branch_results(net, ppc, bus_lookup_aranged, bus_pq, suffix="_sc")


def _get_single_bus_results(net, ppc):
    _set_buses_out_of_service(ppc)
    bus_idx = _get_bus_idx(net)
    case = net._options["case"]
    c = ppc["bus"][bus_idx, C_MIN] if case == "min" else ppc["bus"][bus_idx, C_MAX]
    net["res_bus"]["vm_pu"] = np.nan
    net["res_bus_sc"]["vm_pu"] = c - ppc["bus"][bus_idx, VM]
    net["res_bus_sc"]["va_degree"] = ppc["bus"][bus_idx, VA]


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]

    if net["_options"]["fault"] == "1ph":
        net.res_bus_sc["ikss_ka"] = ppc_0["bus"][ppc_index, IKSS1] + ppc["bus"][ppc_index, IKSS2]
        net.res_bus_sc["rk0_ohm"] = ppc_0["bus"][ppc_index, R_EQUIV_OHM]
        net.res_bus_sc["xk0_ohm"] = ppc_0["bus"][ppc_index, X_EQUIV_OHM]
    else:
        net.res_bus_sc["ikss_ka"] = ppc["bus"][ppc_index, IKSS1] + ppc["bus"][ppc_index, IKSS2]
        net.res_bus_sc["skss_mw"] = ppc["bus"][ppc_index, SKSS]
    if net._options["ip"]:
        net.res_bus_sc["ip_ka"] = ppc["bus"][ppc_index, IP]
    if net._options["ith"]:
        net.res_bus_sc["ith_ka"] = ppc["bus"][ppc_index, ITH]

    # Export also equivalent rk, xk on the calculated bus
    net.res_bus_sc["rk_ohm"] = ppc["bus"][ppc_index, R_EQUIV_OHM]
    net.res_bus_sc["xk_ohm"] = ppc["bus"][ppc_index, X_EQUIV_OHM]

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



@@ -103,14 +121,14 @@

        f, t = branch_lookup["line"]
        minmax = np.maximum if case == "max" else np.minimum

        net.res_line_sc["ikss_ka"] = minmax(ppc["internal"]["branch_ikss_f"][f:t, ppc_index].real.reshape(-1, 1),
                                            ppc["internal"]["branch_ikss_t"][f:t, ppc_index].real.reshape(-1, 1))
        net.res_line_sc["ikss_ka"] = minmax(ppc["internal"]["branch_ikss_f"].iloc[f:t,:].loc[:, ppc_index].values.real.reshape(-1, 1),
                                            ppc["internal"]["branch_ikss_t"].iloc[f:t,:].loc[:, ppc_index].values.real.reshape(-1, 1))
        if net._options["ip"]:
            net.res_line_sc["ip_ka"] = minmax(ppc["internal"]["branch_ip_f"][f:t, ppc_index].real.reshape(-1, 1),
                                              ppc["internal"]["branch_ip_t"][f:t, ppc_index].real.reshape(-1, 1))
            net.res_line_sc["ip_ka"] = minmax(ppc["internal"]["branch_ip_f"].iloc[f:t,:].loc[:, ppc_index].values.real.reshape(-1, 1),
                                              ppc["internal"]["branch_ip_t"].iloc[f:t,:].loc[:, ppc_index].values.real.reshape(-1, 1))
        if net._options["ith"]:
            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))
            net.res_line_sc["ith_ka"] = minmax(ppc["internal"]["branch_ith_f"].iloc[f:t,:].loc[:, ppc_index].values.real.reshape(-1, 1),
                                               ppc["internal"]["branch_ith_t"].iloc[f:t,:].loc[:, ppc_index].values.real.reshape(-1, 1))

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

@@ -129,8 +147,8 @@


    if "trafo" in branch_lookup:
        f, t = branch_lookup["trafo"]
        net.res_trafo_sc["ikss_hv_ka"] = ppc["internal"]["branch_ikss_f"][f:t, ppc_index].real.reshape(-1, 1)
        net.res_trafo_sc["ikss_lv_ka"] = ppc["internal"]["branch_ikss_t"][f:t, ppc_index].real.reshape(-1, 1)
        net.res_trafo_sc["ikss_hv_ka"] = ppc["internal"]["branch_ikss_f"].iloc[f:t,:].loc[:, ppc_index].values.real.reshape(-1, 1)
        net.res_trafo_sc["ikss_lv_ka"] = ppc["internal"]["branch_ikss_t"].iloc[f:t,:].loc[:, ppc_index].values.real.reshape(-1, 1)


def _get_trafo3w_results(net, ppc):

@@ -157,6 +175,6 @@

        hv = int(f + (t - f) / 3)
        mv = int(f + 2 * (t - f) / 3)
        lv = t
        net.res_trafo3w_sc["ikss_hv_ka"] = ppc["internal"]["branch_ikss_f"][f:hv, ppc_index].real.reshape(-1, 1)
        net.res_trafo3w_sc["ikss_mv_ka"] = ppc["internal"]["branch_ikss_t"][hv:mv, ppc_index].real.reshape(-1, 1)
        net.res_trafo3w_sc["ikss_lv_ka"] = ppc["internal"]["branch_ikss_t"][mv:lv, ppc_index].real.reshape(-1, 1)

@@ -1349,7 +1349,8 @@

def create_gen(net, bus, p_mw, vm_pu=1., sn_mva=nan, name=None, index=None, max_q_mvar=nan,
               min_q_mvar=nan, min_p_mw=nan, max_p_mw=nan, min_vm_pu=nan, max_vm_pu=nan,
               scaling=1., type=None, slack=False, controllable=nan, vn_kv=nan,
               xdss_pu=nan, rdss_pu=nan, cos_phi=nan, in_service=True):
               xdss_pu=nan, rdss_ohm=nan, cos_phi=nan, pg_percent=nan, power_station_trafo=None,
               in_service=True):
    """
    Adds a generator to the network.


@@ -1388,9 +1389,13 @@


        **xdss_pu** (float, NaN) - Subtransient generator reactance for short-circuit calculation

        **rdss_pu** (float, NaN) - Subtransient generator resistance for short-circuit calculation
        **rdss_ohm** (float, NaN) - Subtransient generator resistance for short-circuit calculation

        **cos_phi** (float, NaN) - Rated cosine phi of the generator for short-circuit calculation
        
        **pg_percent** (float, NaN) - Rated pg (voltage control range) of the generator for short-circuit calculation
        
        **power_station_trafo** (int, None) - Index of the power station transformer for short-circuit calculation

        **in_service** (bool, True) - True for in_service or False for out of service


@@ -1444,18 +1449,14 @@

    _create_column_and_set_value(net, index, max_vm_pu, "max_vm_pu", "gen", default_val=2.)
    _create_column_and_set_value(net, index, min_vm_pu, "min_vm_pu", "gen", default_val=0.)

    # Short circuit calculation limits
    # Short circuit calculation variables
    _create_column_and_set_value(net, index, vn_kv, "vn_kv", "gen")
    _create_column_and_set_value(net, index, cos_phi, "cos_phi", "gen")

    if not isnan(xdss_pu):
        if "xdss_pu" not in net.gen.columns:
            net.gen.loc[:, "xdss_pu"] = pd.Series(dtype=float64)
        if "rdss_pu" not in net.gen.columns:
            net.gen.loc[:, "rdss_pu"] = pd.Series(dtype=float64)
        net.gen.at[index, "xdss_pu"] = float(xdss_pu)

    _create_column_and_set_value(net, index, rdss_pu, "rdss_pu", "gen")
    _create_column_and_set_value(net, index, cos_phi, "cos_phi", "gen") 
    _create_column_and_set_value(net, index, xdss_pu, "xdss_pu", "gen")
    _create_column_and_set_value(net, index, rdss_ohm, "rdss_ohm", "gen")
    _create_column_and_set_value(net, index, pg_percent, "pg_percent", "gen")
    _create_column_and_set_value(net, index, power_station_trafo,
                                 "power_station_trafo", "gen")

    return index


@@ -1463,7 +1464,8 @@

def create_gens(net, buses, p_mw, vm_pu=1., sn_mva=nan, name=None, index=None, max_q_mvar=None,
                min_q_mvar=None, min_p_mw=None, max_p_mw=None, min_vm_pu=None, max_vm_pu=None,
                scaling=1., type=None, slack=False, controllable=None, vn_kv=None,
                xdss_pu=None, rdss_pu=None, cos_phi=None, in_service=True, **kwargs):
                xdss_pu=None, rdss_ohm=None, cos_phi=None, pg_percent=None, power_station_trafo=None,
                in_service=True, **kwargs):
    """
    Adds generators to the specified buses network.


@@ -1507,11 +1509,16 @@

        **xdss_pu** (list of float, NaN) - Subtransient generator reactance for short-circuit \
            calculation

        **rdss_pu** (list of float, NaN) - Subtransient generator resistance for short-circuit \
        **rdss_ohm** (list of float, NaN) - Subtransient generator resistance for short-circuit \
            calculation

        **cos_phi** (list of float, NaN) - Rated cosine phi of the generator for short-circuit \
            calculation
            
        **pg_percent** (float, NaN) - Rated pg (voltage control range) of the generator for \
            short-circuit calculation
        
        **power_station_trafo** (int, None) - Index of the power station transformer for short-circuit calculation

        **in_service** (bool, True) - True for in_service or False for out of service


@@ -1558,7 +1565,9 @@

    _add_series_to_entries(entries, index, "vn_kv", vn_kv)
    _add_series_to_entries(entries, index, "cos_phi", cos_phi)
    _add_series_to_entries(entries, index, "xdss_pu", xdss_pu)
    _add_series_to_entries(entries, index, "rdss_pu", rdss_pu)
    _add_series_to_entries(entries, index, "rdss_ohm", rdss_ohm)
    _add_series_to_entries(entries, index, "pg_percent", pg_percent)
    _add_series_to_entries(entries, index, "power_station_trafo", power_station_trafo)
    _add_series_to_entries(entries, index, "controllable", controllable, dtyp=bool_,
                           default_val=False)


@@ -2219,7 +2228,8 @@

                                       tap_phase_shifter=False, in_service=True, name=None,
                                       vector_group=None, index=None, max_loading_percent=nan,
                                       parallel=1, df=1., vk0_percent=nan, vkr0_percent=nan,
                                       mag0_percent=nan, mag0_rx=nan, si0_hv_partial=nan, **kwargs):
                                       mag0_percent=nan, mag0_rx=nan, si0_hv_partial=nan,
                                       pt_percent=nan, **kwargs):
    """
    Creates a two-winding transformer in table net["trafo"].
    The trafo parameters are defined through the standard type library.

@@ -2298,6 +2308,8 @@


        **df** (float) - derating factor: maximal current of transformer in relation to nominal \
            current of transformer (from 0 to 1)
            
        **pt_percent** (float, NaN) - pt (shortcircuit voltage) of the transformer with OLTC for short-circuit calculation

        ** only considered in loadflow if calculate_voltage_angles = True


@@ -2351,6 +2363,7 @@

        _create_column_and_set_value(net, index, si0_hv_partial, "si0_hv_partial", "trafo")
        _create_column_and_set_value(net, index, vector_group, "vector_group", "trafo", dtyp=str,
                                     default_val=None)
    _create_column_and_set_value(net, index, pt_percent, "pt_percent", "trafo")   
    _create_column_and_set_value(net, index, max_loading_percent, "max_loading_percent", "trafo")

    return index

@@ -2364,7 +2377,7 @@

                                        vector_group=None, index=None, max_loading_percent=None,
                                        parallel=1, df=1., vk0_percent=None, vkr0_percent=None,
                                        mag0_percent=None, mag0_rx=None, si0_hv_partial=None,
                                        **kwargs):
                                        pt_percent=nan, **kwargs):
    """
    Creates several two-winding transformers in table net["trafo"].
    The trafo parameters are defined through the standard type library.

@@ -2445,6 +2458,8 @@

        **df** (float) - derating factor: maximal current of transformer in relation to nominal \
            current of transformer (from 0 to 1)

        **pt_percent** (float, NaN) - pt (shortcircuit voltage) of the transformer with OLTC for short-circuit calculation

        ** only considered in loadflow if calculate_voltage_angles = True

    OUTPUT:

@@ -2477,6 +2492,7 @@

    _add_series_to_entries(entries, index, "si0_hv_partial", si0_hv_partial)
    _add_series_to_entries(entries, index, "max_loading_percent", max_loading_percent)
    _add_series_to_entries(entries, index, "vector_group", vector_group, dtyp=str)
    _add_series_to_entries(entries, index, "pt_percent", pt_percent)

    _set_multiple_entries(net, "trafo", index, **entries, **kwargs)


@@ -2595,7 +2611,10 @@

                                         tap_step_percent=nan, tap_step_degree=nan, tap_pos=nan,
                                         tap_neutral=nan, tap_max=nan,
                                         tap_min=nan, name=None, in_service=True, index=None,
                                         max_loading_percent=nan, tap_at_star_point=False):
                                         max_loading_percent=nan, tap_at_star_point=False,
                                         vk0_hv_percent=nan, vk0_mv_percent=nan, vk0_lv_percent=nan,
                                         vkr0_hv_percent=nan, vkr0_mv_percent=nan, vkr0_lv_percent=nan,
                                         vector_group=None):
    """
    Adds a three-winding transformer in table net["trafo3w"].


@@ -2670,6 +2689,20 @@

        **The model currently only supports one tap-changer per 3W Transformer.

        **max_loading_percent (float)** - maximum current loading (only needed for OPF)
        
        **vk0_hv_percent** (float) - zero sequence short circuit voltage from high to medium voltage

        **vk0_mv_percent** (float) - zero sequence short circuit voltage from medium to low voltage

        **vk0_lv_percent** (float) - zero sequence short circuit voltage from high to low voltage

        **vkr0_hv_percent** (float) - zero sequence real part of short circuit voltage from high to medium voltage

        **vkr0_mv_percent** (float) - zero sequence real part of short circuit voltage from medium to low voltage

        **vkr0_lv_percent** (float) - zero sequence real part of short circuit voltage from high to low voltage
        
        **vector_group** (list of String) - Vector group of the transformer3w

    OUTPUT:
        **trafo_id** - The unique trafo_id of the created 3W transformer

@@ -2698,12 +2731,15 @@

               "vkr_hv_percent", "vkr_mv_percent", "vkr_lv_percent", "pfe_kw", "i0_percent",
               "shift_mv_degree", "shift_lv_degree", "tap_side", "tap_step_percent",
               "tap_step_degree", "tap_pos", "tap_neutral", "tap_max", "tap_min", "in_service",
               "name", "std_type", "tap_at_star_point"]
               "name", "std_type", "tap_at_star_point", "vk0_hv_percent", "vk0_mv_percent", "vk0_lv_percent",
               "vkr0_hv_percent", "vkr0_mv_percent", "vkr0_lv_percent", "vector_group"]
    values = [lv_bus, mv_bus, hv_bus, vn_hv_kv, vn_mv_kv, vn_lv_kv, sn_hv_mva, sn_mv_mva, sn_lv_mva,
              vk_hv_percent, vk_mv_percent, vk_lv_percent, vkr_hv_percent, vkr_mv_percent,
              vkr_lv_percent, pfe_kw, i0_percent, shift_mv_degree, shift_lv_degree, tap_side,
              tap_step_percent, tap_step_degree, tap_pos, tap_neutral, tap_max, tap_min,
              bool(in_service), name, None, tap_at_star_point]
              bool(in_service), name, None, tap_at_star_point,
              vk0_hv_percent, vk0_mv_percent, vk0_lv_percent,
              vkr0_hv_percent, vkr0_mv_percent, vkr0_lv_percent, vector_group]

    _set_entries(net, "trafo3w", index, **dict(zip(columns, values)))


@@ -2720,7 +2756,10 @@

                                          tap_step_percent=nan, tap_step_degree=nan, tap_pos=nan,
                                          tap_neutral=nan, tap_max=nan, tap_min=nan, name=None,
                                          in_service=True, index=None, max_loading_percent=None,
                                          tap_at_star_point=False, **kwargs):
                                          tap_at_star_point=False,
                                          vk0_hv_percent=nan, vk0_mv_percent=nan, vk0_lv_percent=nan,
                                          vkr0_hv_percent=nan, vkr0_mv_percent=nan, vkr0_lv_percent=nan,
                                          vector_group=None, **kwargs):
    """
    Adds a three-winding transformer in table net["trafo3w"].


@@ -2799,6 +2838,20 @@


        **max_loading_percent (float)** - maximum current loading (only needed for OPF)

        **vk0_hv_percent** (float) - zero sequence short circuit voltage from high to medium voltage

        **vk0_mv_percent** (float) - zero sequence short circuit voltage from medium to low voltage

        **vk0_lv_percent** (float) - zero sequence short circuit voltage from high to low voltage

        **vkr0_hv_percent** (float) - zero sequence real part of short circuit voltage from high to medium voltage

        **vkr0_mv_percent** (float) - zero sequence real part of short circuit voltage from medium to low voltage

        **vkr0_lv_percent** (float) - zero sequence real part of short circuit voltage from high to low voltage
        
        **vector_group** (list of String) - Vector group of the transformer3w

    OUTPUT:
        **trafo_id** - List of trafo_ids of the created 3W transformers


@@ -2836,7 +2889,11 @@

               "tap_step_degree": tap_step_degree, "tap_pos": tp_pos, "tap_neutral": tp_neutral,
               "tap_max": tap_max, "tap_min": tap_min,
               "in_service": array(in_service).astype(bool_), "name": name,
               "tap_at_star_point": array(tap_at_star_point).astype(bool_), "std_type": None}
               "tap_at_star_point": array(tap_at_star_point).astype(bool_), "std_type": None,
               "vk0_hv_percent":vk0_hv_percent, "vk0_mv_percent":vk0_mv_percent,
               "vk0_lv_percent":vk0_lv_percent, "vkr0_hv_percent":vkr0_hv_percent,
               "vkr0_mv_percent":vkr0_mv_percent, "vkr0_lv_percent":vkr0_lv_percent,
               "vector_group": vector_group}

    _add_series_to_entries(entries, index, "max_loading_percent", max_loading_percent)


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

import pandapower.auxiliary as aux
from pandapower.build_branch import _switch_branches, _branches_with_oos_buses, _build_branch_ppc
from pandapower.build_bus import _build_bus_ppc, _calc_pq_elements_and_add_on_ppc, \
_calc_shunts_and_add_on_ppc, _add_gen_impedances_ppc, _add_motor_impedances_ppc
_calc_shunts_and_add_on_ppc, _add_ext_grid_sc_impedance, _add_motor_impedances_ppc
from pandapower.build_gen import _build_gen_ppc, _check_voltage_setpoints_at_same_bus, \
    _check_voltage_angles_at_same_bus, _check_for_reference_bus
from pandapower.opf.make_objective import _make_objective

@@ -28,7 +28,7 @@

           and fill it in the branch matrix.
           Order: 1st: Line values, 2nd: Trafo values
        5. if opf: make opf objective (gencost)
        6. convert internal ppci format for pypower powerflow / 
        6. convert internal ppci format for pypower powerflow /
        opf without out of service elements and rearanged buses

    INPUT:

@@ -37,7 +37,7 @@

        ( 0 - Zero Sequence
          1 - Positive Sequence
          2 - Negative Sequence
        ) 
        )

    OUTPUT:
        **ppc** - The simple matpower format network. Which consists of:

@@ -56,7 +56,7 @@

                              , "gen_is": np.array([], dtype=bool)
                              }
        **ppci** - The "internal" pypower format network for PF calculations
        

    """
    # select elements in service (time consuming, so we do it once)
    net["_is_elements"] = aux._select_is_elements_numba(net, sequence=sequence)

@@ -77,12 +77,13 @@

        # Calculates ppc0 branch impedances from branch elements
        _build_branch_ppc_zero(net, ppc)
    else:
        # Calculates ppc1/ppc2 branch impedances from branch elements  
        # Calculates ppc1/ppc2 branch impedances from branch elements
        _build_branch_ppc(net, ppc)

    # Adds P and Q for loads / sgens in ppc['bus'] (PQ nodes)
    if mode == "sc":
        _add_gen_impedances_ppc(net, ppc)
        _add_ext_grid_sc_impedance(net, ppc)
        # Generator impedance are seperately added in sc module
        _add_motor_impedances_ppc(net, ppc)
    else:
        _calc_pq_elements_and_add_on_ppc(net, ppc, sequence=sequence)

@@ -119,7 +120,7 @@


    aux._replace_nans_with_default_limits(net, ppc)

    # generates "internal" ppci format (for powerflow calc) 
    # generates "internal" ppci format (for powerflow calc)
    # from "external" ppc format and updates the bus lookup
    # Note: Also reorders buses and gens in ppc
    ppci = _ppc2ppci(ppc, net)

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

    lookup = net._pd2ppc_lookups["branch"]
    mode = net._options["mode"]
    ppc["branch"] = np.zeros(shape=(length, branch_cols), dtype=np.complex128)
    # Check if this should be moved to somewhere else
    if mode == "sc":
        from pandapower.shortcircuit.idx_brch import branch_cols_sc
        branch_sc = np.empty(shape=(length, branch_cols_sc), dtype=float)

@@ -132,7 +133,7 @@

    length_km = line["length_km"].values
    parallel = line["parallel"].values
    base_kv = ppc["bus"][from_bus, BASE_KV]
    baseR = np.square(base_kv) / (3 * net.sn_mva) if mode == "pf_3ph" else   np.square(base_kv) / net.sn_mva 
    baseR = np.square(base_kv) / (3 * net.sn_mva) if mode == "pf_3ph" else   np.square(base_kv) / net.sn_mva
    branch[f:t, F_BUS] = from_bus
    branch[f:t, T_BUS] = to_bus
    branch[f:t, BR_R] = line["r_ohm_per_km"].values * length_km / baseR / parallel

@@ -207,7 +208,7 @@

        return trafo_df[par].values


def _calc_branch_values_from_trafo_df(net, ppc, trafo_df=None):
def _calc_branch_values_from_trafo_df(net, ppc, trafo_df=None, seq=1):
    """
    Calculates the MAT/PYPOWER-branch-attributes from the pandapower trafo dataframe.


@@ -250,15 +251,15 @@

    vn_trafo_hv, vn_trafo_lv, shift = _calc_tap_from_dataframe(net, trafo_df)
    ratio = _calc_nominal_ratio_from_dataframe(ppc, trafo_df, vn_trafo_hv, vn_trafo_lv,
                                               bus_lookup)
    r, x, y = _calc_r_x_y_from_dataframe(net, trafo_df, vn_trafo_lv, vn_lv, ppc)
    r, x, y = _calc_r_x_y_from_dataframe(net, trafo_df, vn_trafo_lv, vn_lv, ppc, seq=seq)
    return r, x, y, ratio, shift


def _calc_r_x_y_from_dataframe(net, trafo_df, vn_trafo_lv, vn_lv, ppc):
def _calc_r_x_y_from_dataframe(net, trafo_df, vn_trafo_lv, vn_lv, ppc, seq=1):
    mode = net["_options"]["mode"]
    trafo_model = net["_options"]["trafo_model"]

    r, x = _calc_r_x_from_dataframe(mode,trafo_df, vn_lv, vn_trafo_lv, net.sn_mva)
    r, x = _calc_r_x_from_dataframe(mode,trafo_df, vn_lv, vn_trafo_lv, net.sn_mva, seq=seq)
    if mode == "sc":
        y = 0
        if isinstance(trafo_df, pd.DataFrame):  # 2w trafo is dataframe, 3w trafo is dict

@@ -271,6 +272,7 @@

            x *= kt
    else:
        y = _calc_y_from_dataframe(mode,trafo_df, vn_lv, vn_trafo_lv, net.sn_mva)

    if trafo_model == "pi":
        return r, x, y
    elif trafo_model == "t":

@@ -311,14 +313,14 @@

        **subsceptance** (1d array, np.complex128) - The subsceptance in pu in
        the form (-b_img, -b_real)
    """
   

    baseR = np.square(vn_lv) / (3*sn_mva) if mode == 'pf_3ph' else np.square(vn_lv) / sn_mva
    vn_lv_kv = get_trafo_values(trafo_df, "vn_lv_kv")
    pfe = get_trafo_values(trafo_df, "pfe_kw") * 1e-3
    parallel = get_trafo_values(trafo_df, "parallel")

    ### Calculate subsceptance ###
   

    vnl_squared = (vn_lv_kv ** 2)/3 if mode == 'pf_3ph' else  vn_lv_kv **2
    b_real = pfe / vnl_squared * baseR
    i0 = get_trafo_values(trafo_df, "i0_percent")

@@ -403,15 +405,21 @@

    array[mask] = value
    return array

def _calc_r_x_from_dataframe(mode,trafo_df, vn_lv, vn_trafo_lv, sn_mva):
def _calc_r_x_from_dataframe(mode,trafo_df, vn_lv, vn_trafo_lv, sn_mva, seq=1):
    """
    Calculates (Vectorized) the resitance and reactance according to the
    transformer values

    """
    parallel = get_trafo_values(trafo_df, "parallel")
    vk_percent = get_trafo_values(trafo_df, "vk_percent")
    vkr_percent = get_trafo_values(trafo_df, "vkr_percent")
    if seq == 1:
        vk_percent = get_trafo_values(trafo_df, "vk_percent")
        vkr_percent = get_trafo_values(trafo_df, "vkr_percent")
    elif seq == 0:
        vk_percent = get_trafo_values(trafo_df, "vk0_percent")
        vkr_percent = get_trafo_values(trafo_df, "vkr0_percent")
    else:
        raise UserWarning("Unsupported sequence")
    tap_lv = np.square(vn_trafo_lv / vn_lv) * (3* sn_mva)  if mode == 'pf_3ph' else\
    np.square(vn_trafo_lv / vn_lv) * sn_mva  # adjust for low voltage side voltage converter
    sn_trafo_mva = get_trafo_values(trafo_df, "sn_mva")

@@ -791,7 +799,7 @@

    _is_elements["line"] = net["line"][net["line"]["in_service"].values.astype(bool)]


def _trafo_df_from_trafo3w(net):
def _trafo_df_from_trafo3w(net, seq=1):
    nr_trafos = len(net["trafo3w"])
    trafo2 = dict()
    sides = ["hv", "mv", "lv"]

@@ -799,7 +807,14 @@

    loss_side = net._options["trafo3w_losses"].lower()
    nr_trafos = len(net["trafo3w"])
    t3 = net["trafo3w"]
    _calculate_sc_voltages_of_equivalent_transformers(t3, trafo2, mode)
    if seq==1:
        _calculate_sc_voltages_of_equivalent_transformers(t3, trafo2, mode)
    elif seq==0:
        if mode != "sc":
            raise NotImplementedError("0 seq impedance calculation only implemented for short-circuit calculation!")
        _calculate_sc_voltages_of_equivalent_transformers_zero_sequence(t3, trafo2,)
    else:
        raise UserWarning("Unsupported sequence for trafo3w convertion")
    _calculate_3w_tap_changers(t3, trafo2, sides)
    zeros = np.zeros(len(net.trafo3w))
    aux_buses = net._pd2ppc_lookups["aux"]["trafo3w"]

@@ -842,6 +857,32 @@

842	857		t2["sn_mva"] = {"hv": sn[0, :], "mv": sn[1, :], "lv": sn[2, :]}
843	858
844	859
	860	+	def _calculate_sc_voltages_of_equivalent_transformers_zero_sequence(t3, t2):
	861	+	vk_3w = np.stack([t3.vk_hv_percent.values, t3.vk_mv_percent.values, t3.vk_lv_percent.values])
	862	+	vkr_3w = np.stack([t3.vkr_hv_percent.values, t3.vkr_mv_percent.values, t3.vkr_lv_percent.values])
	863	+	vk0_3w = np.stack([t3.vk0_hv_percent.values, t3.vk0_mv_percent.values, t3.vk0_lv_percent.values])
	864	+	vkr0_3w = np.stack([t3.vkr0_hv_percent.values, t3.vkr0_mv_percent.values, t3.vkr0_lv_percent.values])
	865	+	sn = np.stack([t3.sn_hv_mva.values, t3.sn_mv_mva.values, t3.sn_lv_mva.values])
	866	+
	867	+	vk0_2w_delta = z_br_to_bus_vector(vk0_3w, sn)
	868	+	vkr0_2w_delta = z_br_to_bus_vector(vkr0_3w, sn)
	869	+
	870	+	# Only for "sc", calculated with positive sequence value
	871	+	kt = _transformer_correction_factor(vk_3w, vkr_3w, sn, 1.1)
	872	+	vk0_2w_delta *= kt
	873	+	vkr0_2w_delta *= kt
	874	+
	875	+	vki0_2w_delta = np.sqrt(vk0_2w_delta 2 - vkr0_2w_delta 2)
	876	+	vkr0_2w = wye_delta_vector(vkr0_2w_delta, sn)
	877	+	vki0_2w = wye_delta_vector(vki0_2w_delta, sn)
	878	+	vk0_2w = np.sign(vki0_2w) * np.sqrt(vki0_2w 2 + vkr0_2w 2)
	879	+	if np.any(vk0_2w == 0):
	880	+	raise UserWarning("Equivalent transformer with zero impedance!")
	881	+	t2["vk0_percent"] = {"hv": vk0_2w[0, :], "mv": vk0_2w[1, :], "lv": vk0_2w[2, :]}
	882	+	t2["vkr0_percent"] = {"hv": vkr0_2w[0, :], "mv": vkr0_2w[1, :], "lv": vkr0_2w[2, :]}
	883	+	t2["sn_mva"] = {"hv": sn[0, :], "mv": sn[1, :], "lv": sn[2, :]}
	884	+
	885	+
845	886		def z_br_to_bus_vector(z, sn):
846	887		return sn[0, :] * np.array([z[0, :] / sn[[0, 1], :].min(axis=0), z[1, :] /
847	888		sn[[1, 2], :].min(axis=0), z[2, :] / sn[[0, 2], :].min(axis=0)])

@@ -12,9 +12,10 @@

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

from pandapower.shortcircuit.idx_bus import KAPPA, R_EQUIV, X_EQUIV
from pandapower.shortcircuit.idx_bus import KAPPA, R_EQUIV, X_EQUIV, GS_P, BS_P, K_G
from pandapower.shortcircuit.impedance import _calc_ybus, _calc_zbus, _calc_rx


def _add_kappa_to_ppc(net, ppc):
    if not net._options["kappa"]:
        return

@@ -30,9 +31,11 @@

        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

@@ -40,11 +43,17 @@

        fc = 24
    else:
        raise ValueError("Frequency has to be 50 Hz or 60 Hz according to the standard")

    ppc_c = copy.deepcopy(ppc)
    ppc_c["branch"][:, BR_X] *= fc / net.f_hz

    zero_conductance = np.where(ppc["bus"][:,GS] == 0)
    ppc["bus"][zero_conductance, BS] *= net.f_hz / fc
    # Generator peak admittance application
    ppc_c["bus"][np.ix_(~np.isnan(ppc_c["bus"][:, GS_P]), [BS, GS])] =\
        ppc_c["bus"][np.ix_(~np.isnan(ppc_c["bus"][:, GS_P]), [BS_P, GS_P])]

    # # TODO: Check this
    # zero_conductance = np.where(ppc["bus"][:,GS] == 0)
    # ppc_c["bus"][zero_conductance, BS] *= net.f_hz / fc

    conductance = np.where(ppc["bus"][:,GS] != 0)
    z_shunt = 1 / (ppc_c["bus"][conductance, GS] + 1j * ppc_c["bus"][conductance, BS])

@@ -59,10 +68,11 @@

        # Factorization Ybus once
        ppc_c["internal"]["ybus_fact"] = factorized(ppc_c["internal"]["Ybus"])

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


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

@@ -90,6 +100,7 @@

    rx_equiv = ppc["bus"][:, R_EQUIV] / ppc["bus"][:, X_EQUIV]
    return np.clip(kappa_korr * _kappa(rx_equiv), 1, kappa_max)


def nxgraph_from_ppc(net, ppc):
    bus_lookup = net._pd2ppc_lookups["bus"]
    mg = nx.MultiGraph()

@@ -100,7 +111,15 @@

    vs_buses_pp = list(set(net["ext_grid"][net._is_elements["ext_grid"]].bus.values) |
                       set(net["gen"][net._is_elements["gen"]].bus))
    vs_buses = bus_lookup[vs_buses_pp]
    z = 1 / (ppc["bus"][vs_buses, GS] + ppc["bus"][vs_buses, BS] * 1j)
    mg.add_edges_from(("earth", int(bus), {"r": z.real, "x": z.imag})
                        for bus, z in zip(vs_buses, z))
    gen_vs_buses = vs_buses[~np.isnan(ppc["bus"][vs_buses, K_G])]
    non_gen_vs_buses = vs_buses[np.isnan(ppc["bus"][vs_buses, K_G])]
    if np.size(gen_vs_buses):
        z = 1 / (ppc["bus"][gen_vs_buses, GS_P] + ppc["bus"][gen_vs_buses, BS_P] * 1j)
        mg.add_edges_from(("earth", int(bus), {"r": z.real, "x": z.imag})
                            for bus, z in zip(gen_vs_buses, z))

    if np.size(non_gen_vs_buses):
        z = 1 / (ppc["bus"][non_gen_vs_buses, GS] + ppc["bus"][non_gen_vs_buses, BS] * 1j)
        mg.add_edges_from(("earth", int(bus), {"r": z.real, "x": z.imag})
                            for bus, z in zip(non_gen_vs_buses, z))
    return mg

@@ -10,103 +10,75 @@

from scipy.sparse.linalg import inv as inv_sparse
from scipy.linalg import inv


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, bus):
def _calc_rx(net, ppci, bus_idx):
    # 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"][bus_idx, BASE_KV]) / ppc["baseMVA"]
        base_r = np.square(ppci["bus"][bus_idx, BASE_KV]) / ppci["baseMVA"]
        fault_impedance = (r_fault + x_fault * 1j) / base_r
    else:
        fault_impedance = 0 + 0j

    if net["_options"]["inverse_y"]:
        Zbus = ppc["internal"]["Zbus"]
        Zbus = ppci["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
        z_equiv = _calc_zbus_diag(net, ppci, bus_idx) + fault_impedance
    ppci["bus"][bus_idx, R_EQUIV] = z_equiv.real
    ppci["bus"][bus_idx, X_EQUIV] = z_equiv.imag


def _calc_ybus(ppc):
    Ybus, Yf, Yt = makeYbus(ppc["baseMVA"], ppc["bus"],  ppc["branch"])
def _calc_ybus(ppci):
    Ybus, Yf, Yt = makeYbus(ppci["baseMVA"], ppci["bus"], ppci["branch"])
    if np.isnan(Ybus.data).any():
        raise ValueError("nan value detected in Ybus matrix - check calculation parameters for nan values")
    ppc["internal"]["Yf"] = Yf
    ppc["internal"]["Yt"] = Yt
    ppc["internal"]["Ybus"] = Ybus
    ppci["internal"]["Yf"] = Yf
    ppci["internal"]["Yt"] = Yt
    ppci["internal"]["Ybus"] = Ybus


def _calc_zbus(net, ppc):
def _calc_zbus(net, ppci):
    try:
        Ybus = ppc["internal"]["Ybus"]
        Ybus = ppci["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()
                ppci["internal"]["Zbus"] = inv_sparse(Ybus).toarray()
        else:
            ppc["internal"]["Zbus"] = inv(Ybus.toarray())
            ppci["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]
def _calc_zbus_diag(net, ppci, bus_idx=None):
    ybus_fact = ppci["internal"]["ybus_fact"]
    n_ppci_bus = ppci["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)
    if bus_idx is None:
        diagZ = np.zeros(n_ppci_bus, dtype=np.complex)
        for i in range(n_ppci_bus):
            b = np.zeros(n_ppci_bus, dtype=np.complex)
            b[i] = 1 + 0j
            diagZ[i] = ybus_fact(b)[i]
        ppc["internal"]["diagZ"] = diagZ
        ppci["internal"]["diagZ"] = diagZ
        return diagZ
    else:
        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]
        diagZ = np.zeros(bus_idx.shape[0], dtype=np.complex)
        for ix, b in enumerate(bus_idx):
            rhs = np.zeros(n_ppci_bus, dtype=np.complex)
            rhs[b] = 1 + 0j
            diagZ[ix] = ybus_fact(rhs)[b]
        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

@@ -201,11 +201,16 @@

            # 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
            # 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
            if "br_res_ks_ppci_bus" in result["internal"]:
                br_res_ks_ppci_bus = np.r_[result["internal"]["br_res_ks_ppci_bus"], [-1]]
            else:
                br_res_ks_ppci_bus = np.r_[np.arange(value.shape[1]), [-1]]
            ppc['internal'][key] = pd.DataFrame(data=update_matrix, columns=br_res_ks_ppci_bus)
        else:
            ppc["internal"][key] = value


Files		•	•	Coverage
pandapower	15,810	13,903	1,907	87.94%
setup.py	14	0	14	0.00%
Project Totals (163 files)	15,824	13,903	1,921	87.86%

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

328	251		ikss_all_t = abs(ikss1_all_t)
329	252
330	253		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]
	254	+	ppci["internal"]["branch_ikss_f"] = ikss_all_f / baseI[fb, None]
	255	+	ppci["internal"]["branch_ikss_t"] = ikss_all_t / baseI[tb, None]
333	256		else:
334	257		ikss_all_f[abs(ikss_all_f) < 1e-10] = np.nan
335	258		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])
	259	+	ppci["branch"][:, IKSS_F] = np.nan_to_num(minmax(ikss_all_f, axis=1) / baseI[fb])
	260	+	ppci["branch"][:, IKSS_T] = np.nan_to_num(minmax(ikss_all_t, axis=1) / baseI[tb])
338	261
339	262		if net._options["ip"]:
340		-	kappa = ppc["bus"][:, KAPPA]
	263	+	kappa = ppci["bus"][:, KAPPA]
341	264		if current_sources:
342	265		ip_all_f = np.sqrt(2) * (ikss1_all_f * kappa[bus_idx] + ikss2_all_f)
343	266		ip_all_t = np.sqrt(2) * (ikss1_all_t * kappa[bus_idx] + ikss2_all_t)

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,
	581	+	kappa_method, return_all_currents,
581	582		inverse_y):
582	583		"""
583	584		creates dictionary for pf, opf and short circuit calculations from input parameters.

532	527		c = 1.1
533	528		if not "s_sc_%s_mva" % case in eg:
534	529		raise ValueError("short circuit apparent power s_sc_%s_mva needs to be specified for "% case +
535		-	"external grid \n Try: net.ext_grid['s_sc_max_mva'] = 1000" )
	530	+	"external grid \n Try: net.ext_grid['s_sc_max_mva'] = 1000" )
536	531		s_sc = eg["s_sc_%s_mva" % case].values/ppc['baseMVA']
537	532		if not "rx_%s" % case in eg:
538	533		raise ValueError("short circuit R/X rate rx_%s needs to be specified for external grid \n Try: net.ext_grid['rx_max'] = 0.1" %
539		-	case)
	534	+	case)
540	535		rx = eg["rx_%s" % case].values
541	536
542	537		z_grid = c / s_sc

1349	1349		def create_gen(net, bus, p_mw, vm_pu=1., sn_mva=nan, name=None, index=None, max_q_mvar=nan,
1350	1350		min_q_mvar=nan, min_p_mw=nan, max_p_mw=nan, min_vm_pu=nan, max_vm_pu=nan,
1351	1351		scaling=1., type=None, slack=False, controllable=nan, vn_kv=nan,
1352		-	xdss_pu=nan, rdss_pu=nan, cos_phi=nan, in_service=True):
	1352	+	xdss_pu=nan, rdss_ohm=nan, cos_phi=nan, pg_percent=nan, power_station_trafo=None,
	1353	+	in_service=True):
1353	1354		"""
1354	1355		Adds a generator to the network.
1355	1356

1388	1389
1389	1390		xdss_pu (float, NaN) - Subtransient generator reactance for short-circuit calculation
1390	1391
1391		-	rdss_pu (float, NaN) - Subtransient generator resistance for short-circuit calculation
	1392	+	rdss_ohm (float, NaN) - Subtransient generator resistance for short-circuit calculation
1392	1393
1393	1394		cos_phi (float, NaN) - Rated cosine phi of the generator for short-circuit calculation
	1395	+
	1396	+	pg_percent (float, NaN) - Rated pg (voltage control range) of the generator for short-circuit calculation
	1397	+
	1398	+	power_station_trafo (int, None) - Index of the power station transformer for short-circuit calculation
1394	1399
1395	1400		in_service (bool, True) - True for in_service or False for out of service
1396	1401

11	11		#from pandapower.pd2ppc import _ppc2ppci, _init_ppc
12	12		from pandapower.pypower.idx_brch import BR_B, BR_R, BR_X, F_BUS, T_BUS, branch_cols, BR_STATUS, SHIFT, TAP
13	13		from pandapower.pypower.idx_bus import BASE_KV, BS, GS
14		-	from pandapower.build_branch import _calc_tap_from_dataframe, _transformer_correction_factor, _calc_nominal_ratio_from_dataframe
	14	+	from pandapower.build_branch import _calc_tap_from_dataframe, _transformer_correction_factor, _calc_nominal_ratio_from_dataframe,\
	15	+	get_trafo_values, _trafo_df_from_trafo3w, _calc_branch_values_from_trafo_df
15	16		from pandapower.build_branch import _switch_branches, _branches_with_oos_buses, _initialize_branch_lookup, _end_temperature_correction_factor
16	17
17	18		def _pd2ppc_zero(net, sequence=0):

72	73		branch_sc.fill(np.nan)
73	74		ppc["branch"] = np.hstack((ppc["branch"], branch_sc ))
74	75		ppc["branch"][:, :13] = np.array([0, 0, 0, 0, 0, 250, 250, 250, 1, 0, 1, -360, 360])
	76	+
75	77		_add_line_sc_impedance_zero(net, ppc)
76	78		_add_trafo_sc_impedance_zero(net, ppc)
77		-	if "trafo3w" in lookup:
78		-	raise NotImplementedError("Three winding transformers are not implemented for unbalanced calculations")
	79	+	if mode == "sc":
	80	+	_add_trafo3w_sc_impedance_zero(net, ppc)
	81	+	else:
	82	+	if "trafo3w" in lookup:
	83	+	raise NotImplementedError("Three winding transformers are not implemented for unbalanced calculations")
79	84
80	85
81	86		def _add_trafo_sc_impedance_zero(net, ppc, trafo_df=None):

116	121		# Just put pos seq parameter if zero seq parameter is zero
117	122		if not "vkr0_percent" in trafos:
118	123		raise ValueError("Real part of short circuit voltage Vk0(Real) needs to be specified for transformer \
119		-	modelling \n Try : net.trafo[\"vkr0_percent\"] = net.trafo[\"vkr_percent\"]" )
	124	+	modelling \n Try : net.trafo[\"vkr0_percent\"] = net.trafo[\"vkr_percent\"]" )
120	125		vkr0_percent = trafos["vkr0_percent"].values.astype(float) if \
121	126		trafos["vkr0_percent"].values.astype(float).all() != 0. else \
122	127		trafos["vkr_percent"].values.astype(float)

129	134		# and mag0_percent = 10 ... 100 Zm0/ Zsc0
130	135		# --pg 50 DigSilent Power Factory Transformer manual
131	136		raise ValueError("Magnetizing impedance to vk0 ratio needs to be specified for transformer \
132		-	modelling \n Try : net.trafo[\"mag0_percent\"] = 100" )
	137	+	modelling \n Try : net.trafo[\"mag0_percent\"] = 100" )
133	138		mag0_ratio = trafos.mag0_percent.values.astype(float)
134	139		if not "mag0_rx" in trafos:
135	140		raise ValueError("Magnetizing impedance R/X ratio needs to be specified for transformer \
136		-	modelling \n Try : net.trafo[\"mag0_rx\"] = 0 " )
	141	+	modelling \n Try : net.trafo[\"mag0_rx\"] = 0 " )
137	142		mag0_rx = trafos["mag0_rx"].values.astype(float)
138	143		if not "si0_hv_partial" in trafos:
139	144		raise ValueError("Zero sequence short circuit impedance partition towards HV side needs to be specified for transformer \
140		-	modelling \n Try : net.trafo[\"si0_hv_partial\"] = 0.9 " )
	145	+	modelling \n Try : net.trafo[\"si0_hv_partial\"] = 0.9 " )
141	146		si0_hv_partial = trafos.si0_hv_partial.values.astype(float)
142	147		parallel = trafos.parallel.values.astype(float)
143	148		in_service = trafos["in_service"].astype(int)

208	213		YAB = 1 / zc.astype(complex)
209	214		YAN = 1 / za.astype(complex)
210	215		YBN = 1 / zb.astype(complex)
211		-
	216	+
212	217		# YAB_AN = (zc + za) /(zc * za).astype(complex) # Series conn YAB and YAN
213	218		# YAB_BN = (zc + zb) / (zc * zb).astype(complex) # Series conn YAB and YBN
214	219

243	248
244	249		gs_all = np.hstack([gs_all, y.real * in_service])
245	250		bs_all = np.hstack([bs_all, y.imag * in_service])
246		-
	251	+
247	252
248	253		elif vector_group == "YNyn":
249	254		ppc["branch"][ppc_idx, BR_STATUS] = in_service

296	301		ppc["bus"][buses, BS] += bs
297	302		del net.trafo["_ppc_idx"]
298	303
	304	+
299	305		def _add_ext_grid_sc_impedance_zero(net, ppc):
300	306		mode = net["_options"]["mode"]
301	307

326	332
327	333		z_grid = c / s_sc
328	334		if mode == 'pf_3ph':
329		-	z_grid = c / (s_sc/3)
	335	+	z_grid = c / (s_sc/3)
330	336		x_grid = z_grid / np.sqrt(rx ** 2 + 1)
331	337		r_grid = rx * x_grid
332	338		eg["r"] = r_grid

348	354
349	355		def _add_line_sc_impedance_zero(net, ppc):
350	356		branch_lookup = net["_pd2ppc_lookups"]["branch"]
351		-	mode = net["_options"]["mode"]
	357	+	mode = net["_options"]["mode"]
352	358		if not "line" in branch_lookup:
353	359		return
354	360		line = net["line"]

360	366		tb = bus_lookup[line["to_bus"].values]
361	367		baseR = np.square(ppc["bus"][fb, BASE_KV]) / net.sn_mva
362	368		if mode == 'pf_3ph':
363		-	baseR = np.square(ppc["bus"][fb, BASE_KV]) / (3*net.sn_mva)
	369	+	baseR = np.square(ppc["bus"][fb, BASE_KV]) / (3*net.sn_mva)
364	370		f, t = branch_lookup["line"]
365	371		# line zero sequence impedance
366	372		ppc["branch"][f:t, F_BUS] = fb

11	11		from pandapower.pypower.idx_bus import BASE_KV
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		-	from pandapower.shortcircuit.idx_bus import C_MIN, C_MAX, KAPPA, R_EQUIV, IKSS1, IP, ITH, X_EQUIV, IKSS2, IKCV, M
	14	+	from pandapower.shortcircuit.idx_bus import C_MIN, C_MAX, KAPPA, R_EQUIV, IKSS1, IP, ITH,\
	15	+	X_EQUIV, IKSS2, IKCV, M, R_EQUIV_OHM, X_EQUIV_OHM, V_G, K_SG, SKSS
15	16		from pandapower.shortcircuit.impedance import _calc_zbus_diag
16	17
17	18		from pandapower.pypower.pfsoln import pfsoln as pfsoln_pypower
18	19		from pandapower.pf.ppci_variables import _get_pf_variables_from_ppci
19	20
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)
27	21
	22	+	def _calc_ikss(net, ppci, bus_idx):
28	23		fault = net._options["fault"]
29	24		case = net._options["case"]
30		-	c = ppc["bus"][bus_idx, C_MIN] if case == "min" else ppc["bus"][bus_idx, C_MAX]
31		-	ppc["internal"]["baseI"] = ppc["bus"][:, BASE_KV] * np.sqrt(3) / ppc["baseMVA"]
	25	+	c = ppci["bus"][bus_idx, C_MIN] if case == "min" else ppci["bus"][bus_idx, C_MAX]
	26	+	ppci["internal"]["baseI"] = ppci["bus"][:, BASE_KV] * np.sqrt(3) / ppci["baseMVA"]
	27	+
	28	+	# Only for test, should correspondant to PF result
	29	+	baseZ = ppci["bus"][bus_idx, BASE_KV] ** 2 / ppci["baseMVA"]
	30	+	ppci["bus"][bus_idx, R_EQUIV_OHM] = baseZ * ppci["bus"][bus_idx, R_EQUIV]
	31	+	ppci["bus"][bus_idx, X_EQUIV_OHM] = baseZ * ppci["bus"][bus_idx, X_EQUIV]
32	32
33		-	z_equiv = abs(ppc["bus"][bus_idx, R_EQUIV] + ppc["bus"][bus_idx, X_EQUIV] * 1j)
	33	+	z_equiv = abs(ppci["bus"][bus_idx, R_EQUIV] + ppci["bus"][bus_idx, X_EQUIV] * 1j)
34	34		if fault == "3ph":
35		-	ppc["bus"][bus_idx, IKSS1] = c / z_equiv / ppc["bus"][bus_idx, BASE_KV] / np.sqrt(3) * ppc["baseMVA"]
	35	+	ppci["bus"][bus_idx, IKSS1] = c / z_equiv / ppci["bus"][bus_idx, BASE_KV] / np.sqrt(3) * ppci["baseMVA"]
36	36		elif fault == "2ph":
37		-	ppc["bus"][bus_idx, IKSS1] = c / z_equiv / ppc["bus"][bus_idx, BASE_KV] / 2 * ppc["baseMVA"]
38		-	_current_source_current(net, ppc)
	37	+	ppci["bus"][bus_idx, IKSS1] = c / z_equiv / ppci["bus"][bus_idx, BASE_KV] / 2 * ppci["baseMVA"]
39	38
	39	+	if fault == "3ph":
	40	+	ppci["bus"][bus_idx, SKSS] = np.sqrt(3) * ppci["bus"][bus_idx, IKSS1] * ppci["bus"][bus_idx, BASE_KV]
	41	+	elif fault == "2ph":
	42	+	ppci["bus"][bus_idx, SKSS] = ppci["bus"][bus_idx, IKSS1] * ppci["bus"][bus_idx, BASE_KV] / np.sqrt(3)
40	43
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)
	44	+	# Correct voltage of generator bus inside power station
	45	+	if np.any(~np.isnan(ppci["bus"][:, K_SG])):
	46	+	gen_bus_idx = bus_idx[~np.isnan(ppci["bus"][bus_idx, K_SG])]
	47	+	ppci["bus"][gen_bus_idx, IKSS1] *=\
	48	+	(ppci["bus"][gen_bus_idx, V_G] / ppci["bus"][gen_bus_idx, BASE_KV])
	49	+	ppci["bus"][gen_bus_idx, SKSS] *=\
	50	+	(ppci["bus"][gen_bus_idx, V_G] / ppci["bus"][gen_bus_idx, BASE_KV])
48	51
	52	+	_current_source_current(net, ppci)
	53	+
	54	+
	55	+	def _calc_ikss_1ph(net, ppci, ppci_0, bus_idx):
49	56		case = net._options["case"]
50		-	c = ppc["bus"][bus_idx, C_MIN] if case == "min" else ppc["bus"][bus_idx, C_MAX]
51		-	ppc["internal"]["baseI"] = ppc["bus"][:, BASE_KV] * np.sqrt(3) / ppc["baseMVA"]
52		-	ppc_0["internal"]["baseI"] = ppc_0["bus"][:, BASE_KV] * np.sqrt(3) / ppc_0["baseMVA"]
	57	+	c = ppci["bus"][bus_idx, C_MIN] if case == "min" else ppci["bus"][bus_idx, C_MAX]
	58	+	ppci["internal"]["baseI"] = ppci["bus"][:, BASE_KV] * np.sqrt(3) / ppci["baseMVA"]
	59	+	ppci_0["internal"]["baseI"] = ppci_0["bus"][:, BASE_KV] * np.sqrt(3) / ppci_0["baseMVA"]
	60	+
	61	+	z_equiv = abs((ppci["bus"][bus_idx, R_EQUIV] + ppci["bus"][bus_idx, X_EQUIV] * 1j) * 2 +
	62	+	(ppci_0["bus"][bus_idx, R_EQUIV] + ppci_0["bus"][bus_idx, X_EQUIV] * 1j))
	63	+
	64	+	# Only for test, should correspondant to PF result
	65	+	baseZ = ppci["bus"][bus_idx, BASE_KV] ** 2 / ppci["baseMVA"]
	66	+	ppci["bus"][bus_idx, R_EQUIV_OHM] = baseZ * ppci['bus'][bus_idx, R_EQUIV]
	67	+	ppci["bus"][bus_idx, X_EQUIV_OHM] = baseZ * ppci['bus'][bus_idx, X_EQUIV]
	68	+	ppci_0["bus"][bus_idx, R_EQUIV_OHM] = baseZ * ppci_0['bus'][bus_idx, R_EQUIV]
	69	+	ppci_0["bus"][bus_idx, X_EQUIV_OHM] = baseZ * ppci_0['bus'][bus_idx, X_EQUIV]
53	70
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))
	71	+	ppci_0["bus"][bus_idx, IKSS1] = c / z_equiv / ppci_0["bus"][bus_idx, BASE_KV] * np.sqrt(3) * ppci_0["baseMVA"]
	72	+	ppci["bus"][bus_idx, IKSS1] = c / z_equiv / ppci["bus"][bus_idx, BASE_KV] * np.sqrt(3) * ppci["baseMVA"]
56	73
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)
	74	+	_current_source_current(net, ppci)
60	75
61	76
62		-	def _current_source_current(net, ppc):
63		-	ppc["bus"][:, IKCV] = 0
64		-	ppc["bus"][:, IKSS2] = 0
	77	+	def _current_source_current(net, ppci):
	78	+	ppci["bus"][:, IKCV] = 0
	79	+	ppci["bus"][:, IKSS2] = 0
65	80		bus_lookup = net["_pd2ppc_lookups"]["bus"]
66	81		if not False in net.sgen.current_source.values:
67	82		sgen = net.sgen[net._is_elements["sgen"]]

71	86		return
72	87		if any(pd.isnull(sgen.sn_mva)):
73	88		raise ValueError("sn_mva needs to be specified for all sgens in net.sgen.sn_mva")
74		-	baseI = ppc["internal"]["baseI"]
	89	+	baseI = ppci["internal"]["baseI"]
75	90		sgen_buses = sgen.bus.values
76	91		sgen_buses_ppc = bus_lookup[sgen_buses]
77	92		if not "k" in sgen:
78	93		raise ValueError("Nominal to short-circuit current has to specified in net.sgen.k")
79	94		i_sgen_pu = sgen.sn_mva.values / net.sn_mva * sgen.k.values
80	95		buses, ikcv_pu, _ = _sum_by_group(sgen_buses_ppc, i_sgen_pu, i_sgen_pu)
81		-	ppc["bus"][buses, IKCV] = ikcv_pu
	96	+	ppci["bus"][buses, IKCV] = ikcv_pu
82	97		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)
	98	+	Zbus = ppci["internal"]["Zbus"]
	99	+	ppci["bus"][:, IKSS2] = abs(1 / np.diag(Zbus) * np.dot(Zbus, ppci["bus"][:, IKCV] * -1j) / baseI)
85	100		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)
89		-	ppc["bus"][buses, IKCV] /= baseI[buses]
	101	+	ybus_fact = ppci["internal"]["ybus_fact"]
	102	+	diagZ = _calc_zbus_diag(net, ppci)
	103	+	ppci["bus"][:, IKSS2] = abs(ybus_fact(ppci["bus"][:, IKCV] * -1j) / diagZ / baseI)
	104	+	ppci["bus"][buses, IKCV] /= baseI[buses]
90	105
91	106
92		-	def _calc_ip(net, ppc):
93		-	ip = np.sqrt(2) * (ppc["bus"][:, KAPPA] * ppc["bus"][:, IKSS1] + ppc["bus"][:, IKSS2])
94		-	ppc["bus"][:, IP] = ip
	107	+	def _calc_ip(net, ppci):
	108	+	ip = np.sqrt(2) * (ppci["bus"][:, KAPPA] * ppci["bus"][:, IKSS1] + ppci["bus"][:, IKSS2])
	109	+	ppci["bus"][:, IP] = ip
95	110
96	111
97		-	def _calc_ith(net, ppc):
	112	+	def _calc_ith(net, ppci):
98	113		tk_s = net["_options"]["tk_s"]
99		-	kappa = ppc["bus"][:, KAPPA]
	114	+	kappa = ppci["bus"][:, KAPPA]
100	115		f = 50
101	116		n = 1
102	117		m = (np.exp(4 * f * tk_s * np.log(kappa - 1)) - 1) / (2 * f * tk_s * np.log(kappa - 1))
103	118		m[np.where(kappa > 1.99)] = 0
104		-	ppc["bus"][:, M] = m
105		-	ith = (ppc["bus"][:, IKSS1] + ppc["bus"][:, IKSS2]) * np.sqrt(m + n)
106		-	ppc["bus"][:, ITH] = ith
107		-
108		-
109		-	def _calc_ib_generator(net, ppci):
110		-	Zbus = ppci["internal"]["Zbus"]
111		-	baseI = ppci["internal"]["baseI"]
112		-	tk_s = net._options['tk_s']
113		-	c = 1.1
114		-
115		-	z_equiv = ppci["bus"][:, R_EQUIV] + ppci["bus"][:, X_EQUIV] * 1j
116		-	I_ikss = c / z_equiv / ppci["bus"][:, BASE_KV] / np.sqrt(3) * ppci["baseMVA"]
117		-
118		-	# calculate voltage source branch current
119		-	# I_ikss = ppci["bus"][:, IKSS1]
120		-	V_ikss = (I_ikss * baseI) * Zbus
121		-
122		-	gen = net["gen"][net._is_elements["gen"]]
123		-	gen_vn_kv = gen.vn_kv.values
	119	+	ppci["bus"][:, M] = m
	120	+	ith = (ppci["bus"][:, IKSS1] + ppci["bus"][:, IKSS2]) * np.sqrt(m + n)
	121	+	ppci["bus"][:, ITH] = ith
124	122
125		-	gen_buses = ppci['gen'][:, GEN_BUS].astype(np.int64)
126		-	gen_mbase = ppci['gen'][:, MBASE]
127		-	gen_i_rg = gen_mbase / (np.sqrt(3) * gen_vn_kv)
128	123
129		-	gen_buses_ppc, gen_sn_mva, I_rG = _sum_by_group(gen_buses, gen_mbase, gen_i_rg)
	124	+	# TODO: Ib for generation close bus
	125	+	# def _calc_ib_generator(net, ppci):
	126	+	# # Zbus = ppci["internal"]["Zbus"]
	127	+	# # baseI = ppci["internal"]["baseI"]
	128	+	# tk_s = net._options['tk_s']
	129	+	# c = 1.1
130	130
131		-	# shunt admittance of generator buses and generator short circuit current
132		-	# YS = ppci["bus"][gen_buses_ppc, GS] + ppci["bus"][gen_buses_ppc, BS] * 1j
133		-	# I_kG = V_ikss.T[:, gen_buses_ppc] * YS / baseI[gen_buses_ppc]
	131	+	# z_equiv = ppci["bus"][:, R_EQUIV] + ppci["bus"][:, X_EQUIV] * 1j
	132	+	# I_ikss = c / z_equiv / ppci["bus"][:, BASE_KV] / np.sqrt(3) * ppci["baseMVA"]
134	133
135		-	xdss_pu = gen.xdss_pu.values
136		-	rdss_pu = gen.rdss_pu.values
137		-	cosphi = gen.cos_phi.values
138		-	X_dsss = xdss_pu * np.square(gen_vn_kv) / gen_mbase
139		-	R_dsss = rdss_pu * np.square(gen_vn_kv) / gen_mbase
	134	+	# # calculate voltage source branch current
	135	+	# # I_ikss = ppci["bus"][:, IKSS1]
	136	+	# # V_ikss = (I_ikss * baseI) * Zbus
140	137
141		-	K_G = ppci['bus'][gen_buses, BASE_KV] / gen_vn_kv * c / (1 + xdss_pu * np.sin(np.arccos(cosphi)))
142		-	Z_G = (R_dsss + 1j * X_dsss)
	138	+	# gen = net["gen"][net._is_elements["gen"]]
	139	+	# gen_vn_kv = gen.vn_kv.values
143	140
144		-	I_kG = c * ppci['bus'][gen_buses, BASE_KV] / np.sqrt(3) / (Z_G * K_G) * ppci["baseMVA"]
	141	+	# # Check difference ext_grid and gen
	142	+	# gen_buses = ppci['gen'][:, GEN_BUS].astype(np.int64)
	143	+	# gen_mbase = ppci['gen'][:, MBASE]
	144	+	# gen_i_rg = gen_mbase / (np.sqrt(3) * gen_vn_kv)
145	145
146		-	dV_G = 1j * X_dsss * K_G * I_kG
147		-	V_Is = c * ppci['bus'][gen_buses, BASE_KV] / np.sqrt(3)
	146	+	# gen_buses_ppc, gen_sn_mva, I_rG = _sum_by_group(gen_buses, gen_mbase, gen_i_rg)
148	147
149		-	# I_kG_contribution = I_kG.sum(axis=1)
150		-	# ratio_SG_ikss = I_kG_contribution / I_ikss
151		-	# close_to_SG = ratio_SG_ikss > 5e-2
	148	+	# # shunt admittance of generator buses and generator short circuit current
	149	+	# # YS = ppci["bus"][gen_buses_ppc, GS] + ppci["bus"][gen_buses_ppc, BS] * 1j
	150	+	# # I_kG = V_ikss.T[:, gen_buses_ppc] * YS / baseI[gen_buses_ppc]
152	151
153		-	close_to_SG = I_kG / I_rG > 2
	152	+	# xdss_pu = gen.xdss_pu.values
	153	+	# rdss_pu = gen.rdss_pu.values
	154	+	# cosphi = gen.cos_phi.values
	155	+	# X_dsss = xdss_pu * np.square(gen_vn_kv) / gen_mbase
	156	+	# R_dsss = rdss_pu * np.square(gen_vn_kv) / gen_mbase
154	157
155		-	if tk_s == 2e-2:
156		-	mu = 0.84 + 0.26 * np.exp(-0.26 * abs(I_kG) / I_rG)
157		-	elif tk_s == 5e-2:
158		-	mu = 0.71 + 0.51 * np.exp(-0.3 * abs(I_kG) / I_rG)
159		-	elif tk_s == 10e-2:
160		-	mu = 0.62 + 0.72 * np.exp(-0.32 * abs(I_kG) / I_rG)
161		-	elif tk_s >= 25e-2:
162		-	mu = 0.56 + 0.94 * np.exp(-0.38 * abs(I_kG) / I_rG)
163		-	else:
164		-	raise UserWarning('not implemented for other tk_s than 20ms, 50ms, 100ms and >=250ms')
	158	+	# K_G = ppci['bus'][gen_buses, BASE_KV] / gen_vn_kv * c / (1 + xdss_pu * np.sin(np.arccos(cosphi)))
	159	+	# Z_G = (R_dsss + 1j * X_dsss)
165	160
166		-	mu = np.clip(mu, 0, 1)
	161	+	# I_kG = c * ppci['bus'][gen_buses, BASE_KV] / np.sqrt(3) / (Z_G * K_G) * ppci["baseMVA"]
167	162
168		-	I_ikss_G = abs(I_ikss - np.sum((1 - mu) * I_kG, axis=1))
	163	+	# dV_G = 1j * X_dsss * K_G * I_kG
	164	+	# V_Is = c * ppci['bus'][gen_buses, BASE_KV] / np.sqrt(3)
169	165
170		-	# I_ikss_G = I_ikss - np.sum(abs(V_ikss.T[:, gen_buses_ppc]) * (1-mu) * I_kG, axis=1)
	166	+	# # I_kG_contribution = I_kG.sum(axis=1)
	167	+	# # ratio_SG_ikss = I_kG_contribution / I_ikss
	168	+	# # close_to_SG = ratio_SG_ikss > 5e-2
171	169
172		-	I_ikss_G = abs(I_ikss - np.sum(dV_G / V_Is * (1 - mu) * I_kG, axis=1))
	170	+	# close_to_SG = I_kG / I_rG > 2
173	171
174		-	return I_ikss_G
	172	+	# if tk_s == 2e-2:
	173	+	# mu = 0.84 + 0.26 * np.exp(-0.26 * abs(I_kG) / I_rG)
	174	+	# elif tk_s == 5e-2:
	175	+	# mu = 0.71 + 0.51 * np.exp(-0.3 * abs(I_kG) / I_rG)
	176	+	# elif tk_s == 10e-2:
	177	+	# mu = 0.62 + 0.72 * np.exp(-0.32 * abs(I_kG) / I_rG)
	178	+	# elif tk_s >= 25e-2:
	179	+	# mu = 0.56 + 0.94 * np.exp(-0.38 * abs(I_kG) / I_rG)
	180	+	# else:
	181	+	# raise UserWarning('not implemented for other tk_s than 20ms, 50ms, 100ms and >=250ms')
175	182
	183	+	# mu = np.clip(mu, 0, 1)
176	184
177		-	def _calc_single_bus_sc(net, ppc, bus):
178		-	# case = net._options["case"]
179		-	bus_idx = net._pd2ppc_lookups["bus"][bus]
180		-	n = ppc["bus"].shape[0]
181		-	Zbus = ppc["internal"]["Zbus"]
182		-	# Yf = ppc["internal"]["Yf"]
183		-	# Yt = ppc["internal"]["Yf"]
184		-	baseI = ppc["internal"]["baseI"]
185		-	# fb = np.real(ppc["branch"][:, 0]).astype(int)
186		-	# tb = np.real(ppc["branch"][:, 1]).astype(int)
187		-	# c = ppc["bus"][:, C_MIN] if case == "min" else ppc["bus"][:, C_MAX]
	185	+	# I_ikss_G = abs(I_ikss - np.sum((1 - mu) * I_kG, axis=1))
188	186
189		-	# calculate voltage source branch current
190		-	V_ikss = (ppc["bus"][:, IKSS1] * baseI) * Zbus
191		-	V = V_ikss[:, bus_idx]
192		-	# ikss_all_f = np.conj(Yf.dot(V_ikss))
193		-	# ikss_all_t = np.conj(Yt.dot(V_ikss))
194		-	current_sources = any(ppc["bus"][:, IKCV]) > 0
195		-	if current_sources:
196		-	current = np.tile(-ppc["bus"][:, IKCV], (n, 1))
197		-	np.fill_diagonal(current, current.diagonal() + ppc["bus"][:, IKSS2])
198		-	V_source = np.dot((current * baseI), Zbus).T
199		-	V = V + V_source[:, bus_idx]
200		-	# add current source branch current if there is one
201		-	# ppc["branch"][:, IKSS_F] = abs(ikss_all_f[:, bus_idx] / baseI[fb])
202		-	# ppc["branch"][:, IKSS_T] = abs(ikss_all_t[:, bus_idx] / baseI[tb])
203		-	calc_branch_results(net, ppc, V)
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)
212		-
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)
	187	+	# # I_ikss_G = I_ikss - np.sum(abs(V_ikss.T[:, gen_buses_ppc]) * (1-mu) * I_kG, axis=1)
252	188
	189	+	# I_ikss_G = abs(I_ikss - np.sum(dV_G / V_Is * (1 - mu) * I_kG, axis=1))
253	190
254		-	def calc_branch_results(net, ppci, V):
255		-	Ybus = ppci["internal"]["Ybus"]
256		-	Yf = ppci["internal"]["Yf"]
257		-	Yt = ppci["internal"]["Yt"]
258		-	baseMVA, bus, gen, branch, ref, _, _, _, _, _, ref_gens = _get_pf_variables_from_ppci(ppci)
259		-	bus, gen, branch = pfsoln_pypower(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, ref_gens)
260		-	ppci["bus"], ppci["gen"], ppci["branch"] = bus, gen, branch
	191	+	# return I_ikss_G
261	192
262	193
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]]
	194	+	def _calc_branch_currents(net, ppci, bus_idx):
272	195		n_sc_bus = np.shape(bus_idx)[0]
273	196
274	197		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	198		minmax = np.nanmin if case == "min" else np.nanmax
283	199
	200	+	Yf = ppci["internal"]["Yf"]
	201	+	Yt = ppci["internal"]["Yt"]
	202	+	baseI = ppci["internal"]["baseI"]
	203	+	n_bus = ppci["bus"].shape[0]
	204	+	fb = np.real(ppci["branch"][:, 0]).astype(int)
	205	+	tb = np.real(ppci["branch"][:, 1]).astype(int)
	206	+
284	207		# calculate voltage source branch current
285	208		if net["_options"]["inverse_y"]:
286		-	Zbus = ppc["internal"]["Zbus"]
287		-	V_ikss = (ppc["bus"][:, IKSS1] * baseI) * Zbus
	209	+	Zbus = ppci["internal"]["Zbus"]
	210	+	V_ikss = (ppci["bus"][:, IKSS1] * baseI) * Zbus
288	211		V_ikss = V_ikss[:, bus_idx]
289	212		else:
290		-	ybus_fact = ppc["internal"]["ybus_fact"]
	213	+	ybus_fact = ppci["internal"]["ybus_fact"]
291	214		V_ikss = np.zeros((n_bus, n_sc_bus), dtype=np.complex)
292	215		for ix, b in enumerate(bus_idx):
293	216		ikss = np.zeros(n_bus, dtype=np.complex)
294		-	ikss[b] = ppc["bus"][b, IKSS1] * baseI[b]
	217	+	ikss[b] = ppci["bus"][b, IKSS1] * baseI[b]
295	218		V_ikss[:, ix] = ybus_fact(ikss)
296	219
297	220		ikss1_all_f = np.conj(Yf.dot(V_ikss))

300	223		ikss1_all_t[abs(ikss1_all_t) < 1e-10] = 0.
301	224
302	225		# add current source branch current if there is one
303		-	current_sources = any(ppc["bus"][:, IKCV]) > 0
	226	+	current_sources = any(ppci["bus"][:, IKCV]) > 0
304	227		if current_sources:
305		-	current = np.tile(-ppc["bus"][:, IKCV], (n_sc_bus, 1))
	228	+	current = np.tile(-ppci["bus"][:, IKCV], (n_sc_bus, 1))
306	229		for ix, b in enumerate(bus_idx):
307		-	current[ix, b] += ppc["bus"][b, IKSS2]
	230	+	current[ix, b] += ppci["bus"][b, IKSS2]
308	231
309	232		# calculate voltage source branch current
310	233		if net["_options"]["inverse_y"]:
311		-	Zbus = ppc["internal"]["Zbus"]
	234	+	Zbus = ppci["internal"]["Zbus"]
312	235		V = np.dot((current * baseI), Zbus).T
313	236		else:
314		-	ybus_fact = ppc["internal"]["ybus_fact"]
	237	+	ybus_fact = ppci["internal"]["ybus_fact"]
315	238		V = np.zeros((n_bus, n_sc_bus), dtype=np.complex)
316	239		for ix, b in enumerate(bus_idx):
317	240		V[:, ix] = ybus_fact(current[ix, :] * baseI[b])
318	241
319		-	fb = np.real(ppc["branch"][:, 0]).astype(int)
320		-	tb = np.real(ppc["branch"][:, 1]).astype(int)
	242	+	fb = np.real(ppci["branch"][:, 0]).astype(int)
	243	+	tb = np.real(ppci["branch"][:, 1]).astype(int)
321	244		ikss2_all_f = np.conj(Yf.dot(V))
322	245		ikss2_all_t = np.conj(Yt.dot(V))
323	246

346	269		ip_all_t = np.sqrt(2) * ikss1_all_t * kappa[bus_idx]
347	270
348	271		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]
	272	+	ppci["internal"]["branch_ip_f"] = abs(ip_all_f) / baseI[fb, None]
	273	+	ppci["internal"]["branch_ip_t"] = abs(ip_all_t) / baseI[tb, None]
351	274		else:
352	275		ip_all_f[abs(ip_all_f) < 1e-10] = np.nan
353	276		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])
	277	+	ppci["branch"][:, IP_F] = np.nan_to_num(minmax(abs(ip_all_f), axis=1) / baseI[fb])
	278	+	ppci["branch"][:, IP_T] = np.nan_to_num(minmax(abs(ip_all_t), axis=1) / baseI[tb])
356	279
357	280		if net._options["ith"]:
358	281		n = 1
359		-	m = ppc["bus"][bus_idx, M]
	282	+	m = ppci["bus"][bus_idx, M]
360	283		ith_all_f = ikss_all_f * np.sqrt(m + n)
361	284		ith_all_t = ikss_all_t * np.sqrt(m + n)
362	285
363	286		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]
	287	+	ppci["internal"]["branch_ith_f"] = ith_all_f / baseI[fb, None]
	288	+	ppci["internal"]["branch_ith_t"] = ith_all_t / baseI[tb, None]
366	289		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])
	290	+	ppci["branch"][:, ITH_F] = np.nan_to_num(minmax(ith_all_f, axis=1) / baseI[fb])
	291	+	ppci["branch"][:, ITH_T] = np.nan_to_num(minmax(ith_all_t, axis=1) / baseI[fb])
	292	+
	293	+	# Update bus index for branch results
	294	+	if net._options["return_all_currents"]:
	295	+	ppci["internal"]["br_res_ks_ppci_bus"] = bus_idx

17	15		IK_T = start + 7
18	16		IB_F = start + 8
19	17		IB_T = start + 9
	18	+	K_T = start + 10
	19	+	K_ST = start + 11
20	20
21		-	branch_cols_sc = 10

506	506		ppc["bus"][b, GS] = vp
507	507		ppc["bus"][b, BS] = -vq
508	508
509		-
510		-	def _add_gen_impedances_ppc(net, ppc):
511		-	_add_ext_grid_sc_impedance(net, ppc)
512		-	_add_gen_sc_impedance(net, ppc)
513		-
514		-
	509	+	# Short circuit relevant routines
515	510		def _add_ext_grid_sc_impedance(net, ppc):
516	511		from pandapower.shortcircuit.idx_bus import C_MAX, C_MIN
517	512		mode = net._options["mode"]

553	548		ppc["bus"][buses, BS] = bs * ppc['baseMVA']
554	549		return gs * ppc['baseMVA'], bs * ppc['baseMVA']
555	550
556		-
557		-	def _add_gen_sc_impedance(net, ppc):
558		-	from pandapower.shortcircuit.idx_bus import C_MAX
559		-	gen = net["gen"][net._is_elements["gen"]]
560		-	if len(gen) == 0:
561		-	return
562		-	gen_buses = gen.bus.values
563		-	bus_lookup = net["_pd2ppc_lookups"]["bus"]
564		-	gen_buses_ppc = bus_lookup[gen_buses]
565		-	vn_gen = gen.vn_kv.values
566		-	sn_gen = gen.sn_mva.values
567		-
568		-	rdss_pu = gen.rdss_pu.values
569		-	xdss_pu = gen.xdss_pu.values
570		-	gens_without_r = np.isnan(rdss_pu)
571		-	if gens_without_r.any():
572		-	# use the estimations from the IEC standard for generators without defined rdss_pu
573		-	lv_gens = (vn_gen <= 1.) & gens_without_r
574		-	hv_gens = (vn_gen > 1.) & gens_without_r
575		-	large_hv_gens = (sn_gen >= 100) & hv_gens
576		-	small_hv_gens = (sn_gen < 100) & hv_gens
577		-	rdss_pu[lv_gens] = 0.15 * xdss_pu[lv_gens]
578		-	rdss_pu[large_hv_gens] = 0.05 * xdss_pu[large_hv_gens]
579		-	rdss_pu[small_hv_gens] = 0.07 * xdss_pu[small_hv_gens]
580		-
581		-	vn_net = ppc["bus"][gen_buses_ppc, BASE_KV]
582		-	cmax = ppc["bus"][gen_buses_ppc, C_MAX]
583		-	phi_gen = np.arccos(gen.cos_phi.values)
584		-	kg = vn_gen / vn_net * cmax / (1 + xdss_pu * np.sin(phi_gen))
585		-
586		-	z_gen = (rdss_pu + xdss_pu * 1j) * kg / sn_gen
587		-	y_gen = 1 / z_gen
588		-
589		-	buses, gs, bs = _sum_by_group(gen_buses_ppc, y_gen.real, y_gen.imag)
590		-	ppc["bus"][buses, GS] = gs
591		-	ppc["bus"][buses, BS] = bs
592		-
593		-
594	551		def _add_motor_impedances_ppc(net, ppc):
595	552		if net._options["case"] == "min":
596	553		return

20	20		IB = start + 10
21	21		ITH = start + 11
22	22		IK = start + 12
	23	+	R_EQUIV_OHM = start + 13
	24	+	X_EQUIV_OHM = start + 14
	25	+	K_G = start + 15
	26	+	K_SG = start + 16
	27	+	V_G = start + 17
	28	+	PS_TRAFO_IX = start + 18
	29	+	GS_P = start + 19
	30	+	BS_P = start + 20
	31	+	SKSS = start + 21
23	32
24		-	bus_cols_sc = 13

18	18		from pandapower.pd2ppc import _pd2ppc
19	19		from pandapower.pd2ppc_zero import _pd2ppc_zero
20	20		from pandapower.results import _copy_results_ppci_to_ppc
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	+
	22	+	from pandapower.shortcircuit.currents import _calc_ikss,\
	23	+	_calc_ikss_1ph, _calc_ip, _calc_ith, _calc_branch_currents
24	24		from pandapower.shortcircuit.impedance import _calc_zbus, _calc_ybus, _calc_rx
	25	+	from pandapower.shortcircuit.ppc_conversion import _init_ppc, _create_k_updated_ppci, _get_is_ppci_bus
25	26		from pandapower.shortcircuit.kappa import _add_kappa_to_ppc
26		-	from pandapower.shortcircuit.results import _extract_results, _extract_single_results
	27	+	from pandapower.shortcircuit.results import _extract_results, _copy_result_to_ppci_orig
27	28		from pandapower.results import init_results
28	29
29	30
30		-	def calc_sc(net, fault="3ph", case='max', lv_tol_percent=10, topology="auto", ip=False,
	31	+	def calc_sc(net, bus=None,
	32	+	fault="3ph", case='max', lv_tol_percent=10, topology="auto", ip=False,
31	33		ith=False, tk_s=1., kappa_method="C", r_fault_ohm=0., x_fault_ohm=0.,
32	34		branch_results=False, check_connectivity=True, return_all_currents=False,
33		-	bus=None, inverse_y=True):
	35	+	inverse_y=True):
34	36		"""
35	37		Calculates minimal or maximal symmetrical short-circuit currents.
36	38		The calculation is based on the method of the equivalent voltage source

46	48		INPUT:
47	49		net (pandapowerNet) pandapower Network
48	50
	51	+	bus (int, list, np.array, None) defines if short-circuit calculations should only be calculated for defined bus
	52	+
49	53		*fault (str, 3ph) type of fault
50	54
51	55		- "3ph" for three-phase

68	72
69	73		ip (bool, False) if True, calculate aperiodic short-circuit current
70	74
71		-	Ith (bool, False) if True, calculate equivalent thermical short-circuit current Ith
	75	+	ith (bool, False) if True, calculate equivalent thermical short-circuit current Ith
72	76
73	77		topology (str, "auto") define option for meshing (only relevant for ip and ith)
74	78

89	93		return_all_currents (bool, False) applies only if branch_results=True, if True short-circuit currents for
90	94		each (branch, bus) tuple is returned otherwise only the max/min is returned
91	95
92		-	bus (int, list, np.array, None) defines if short-circuit calculations should only be calculated for defined bus
93		-
94	96		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	97
96	98

106	108		"Only 3ph, 2ph and 1ph short-circuit currents implemented")
107	109
108	110		if len(net.gen) and (ip or ith):
109		-	logger.warning("aperiodic and thermal short-circuit currents are only implemented for "
	111	+	logger.warning("aperiodic, thermal short-circuit currents are only implemented for "
110	112		"faults far from generators!")
111	113
112	114		if case not in ['max', 'min']:
113	115		raise ValueError('case can only be "min" or "max" for minimal or maximal short "\
114	116		"circuit current')
	117	+
115	118		if topology not in ["meshed", "radial", "auto"]:
116	119		raise ValueError(
117	120		'specify network structure as "meshed", "radial" or "auto"')

120	123		logger.warning("Branch results are in beta mode and might not always be reliable, "
121	124		"especially for transformers")
122	125
123		-	# Convert bus to numpy array for better performance
124		-	if isinstance(bus, int):
	126	+	# Convert bus to numpy array
	127	+	if bus is None:
	128	+	bus = net.bus.index.values
	129	+	elif isinstance(bus, int):
125	130		bus = np.array([bus])
126		-	elif isinstance(bus, list):
	131	+	elif not isinstance(bus, np.ndarray):
127	132		bus = np.array(bus)
128	133
129	134		kappa = ith or ip
130	135		net["_options"] = {}
131	136		_add_ppc_options(net, calculate_voltage_angles=False, trafo_model="pi",
132		-	check_connectivity=check_connectivity, mode="sc", switch_rx_ratio=2,
133		-	init_vm_pu="flat", init_va_degree="flat", enforce_q_lims=False,
134		-	recycle=None)
	137	+	check_connectivity=check_connectivity, mode="sc", switch_rx_ratio=2,
	138	+	init_vm_pu="flat", init_va_degree="flat", enforce_q_lims=False,
	139	+	recycle=None)
135	140		_add_sc_options(net, fault=fault, case=case, lv_tol_percent=lv_tol_percent, tk_s=tk_s,
136	141		topology=topology, r_fault_ohm=r_fault_ohm, kappa_method=kappa_method,
137	142		x_fault_ohm=x_fault_ohm, kappa=kappa, ip=ip, ith=ith,
138	143		branch_results=branch_results, return_all_currents=return_all_currents,
139	144		inverse_y=inverse_y)
140	145		init_results(net, "sc")
	146	+
141	147		if fault in ("2ph", "3ph"):
142	148		_calc_sc(net, bus)
143	149		elif fault == "1ph":

146	152		raise ValueError("Invalid fault %s" % fault)
147	153
148	154
149		-	def calc_single_sc(net, bus, fault="3ph", case='max', lv_tol_percent=10,
150		-	check_connectivity=True, inverse_y=True):
151		-	"""
152		-	Calculates minimal or maximal symmetrical short-circuit currents.
153		-	The calculation is based on the method of the equivalent voltage source
154		-	according to DIN/IEC EN 60909.
155		-	The initial short-circuit alternating current ikss is the basis of the short-circuit
156		-	calculation and is therefore always calculated.
157		-	Other short-circuit currents can be calculated from ikss with the conversion factors defined
158		-	in DIN/IEC EN 60909.
159		-
160		-	The output is stored in the net.res_bus_sc table as a short_circuit current
161		-	for each bus.
162		-
163		-	INPUT:
164		-	net (pandapowerNet) pandapower Network
165		-
166		-	*fault (str, 3ph) type of fault
167		-
168		-	- "3ph" for three-phase
169		-
170		-	- "2ph" for two-phase short-circuits
171		-
172		-	case (str, "max")
173		-
174		-	- "max" for maximal current calculation
175		-
176		-	- "min" for minimal current calculation
177		-
178		-	lv_tol_percent (int, 10) voltage tolerance in low voltage grids
179		-
180		-	- 6 for 6% voltage tolerance
181		-
182		-	- 10 for 10% voltage olerance
183		-
184		-	r_fault_ohm (float, 0) fault resistance in Ohm
185		-
186		-	x_fault_ohm (float, 0) fault reactance in Ohm
187		-
188		-	OUTPUT:
189		-
190		-	EXAMPLE:
191		-	calc_sc(net)
	155	+	def _calc_current(net, ppci_orig, bus):
	156	+	# Select required ppci bus
	157	+	ppci_bus = _get_is_ppci_bus(net, bus)
192	158
193		-	print(net.res_bus_sc)
194		-	"""
195		-	if fault not in ["3ph", "2ph"]:
196		-	raise NotImplementedError("Only 3ph and 2ph short-circuit currents implemented")
197		-
198		-	if case not in ['max', 'min']:
199		-	raise ValueError('case can only be "min" or "max" for minimal or maximal short "\
200		-	"circuit current')
201		-	net["_options"] = {}
202		-	_add_ppc_options(net, calculate_voltage_angles=False, trafo_model="pi",
203		-	check_connectivity=check_connectivity, mode="sc", switch_rx_ratio=2,
204		-	init_vm_pu="flat", init_va_degree="flat", enforce_q_lims=False,
205		-	recycle=None)
206		-	_add_sc_options(net, fault=fault, case=case, lv_tol_percent=lv_tol_percent, tk_s=1.,
207		-	topology="auto", r_fault_ohm=0., kappa_method="C",
208		-	x_fault_ohm=0., kappa=False, ip=False, ith=False,
209		-	branch_results=True, return_all_currents=False,
210		-	inverse_y=inverse_y)
211		-	init_results(net, "sc")
212		-	if fault in ("2ph", "3ph"):
213		-	_calc_sc_single(net, bus)
214		-	elif fault == "1ph":
215		-	raise NotImplementedError("1ph short-circuits are not yet implemented")
216		-	else:
217		-	raise ValueError("Invalid fault %s" % fault)
	159	+	# update ppci
	160	+	non_ps_gen_ppci_bus, non_ps_gen_ppci, ps_gen_bus_ppci_dict =\
	161	+	_create_k_updated_ppci(net, ppci_orig, ppci_bus=ppci_bus)
218	162
	163	+	# For each ps_gen_bus one unique ppci is required
	164	+	ps_gen_ppci_bus = list(ps_gen_bus_ppci_dict.keys())
219	165
220		-	def _calc_sc_single(net, bus):
221		-	_add_auxiliary_elements(net)
222		-	ppc, ppci = _pd2ppc(net)
223		-	_calc_ybus(ppci)
	166	+	for calc_bus in ps_gen_ppci_bus+[non_ps_gen_ppci_bus]:
	167	+	if isinstance(calc_bus, np.ndarray):
	168	+	# Use ppci for general bus
	169	+	this_ppci, this_ppci_bus = non_ps_gen_ppci, calc_bus
	170	+	else:
	171	+	# Use specific ps_gen_bus ppci
	172	+	this_ppci, this_ppci_bus = ps_gen_bus_ppci_dict[calc_bus], np.array([calc_bus])
224	173
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"])
	174	+	_calc_ybus(this_ppci)
	175	+	if net["_options"]["inverse_y"]:
	176	+	_calc_zbus(net, this_ppci)
	177	+	else:
	178	+	# Factorization Ybus once
	179	+	this_ppci["internal"]["ybus_fact"] = factorized(this_ppci["internal"]["Ybus"])
233	180
234		-	_calc_rx(net, ppci, bus)
235		-	_calc_ikss(net, ppci, bus)
236		-	_calc_single_bus_sc_no_y_inv(net, ppci, bus)
	181	+	_calc_rx(net, this_ppci, this_ppci_bus)
	182	+	_calc_ikss(net, this_ppci, this_ppci_bus)
	183	+	_add_kappa_to_ppc(net, this_ppci)
	184	+	if net["_options"]["ip"]:
	185	+	_calc_ip(net, this_ppci)
	186	+	if net["_options"]["ith"]:
	187	+	_calc_ith(net, this_ppci)
237	188
238		-	# Delete factorization object
239		-	ppci["internal"].pop("ybus_fact")
	189	+	if net._options["branch_results"]:
	190	+	_calc_branch_currents(net, this_ppci, this_ppci_bus)
240	191
241		-	ppc = _copy_results_ppci_to_ppc(ppci, ppc, "sc")
242		-	_extract_single_results(net, ppc)
243		-	_clean_up(net)
	192	+	_copy_result_to_ppci_orig(ppci_orig, this_ppci, this_ppci_bus,
	193	+	calc_options=net._options)
244	194
245	195
246	196		def _calc_sc(net, bus):
247		-	_add_auxiliary_elements(net)
248		-	ppc, ppci = _pd2ppc(net)
249		-	_calc_ybus(ppci)
	197	+	ppc, ppci = _init_ppc(net)
250	198
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		-
264		-	if net["_options"]["ip"]:
265		-	_calc_ip(net, ppci)
266		-	if net["_options"]["ith"]:
267		-	_calc_ith(net, ppci)
268		-
269		-	if net._options["branch_results"]:
270		-	_calc_branch_currents(net, ppci, bus)
	199	+	_calc_current(net, ppci, bus)
271	200
272	201		ppc = _copy_results_ppci_to_ppc(ppci, ppc, "sc")
273	202		_extract_results(net, ppc, ppc_0=None, bus=bus)

299	228		ppci["internal"]["ybus_fact"] = factorized(ppci["internal"]["Ybus"])
300	229		ppci_0["internal"]["ybus_fact"] = factorized(ppci_0["internal"]["Ybus"])
301	230
302		-	_calc_rx(net, ppci, bus=bus)
	231	+	ppci_bus = _get_is_ppci_bus(net, bus)
	232	+	_calc_rx(net, ppci, ppci_bus)
303	233		_add_kappa_to_ppc(net, ppci)
304	234
305		-	_calc_rx(net, ppci_0, bus=bus)
306		-	_calc_ikss_1ph(net, ppci, ppci_0, bus=bus)
	235	+	_calc_rx(net, ppci_0, ppci_bus)
	236	+	_calc_ikss_1ph(net, ppci, ppci_0, ppci_bus)
307	237
308	238		if net._options["branch_results"]:
309		-	_calc_branch_currents(net, ppci, bus=bus)
	239	+	_calc_branch_currents(net, ppci, ppci_bus)
	240	+
310	241		ppc_0 = _copy_results_ppci_to_ppc(ppci_0, ppc_0, "sc")
311	242		ppc = _copy_results_ppci_to_ppc(ppci, ppc, "sc")
312	243		_extract_results(net, ppc, ppc_0, bus=bus)

35	71		_get_trafo_results(net, ppc)
36	72		_get_trafo3w_results(net, ppc)
37	73
38		-	def _extract_single_results(net, ppc):
39		-	for element in ["line", "trafo"]:
40		-	net["res_%s_sc"%element] = pd.DataFrame(np.nan, index=net[element].index,
41		-	columns=net["_empty_res_%s"%element].columns,
42		-	dtype='float')
43		-	_get_single_bus_results(net, ppc)
44		-	net["_options"]["ac"] = True
45		-	net["_options"]["trafo_loading"] = "current"
46		-	bus_lookup_aranged = _get_aranged_lookup(net)
47		-	bus_pq = np.zeros(shape=(len(net["bus"].index), 2), dtype=np.float)
48		-	_get_branch_results(net, ppc, bus_lookup_aranged, bus_pq, suffix="_sc")
49		-
50		-
51		-	def _get_single_bus_results(net, ppc):
52		-	_set_buses_out_of_service(ppc)
53		-	bus_idx = _get_bus_idx(net)
54		-	case = net._options["case"]
55		-	c = ppc["bus"][bus_idx, C_MIN] if case == "min" else ppc["bus"][bus_idx, C_MAX]
56		-	net["res_bus"]["vm_pu"] = np.nan
57		-	net["res_bus_sc"]["vm_pu"] = c - ppc["bus"][bus_idx, VM]
58		-	net["res_bus_sc"]["va_degree"] = ppc["bus"][bus_idx, VA]
59		-
60	74
61	75		def _get_bus_results(net, ppc, ppc_0, bus):
62		-	if bus is None:
63		-	bus = slice(None)
64		-
65	76		bus_lookup = net._pd2ppc_lookups["bus"]
66	77		ppc_index = bus_lookup[net.bus.index]
67	78
68	79		if net["_options"]["fault"] == "1ph":
69	80		net.res_bus_sc["ikss_ka"] = ppc_0["bus"][ppc_index, IKSS1] + ppc["bus"][ppc_index, IKSS2]
	81	+	net.res_bus_sc["rk0_ohm"] = ppc_0["bus"][ppc_index, R_EQUIV_OHM]
	82	+	net.res_bus_sc["xk0_ohm"] = ppc_0["bus"][ppc_index, X_EQUIV_OHM]
70	83		else:
71	84		net.res_bus_sc["ikss_ka"] = ppc["bus"][ppc_index, IKSS1] + ppc["bus"][ppc_index, IKSS2]
	85	+	net.res_bus_sc["skss_mw"] = ppc["bus"][ppc_index, SKSS]
72	86		if net._options["ip"]:
73	87		net.res_bus_sc["ip_ka"] = ppc["bus"][ppc_index, IP]
74	88		if net._options["ith"]:
75	89		net.res_bus_sc["ith_ka"] = ppc["bus"][ppc_index, ITH]
76	90
	91	+	# Export also equivalent rk, xk on the calculated bus
	92	+	net.res_bus_sc["rk_ohm"] = ppc["bus"][ppc_index, R_EQUIV_OHM]
	93	+	net.res_bus_sc["xk_ohm"] = ppc["bus"][ppc_index, X_EQUIV_OHM]
	94	+
77	95		net.res_bus_sc = net.res_bus_sc.loc[bus, :]
78	96
79	97

103	121		f, t = branch_lookup["line"]
104	122		minmax = np.maximum if case == "max" else np.minimum
105	123
106		-	net.res_line_sc["ikss_ka"] = minmax(ppc["internal"]["branch_ikss_f"][f:t, ppc_index].real.reshape(-1, 1),
107		-	ppc["internal"]["branch_ikss_t"][f:t, ppc_index].real.reshape(-1, 1))
	124	+	net.res_line_sc["ikss_ka"] = minmax(ppc["internal"]["branch_ikss_f"].iloc[f:t,:].loc[:, ppc_index].values.real.reshape(-1, 1),
	125	+	ppc["internal"]["branch_ikss_t"].iloc[f:t,:].loc[:, ppc_index].values.real.reshape(-1, 1))
108	126		if net._options["ip"]:
109		-	net.res_line_sc["ip_ka"] = minmax(ppc["internal"]["branch_ip_f"][f:t, ppc_index].real.reshape(-1, 1),
110		-	ppc["internal"]["branch_ip_t"][f:t, ppc_index].real.reshape(-1, 1))
	127	+	net.res_line_sc["ip_ka"] = minmax(ppc["internal"]["branch_ip_f"].iloc[f:t,:].loc[:, ppc_index].values.real.reshape(-1, 1),
	128	+	ppc["internal"]["branch_ip_t"].iloc[f:t,:].loc[:, ppc_index].values.real.reshape(-1, 1))
111	129		if net._options["ith"]:
112		-	net.res_line_sc["ith_ka"] = minmax(ppc["internal"]["branch_ith_f"][f:t, ppc_index].real.reshape(-1, 1),
113		-	ppc["internal"]["branch_ith_t"][f:t, ppc_index].real.reshape(-1, 1))
	130	+	net.res_line_sc["ith_ka"] = minmax(ppc["internal"]["branch_ith_f"].iloc[f:t,:].loc[:, ppc_index].values.real.reshape(-1, 1),
	131	+	ppc["internal"]["branch_ith_t"].iloc[f:t,:].loc[:, ppc_index].values.real.reshape(-1, 1))
114	132
115	133		def _get_trafo_results(net, ppc):
116	134		branch_lookup = net._pd2ppc_lookups["branch"]

129	147
130	148		if "trafo" in branch_lookup:
131	149		f, t = branch_lookup["trafo"]
132		-	net.res_trafo_sc["ikss_hv_ka"] = ppc["internal"]["branch_ikss_f"][f:t, ppc_index].real.reshape(-1, 1)
133		-	net.res_trafo_sc["ikss_lv_ka"] = ppc["internal"]["branch_ikss_t"][f:t, ppc_index].real.reshape(-1, 1)
	150	+	net.res_trafo_sc["ikss_hv_ka"] = ppc["internal"]["branch_ikss_f"].iloc[f:t,:].loc[:, ppc_index].values.real.reshape(-1, 1)
	151	+	net.res_trafo_sc["ikss_lv_ka"] = ppc["internal"]["branch_ikss_t"].iloc[f:t,:].loc[:, ppc_index].values.real.reshape(-1, 1)
134	152
135	153
136	154		def _get_trafo3w_results(net, ppc):

157	175		hv = int(f + (t - f) / 3)
158	176		mv = int(f + 2 * (t - f) / 3)
159	177		lv = t
160		-	net.res_trafo3w_sc["ikss_hv_ka"] = ppc["internal"]["branch_ikss_f"][f:hv, ppc_index].real.reshape(-1, 1)
161		-	net.res_trafo3w_sc["ikss_mv_ka"] = ppc["internal"]["branch_ikss_t"][hv:mv, ppc_index].real.reshape(-1, 1)
162		-	net.res_trafo3w_sc["ikss_lv_ka"] = ppc["internal"]["branch_ikss_t"][mv:lv, ppc_index].real.reshape(-1, 1)