Current Source Elements

[kostek184-coder]

How to model a current source such as IBR?

[nitbharambe]

Hi @kostek184-coder , we can model it as a sym_gen or sym_load in a steady state power flow or state estimation calculation. You can also specify it to be of constant power/current/impedance(ZIP) if required.

IBR related modelling depend on context of the calculation and features you would like to use. Let me know if you are looking for something different than a simple load or generator or a study other than steady state powerflow.

[kostek184-coder]

I would like to use them to short circuit calculation. This model could be the same as in panda power with definition of current angle of the current injection

https://pandapower.readthedocs.io/en/latest/shortcircuit/current_source.html

[nitbharambe]

@kostek184-coder Ok the context is clear now! So including current contributions from DER for short circuit calculations. I see there wouldn't be a direct way to simulate this in PGM at the moment.
However, this seems like a useful thing to support / enable via PGM. I shall get back to you soon on what can be done.

[nitbharambe]

Looking at the documentation of pandapower, in short, there can be a way around substituting a DER with a source component from PGM. I was planning to demonstrate this using a script. But I have limited access currently, so it might take up more time.

[nitbharambe]

Hi @kostek184-coder , there is not easy and direct way to do so via PGM at the moment. Feel free to create a feature request for such contributions.

The workaround way involving current injection via source component is possible but kind of complicated to implement. Hence do so if there is urgency:

Source as current injection demo

It involves first calculating short circuit current. Then we conduct a power flow generators modelled as a source with sk = sn_kva * k / scaling_factor and u_ref = 1.0 * scaling_factor where scaling_factor is any large number. By virtue of superposition theorem, we can then add all currents of short circuit and load flow to find resultant fault current.

TLDR: Attaching code of the "hack" below just for information but I would recommend avoiding it.

# /// script
# requires-python = ">=3.11"
# dependencies = [
#   "power-grid-model>=1.9",
#   "pandapower>=2.13",
#   "numpy",
# ]
# ///
"""Three-phase short-circuit comparison: power-grid-model vs pandapower

Topology
--------
  PGM:  source ── node0 ── line0 ── node1 [sym_gen] ── line1 ── node2 [fault]
  PP:   ext_grid ─ bus0 ── line0 ── bus1  [sgen]    ── line1 ── bus2  (calc_sc)

IEC 60909 source-impedance convention
--------------------------------------
  pandapower:  Z_src = c · Un² / Sk          (IEC 60909 eq. 18)
  PGM native:  Z_src = Un² / sk  (c applied to driving voltage E = c·u_ref·Un only)

  To align: pass sk_pgm = Sk / c_max so both tools use the same Thevenin impedance.
"""

import numpy as np
import pandapower as pp
import pandapower.shortcircuit as sc
from power_grid_model import CalculationMethod, PowerGridModel, initialize_array
from power_grid_model.enum import FaultPhase, FaultType, ShortCircuitVoltageScaling

# ── Shared grid parameters ─────────────────────────────────────────────────────
U_NOM_KV = 10.0
U_NOM_V  = U_NOM_KV * 1e3
FREQ_HZ  = 50.0

C_MAX    = 1.1     # IEC 60909 voltage factor, Un > 1 kV
SK_MVA   = 100.0   # Source short-circuit power [MVA]
RX_SRC   = 0.1     # Source R/X ratio

# Line 0: node0 → node1
R0_OHM = 0.1
X0_OHM = 0.1

# Line 1: node1 → node2
R1_OHM = 0.1
X1_OHM = 0.1

# Generator at node1/bus1  (topology only — no SC current-source contribution)
P_GEN_MW = 1.0
K_GEN = 0.7


# ── PGM ────────────────────────────────────────────────────────────────────────

def run_pgm() -> dict:
    node = initialize_array("input", "node", 3)
    node["id"]      = [0, 1, 2]
    node["u_rated"] = U_NOM_V

    source = initialize_array("input", "source", 1)
    source["id"]       = [10]
    source["node"]     = [0]
    source["status"]   = [1]
    source["u_ref"]    = [1.0]
    source["sk"]       = [SK_MVA * 1e6 / C_MAX]
    source["rx_ratio"] = [RX_SRC]

    line = initialize_array("input", "line", 2)
    line["id"]          = [20, 21]
    line["from_node"]   = [0,      1]
    line["to_node"]     = [1,      2]
    line["from_status"] = [1,      1]
    line["to_status"]   = [1,      1]
    line["r1"]          = [R0_OHM, R1_OHM]
    line["x1"]          = [X0_OHM, X1_OHM]
    line["c1"]          = [0.0,    0.0]
    line["tan1"]        = [0.0,    0.0]

    fault = initialize_array("input", "fault", 1)
    fault["id"]           = [40]
    fault["status"]       = [1]
    fault["fault_object"] = [2]              # node id of faulted bus
    fault["fault_type"]   = [FaultType.three_phase]
    fault["fault_phase"]  = [FaultPhase.abc]
    fault["r_f"]          = [0.0]
    fault["x_f"]          = [0.0]

    input_data = {"node": node, "source": source, "line": line, "sym_gen": sym_gen, "fault": fault}
    model = PowerGridModel(input_data)

    sc_output = model.calculate_short_circuit(
        calculation_method=CalculationMethod.iec60909,
        short_circuit_voltage_scaling=ShortCircuitVoltageScaling.maximum,
    )

    # Some calculation
    Z_src_abs = C_MAX * U_NOM_V**2 / (SK_MVA * 1e6)       # |Z_src| [Ω]
    X_src = Z_src_abs / np.sqrt(1.0 + RX_SRC**2)
    R_src = RX_SRC * X_src
    Z_src = R_src + 1j * X_src
    Y_src = 1.0 / Z_src

    SK_GEN = K_GEN * P_GEN_MW * 1e6
    # We make internal impedance really small while adusting the voltage reference to keep the current the same.
    SCALING_FACTOR = 1.0e3 

    # Souce representing sym_gen current contribution to fault
    source_2 = initialize_array("input", "source", 1)
    source_2["id"]       = [10]
    source_2["node"]     = [1]
    source_2["status"]   = [1]
    source_2["u_ref"]    = [1.1 * SCALING_FACTOR] 
    source_2["u_ref_angle"] = sc_output["node"]["u_angle"][1, 0]
    source_2["sk"]       = [SK_GEN / C_MAX / SCALING_FACTOR]
    source_2["rx_ratio"] = [RX_SRC]

    # Shunt representing source internal impedance at node0 and a shunt representing fault at node2
    shunt_2 = initialize_array("input", "shunt", 2)
    shunt_2["id"]     = [50, 51]
    shunt_2["node"]   = [0, 2]
    shunt_2["status"] = [1, 1]
    shunt_2["g1"] = [Y_src.real, 1e10]
    shunt_2["b1"] = [Y_src.imag, 1e10]
    
    input_data_2 = {"node": node, "source": source_2, "line": line, "shunt": shunt_2}

    model_2 = PowerGridModel(input_data_2)
    lf_output = model_2.calculate_power_flow(calculation_method=CalculationMethod.newton_raphson)

    sc_output["fault"]["i_f"] += lf_output["line"]["i_to"][1]  # add sym_gen contribution to fault current

    return sc_output

# ── pandapower ─────────────────────────────────────────────────────────────────

def run_pandapower() -> tuple:
    net = pp.create_empty_network(f_hz=FREQ_HZ, sn_mva=100.0)

    b0 = pp.create_bus(net, vn_kv=U_NOM_KV, name="bus0")
    b1 = pp.create_bus(net, vn_kv=U_NOM_KV, name="bus1")
    b2 = pp.create_bus(net, vn_kv=U_NOM_KV, name="bus2")

    pp.create_ext_grid(
        net, bus=b0, vm_pu=1.0,
        s_sc_max_mva=SK_MVA,  # pandapower applies Z = c·Un²/Sk internally
        rx_max=RX_SRC,
    )

    # length_km=1 so r_ohm_per_km == total R in ohms
    pp.create_line_from_parameters(
        net, from_bus=b0, to_bus=b1, length_km=1.0,
        r_ohm_per_km=R0_OHM, x_ohm_per_km=X0_OHM,
        c_nf_per_km=0.0, max_i_ka=99.0,
    )
    pp.create_line_from_parameters(
        net, from_bus=b1, to_bus=b2, length_km=1.0,
        r_ohm_per_km=R1_OHM, x_ohm_per_km=X1_OHM,
        c_nf_per_km=0.0, max_i_ka=99.0,
    )

    # sgen: current_source=False → topology element only, no SC current contribution
    pp.create_sgen(net, bus=b1, p_mw=P_GEN_MW, q_mvar=0.0, sn_mva=P_GEN_MW,
                   current_source=True, name="sgen", k=K_GEN)

    sc.calc_sc(net, fault="3ph", case="max", bus=b2)

    return net, b2

# ── Analytical (IEC 60909 equivalent-voltage-source method) ───────────────────

def analytical_ikss() -> float:
    """
    IEC 60909 §4.3.1 — three-phase bolted fault at node2.

    Step 1: source Thevenin impedance (IEC 60909 eq. 18)
        Z_src = c · Un² / Sk
    Step 2: series sum of all impedances from source to fault
        Z_total = Z_src + Z_line0 + Z_line1
    Step 3: initial symmetrical short-circuit current (IEC 60909 eq. 29)
        Ikss = c · Un / (√3 · |Z_total|)   [A]

    The sym_gen at node1 is a PQ load-flow element with no SC current-source
    model, so it does not appear in the Thevenin equivalent.
    """
    Z_src_abs = C_MAX * U_NOM_V**2 / (SK_MVA * 1e6)       # |Z_src| [Ω]
    X_src = Z_src_abs / np.sqrt(1.0 + RX_SRC**2)
    R_src = RX_SRC * X_src
    Z_src = R_src + 1j * X_src

    Z_line0 = R0_OHM + 1j * X0_OHM
    Z_line1 = R1_OHM + 1j * X1_OHM
    Z_tot = Z_src + Z_line0 + Z_line1
    Z_tot_abs = np.abs(Z_tot)

    i_fault_source = C_MAX * U_NOM_V / (np.sqrt(3.0) * Z_tot_abs)           # [A]

    Z_gen_to_source_abs = np.abs(Z_src + Z_line0)
    fault_contribution_ratio = Z_gen_to_source_abs / Z_tot_abs
    i_fault_gen = + fault_contribution_ratio * K_GEN * P_GEN_MW * 1e6 / (np.sqrt(3.0) * U_NOM_V)

    i_fault = i_fault_source + i_fault_gen

    return i_fault


# ── Comparison ─────────────────────────────────────────────────────────────────

def compare(pgm_out: dict, net, fault_bus: int, ikss_analytical_A: float) -> None:
    ikss_pgm = float(pgm_out["fault"]["i_f"][0][0])           # phase-A [A]
    ikss_pp  = float(net.res_bus_sc.at[fault_bus, "ikss_ka"]) * 1e3  # kA → A

    ref  = ikss_analytical_A
    rows = [
        ("power-grid-model", ikss_pgm),
        ("pandapower",       ikss_pp),
        ("analytical",       ref),
    ]

    sep = "─" * 52
    print(f"\nFault current — 3-phase bolted, IEC 60909 max (c={C_MAX})")
    print(sep)
    print(f"  {'Tool':<22} {'Ikss [A]':>10}  {'Δ vs analytical':>14}")
    print(f"  {sep}")
    for label, val in rows:
        diff = f"{(val - ref) / ref * 100:>+.3f}%" if label != "analytical" else "reference"
        print(f"  {label:<22} {val:>10.1f}  {diff:>14}")
    print(f"  {sep}")


# ── Main ───────────────────────────────────────────────────────────────────────

def main() -> None:
    pgm_out  = run_pgm()
    net, pp_fault_bus = run_pandapower()
    ikss_an = analytical_ikss()
    compare(pgm_out, net, pp_fault_bus, ikss_an)


if __name__ == "__main__":
    main()

gives

Fault current — 3-phase bolted, IEC 60909 max (c=1.1)
────────────────────────────────────────────────────
  Tool                     Ikss [A]  Δ vs analytical
  ────────────────────────────────────────────────────
  power-grid-model           4808.3         -0.000%
  pandapower                 4808.3         +0.000%
  analytical                 4808.3       reference
  ────────────────────────────────────────────────────