Feature/fixing transformer kt min case

pandapower/build_branch.py

@@ -476,13 +476,11 @@ def _calc_r_x_y_from_dataframe(net, trafo_df, vn_trafo_lv, vn_lv, ppc, sequence=
         else:
             g, b = 0, 0  # why for sc are we assigning y directly as 0?
         if isinstance(trafo_df, pd.DataFrame):  # 2w trafo is dataframe, 3w trafo is dict
-            bus_lookup = net._pd2ppc_lookups["bus"]
-            cmax = ppc["bus"][bus_lookup[net.trafo.lv_bus.values], C_MAX]
-            # todo: kt is only used for case = max and only for network transformers! (IEC 60909-0:2016 section 6.3.3)
-            # kt is only calculated for network transformers (IEC 60909-0:2016 section 6.3.3)
             if not net._options.get("use_pre_fault_voltage", False):
-                kt = _transformer_correction_factor(
-                    trafo_df, trafo_df.vk_percent, trafo_df.vkr_percent, trafo_df.sn_mva, cmax)
+                bus_lookup = net._pd2ppc_lookups["bus"]
+                cmax = ppc["bus"][bus_lookup[net.trafo.lv_bus.values], C_MAX]
+                case = net._options["case"]
+                kt = _transformer_correction_factor(trafo_df, trafo_df.vk_percent, trafo_df.vkr_percent, trafo_df.sn_mva, cmax, case)
                 r *= kt
                 x *= kt
     else:
@@ -1321,7 +1319,7 @@ def _end_temperature_correction_factor(net, short_circuit=False, dc=False):
     return r_correction_for_temperature
 
 
-def _transformer_correction_factor(trafo_df, vk, vkr, sn, cmax):
+def _transformer_correction_factor(trafo_df, vk, vkr, sn, cmax, case):
     """
         2W-Transformer impedance correction factor in short circuit calculations,
         based on the IEC 60909-0:2016 standard.
@@ -1330,6 +1328,7 @@ def _transformer_correction_factor(trafo_df, vk, vkr, sn, cmax):
             vkr: real-part of transformer short-circuit voltage, percent
             sn: transformer rating, kVA
             cmax: voltage factor to account for maximum worst-case currents, based on the lv side
+            case: short-circuit calculation case (str, "min"/"max")
 
         Returns:
             kt: transformer impedance correction factor for short-circuit calculations
@@ -1338,7 +1337,12 @@ def _transformer_correction_factor(trafo_df, vk, vkr, sn, cmax):
     ----------
     trafo_df
 
-        """
+    """
+
+    # The transformer correction factor shall only be applied in the max case according to
+    # norm IEC 60909-0:2016 section 6.3.3
+    if case != "max":
+        return np.ones(len(trafo_df))
 
     if "power_station_unit" in trafo_df.columns:
         power_station_unit = trafo_df.power_station_unit.fillna(False).values.astype(bool)
@@ -1370,6 +1374,7 @@ def _trafo_df_from_trafo3w(net, sequence=1):
     loss_side = t3.loss_side.values if "loss_side" in t3.columns else np.full(len(t3),
                                                                               net._options["trafo3w_losses"].lower())
     nr_trafos = len(net["trafo3w"])
+
     if sequence == 1:
         if 'tap_dependency_table' in t3:
             mode_tmp = "type_c" if mode == "sc" and net._options.get("use_pre_fault_voltage", False) else mode
@@ -1382,7 +1387,8 @@ def _trafo_df_from_trafo3w(net, sequence=1):
         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,)
+        case = net._options.get("case", 'max')
+        _calculate_sc_voltages_of_equivalent_transformers_zero_sequence(t3, trafo2, case=case)
     else:
         raise UserWarning("Unsupported sequence for trafo3w convertion")
     _calculate_3w_tap_changers(t3, trafo2, sides)
@@ -1416,8 +1422,7 @@ def _trafo_df_from_trafo3w(net, sequence=1):
     return {var: np.concatenate([trafo2[var][side] for side in sides]) for var in trafo2.keys()}
 
 
-def _calculate_sc_voltages_of_equivalent_transformers(
-        t3, t2, mode, characteristic=None, net=None):
+def _calculate_sc_voltages_of_equivalent_transformers(t3, t2, mode, characteristic=None, net=None):
     if "tap_dependency_table" in t3:
         tap_dependency_table = get_trafo_values(t3, "tap_dependency_table")
         tap_dependency_table = np.array(
@@ -1444,7 +1449,8 @@ def _calculate_sc_voltages_of_equivalent_transformers(
     vk_2w_delta = z_br_to_bus_vector(vk_3w, sn)
     vkr_2w_delta = z_br_to_bus_vector(vkr_3w, sn)
     if mode == "sc":
-        kt = _transformer_correction_factor(t3, vk_3w, vkr_3w, sn, 1.1)
+        case = net._options.get("case", 'max')
+        kt = _transformer_correction_factor(t3, vk_3w, vkr_3w, sn, 1.1, case)
         vk_2w_delta *= kt
         vkr_2w_delta *= kt
     vki_2w_delta = np.sqrt(vk_2w_delta ** 2 - vkr_2w_delta ** 2)
@@ -1458,7 +1464,7 @@ def _calculate_sc_voltages_of_equivalent_transformers(
     t2["sn_mva"] = {"hv": sn[0, :], "mv": sn[1, :], "lv": sn[2, :]}
 
 
-def _calculate_sc_voltages_of_equivalent_transformers_zero_sequence(t3, t2):
+def _calculate_sc_voltages_of_equivalent_transformers_zero_sequence(t3, t2, case):
     vk_3w = np.stack([t3.vk_hv_percent.values, t3.vk_mv_percent.values, t3.vk_lv_percent.values])
     vkr_3w = np.stack([t3.vkr_hv_percent.values, t3.vkr_mv_percent.values, t3.vkr_lv_percent.values])
     vk0_3w = np.stack([t3.vk0_hv_percent.values, t3.vk0_mv_percent.values, t3.vk0_lv_percent.values])
@@ -1469,7 +1475,7 @@ def _calculate_sc_voltages_of_equivalent_transformers_zero_sequence(t3, t2):
     vkr0_2w_delta = z_br_to_bus_vector(vkr0_3w, sn)
 
     # Only for "sc", calculated with positive sequence value
-    kt = _transformer_correction_factor(t3, vk_3w, vkr_3w, sn, 1.1)
+    kt = _transformer_correction_factor(t3, vk_3w, vkr_3w, sn, 1.1, case)
     vk0_2w_delta *= kt
     vkr0_2w_delta *= kt
 
@@ -1522,8 +1528,7 @@ def _calculate_3w_tap_changers(t3, t2, sides):
 
         # t3 trafos with tap changer at star points
         if any_at_star_point & np.any(mask_star_point := (tap_mask & at_star_point)):
-            t = (tap_arrays["tap_step_percent"][side][mask_star_point] *
-                 np.exp(1j * np.deg2rad(tap_arrays["tap_step_degree"][side][mask_star_point])))
+            t = tap_arrays["tap_step_percent"][side][mask_star_point] * np.exp(1j * np.deg2rad(tap_arrays["tap_step_degree"][side][mask_star_point]))
             tap_pos = tap_arrays["tap_pos"][side][mask_star_point]
             tap_neutral = tap_arrays["tap_neutral"][side][mask_star_point]
             t_corrected = 100 * t / (100 + (t * (tap_pos-tap_neutral)))

pandapower/build_bus.py

@@ -963,12 +963,12 @@ def _add_ext_grid_sc_impedance(net, ppc):
     else:
         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 "
-                          "external grid \n Try: net.ext_grid['s_sc_max_mva'] = 1000") % case)
+        raise ValueError(f"short circuit apparent power s_sc_{case}_mva needs to be specified for external grid \n"
+                         f" Try: net.ext_grid['s_sc_{case}_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)
+        raise ValueError(f"short circuit R/X rate rx_{case} needs to be specified for external grid \n"
+                         f" Try: net.ext_grid['rx_{case}'] = 0.1")
     rx = eg["rx_%s" % case].values
 
     z_grid = c / s_sc

pandapower/pd2ppc_zero.py

@@ -285,9 +285,8 @@ def _add_trafo_sc_impedance_zero(net, ppc, trafo_df=None, k_st=None):
 
         if mode == "sc":  # or trafo_model == "pi":
             cmax = net._ppc["bus"][lv_buses_ppc, C_MAX]
-            kt = _transformer_correction_factor(
-                trafos, vk_percent, vkr_percent, sn_trafo_mva, cmax
-            )
+            case = net._options["case"]
+            kt = _transformer_correction_factor(trafos, vk_percent, vkr_percent, sn_trafo_mva, cmax, case)
             z0_k *= kt
 
             # different formula must be applied for power station unit transformers:

pandapower/shortcircuit/ppc_conversion.py

@@ -70,7 +70,8 @@ def _add_kt(net, ppc):
         f, t = net["_pd2ppc_lookups"]["branch"]["trafo"]
         trafo_df = net["trafo"]
         cmax = ppc["bus"][bus_lookup[get_trafo_values(trafo_df, "lv_bus")], C_MAX]
-        kt = _transformer_correction_factor(trafo_df, trafo_df.vk_percent, trafo_df.vkr_percent, trafo_df.sn_mva, cmax)
+        case = net._options["case"] if net._options.get("mode", "") == "sc" else None
+        kt = _transformer_correction_factor(trafo_df, trafo_df.vk_percent, trafo_df.vkr_percent, trafo_df.sn_mva, cmax, case)
         branch[f:t, K_T] = kt
 
 

pandapower/test/shortcircuit/test_1ph.py

@@ -119,16 +119,26 @@ def test_1ph_shortcircuit_3w():
 
 def test_1ph_shortcircuit_min():
     results = {
-        "Yy": [0.52209346201, 0.66632662571, 0.66756160176, 0.72517293174],
-        "Yyn": [0.52209346201, 2.4135757259, 1.545054139, 0.99373917957],
-        "Yd": [0.52209346201, 0.66632662571, 0.66756160176, 0.72517293174],
-        "YNy": [0.62316686505, 0.66632662571, 0.66756160176, 0.72517293174],
-        "YNyn": [0.620287259, 2.9155736491, 1.7561556936, 1.0807305212],
-        "YNd": [0.75434229157, 0.66632662571, 0.66756160176, 0.72517293174],
-        "Dy": [0.52209346201, 0.66632662571, 0.66756160176, 0.72517293174],
-        "Dyn": [0.52209346201, 3.4393798093, 1.9535982949, 1.1558364456],
-        "Dd": [0.52209346201, 0.66632662571, 0.66756160176, 0.72517293174]
-    }
+                 "Yy":  [0.52209345151, 0.66680074537, 0.66803609106, 0.72571559085]
+                ,"Yyn": [0.52209345151, 2.40601861813, 1.54211362172, 0.99258179645]
+                ,"Yd":  [0.52209345151, 0.66680074537, 0.66803609106, 0.72571559085]
+                ,"YNy": [0.62305132639, 0.66680074537, 0.66803609106, 0.72571559085]
+                ,"YNyn":[0.62019004941, 2.90615348290, 1.75300032937, 1.07961489722]
+                ,"YNd": [0.75372631070, 0.66680074537, 0.66803609106, 0.72571559085]
+                ,"Dy":  [0.52209345151, 0.66680074537, 0.66803609106, 0.72571559085]
+                ,"Dyn": [0.52209345151, 3.42429321754, 1.94906060941, 1.15434123431]
+                ,"Dd":  [0.52209345151, 0.66680074537, 0.66803609106, 0.72571559085]
+               }
+#        "Yy": [0.52209346201, 0.66632662571, 0.66756160176, 0.72517293174],
+#        "Yyn": [0.52209346201, 2.4135757259, 1.545054139, 0.99373917957],
+#        "Yd": [0.52209346201, 0.66632662571, 0.66756160176, 0.72517293174],
+#        "YNy": [0.62316686505, 0.66632662571, 0.66756160176, 0.72517293174],
+#        "YNyn": [0.620287259, 2.9155736491, 1.7561556936, 1.0807305212],
+#        "YNd": [0.75434229157, 0.66632662571, 0.66756160176, 0.72517293174],
+#        "Dy": [0.52209346201, 0.66632662571, 0.66756160176, 0.72517293174],
+#        "Dyn": [0.52209346201, 3.4393798093, 1.9535982949, 1.1558364456],
+#        "Dd": [0.52209346201, 0.66632662571, 0.66756160176, 0.72517293174]
+#    }
 
     for inv_y in (False, True):
         for vc, result in results.items():
@@ -338,7 +348,7 @@ def test_iec60909_example_4_two_trafo3w_two_earth():
     # r = 0.59621204768
     # x = 6.0598429694
     ikss_pf_max = [26.0499, 20.9472, 9.1722, 18.7457]
-    ikss_pf_min = [3.9622, 8.3914, 5.0248, 6.9413]
+    ikss_pf_min = [3.9447, 8.2164, 4.9642, 6.8273]
 
     calc_sc(net, fault="1ph", case="max")
     assert np.allclose(net.res_bus_sc.ikss_ka.values[:4], np.array(ikss_pf_max), atol=1e-4)

pandapower/test/shortcircuit/test_all_currents.py

@@ -707,61 +707,98 @@ def test_trafo_3w():
     pass
 
 
-def test_trafo_impedance():
+# Todo: "min" case does not work, since parameters are missing.
+@pytest.mark.parametrize("trafo_impedance_case", ["max", "min"])
+def test_trafo_impedance(trafo_impedance_case):
+    case = trafo_impedance_case
+
     net = create_empty_network(sn_mva=0.16)
-    create_bus(net, 20)
-    create_buses(net, 2, 0.4)
-    create_ext_grid(net, 0, s_sc_max_mva=346.4102, rx_max=0.1)
+    b0 = create_bus(net, 20)
+    b1, b2 = create_buses(net, 2, 0.4)
+    create_ext_grid(net, 0, s_sc_max_mva=100, s_sc_min_mva=80, rx_max=0.1, rx_min=0.1)
     v_lv = 410
-    create_transformer_from_parameters(net, 0, 1, 0.4, 20, v_lv / 1e3, 1.15, 4, 0, 0)
-    create_line_from_parameters(net, 1, 2, 0.004, 0.208, 0.068, 0, 1, parallel=2)
+    create_transformer_from_parameters(
+        net=net,
+        hv_bus=0,
+        lv_bus=1,
+        sn_mva=1,
+        vn_hv_kv=20,
+        vn_lv_kv=0.4,
+        vkr_percent=0.96,
+        vk_percent=6,
+        pfe_kw=0,
+        i0_percent=0
+    )
+    create_line_from_parameters(
+        net=net,
+        from_bus=1,
+        to_bus=2,
+        length_km=0.5,
+        r_ohm_per_km=0.208,
+        x_ohm_per_km=0.068,
+        c_nf_per_km=1,
+        max_i_ka=1,
+        parallel=1,
+        endtemp_degree=20
+    )
     # create_load(net, 2, 0.1)
 
-    runpp(net)
-
-    calc_sc(net, case='max', lv_tol_percent=6., bus=2, branch_results=True, use_pre_fault_voltage=False)
-
+    # runpp(net)
+
     # trafo:
-    z_tlv = 4 / 100 * v_lv ** 2 / (400 * 1e3)
-    r_tlv = 4600 * v_lv ** 2 / ((400 * 1e3) ** 2)
-    x_tlv = np.sqrt(z_tlv ** 2 - r_tlv ** 2)
-    z_tlv = r_tlv + 1j * x_tlv
-    x_t = x_tlv * 400 * 1e3 / (v_lv ** 2)
-    k_t = 0.95 * 1.05 / (1 + 0.6 * x_t)
-    z_tk = k_t * z_tlv
+    # v_lv = 410
+    # z_tlv = 4 / 100 * v_lv ** 2 / (400 * 1e3)
+    # r_tlv = 4600 * v_lv ** 2 / ((400 * 1e3) ** 2)
+    # x_tlv = np.sqrt(z_tlv ** 2 - r_tlv ** 2)
+    # z_tlv = r_tlv + 1j * x_tlv
+    # x_t = x_tlv * 400 * 1e3 / (v_lv ** 2)
+    # k_t = 0.95 * 1.05 / (1 + 0.6 * x_t)
+    # z_tk = k_t * z_tlv
+
+    calc_sc(net, case=case, lv_tol_percent=6., bus=2, branch_results=True, use_pre_fault_voltage=False)
+
+    # min test case
+    if trafo_impedance_case == "min":
+        assert np.allclose(net.res_bus_sc.ikss_ka, 1.906175, rtol=0, atol=1e-3)
+        assert np.allclose(net.res_bus_sc.skss_mw, 1.320637, rtol=0, atol=1e-5)
+        assert np.allclose(net.res_bus_sc.rk_ohm, 0.106, rtol=0, atol=1e-3)
+        assert np.allclose(net.res_bus_sc.xk_ohm, 0.045, rtol=0, atol=1e-3)
+
+    # max test case
+    else:
+        assert np.allclose(net.res_bus_sc.ikss_ka, 2.112413, rtol=0, atol=1e-3)
+        assert np.allclose(net.res_bus_sc.skss_mw, 1.463523, rtol=0, atol=1e-5)
+        assert np.allclose(net.res_bus_sc.rk_ohm, 0.106, rtol=0, atol=1e-3)
+        assert np.allclose(net.res_bus_sc.xk_ohm, 0.045, rtol=0, atol=1e-3)
+
+
+    # ppci = net.ppci
+    # tap = ppci["branch"][:, TAP].real
+    # ikss1 = ppci["bus"][:, IKSS1] * np.exp(1j * np.deg2rad(ppci["bus"][:, PHI_IKSS1_DEGREE]))
 
     # line:
-    z_l = 0.416 * 1e-3 + 1j * 0.136 * 1e-3  # Ohm
-
-    # assert np.allclose(net.res_bus_sc.rk_ohm * 1e3, 5.18, rtol=0, atol=1e-6)
-    # assert np.allclose(net.res_bus_sc.xk_ohm * 1e3, 16.37, rtol=0, atol=1e-6)
-    assert np.allclose(net._ppc['branch'][:, BR_R].real, [0.416 * 1e-3, z_tk.real], rtol=0, atol=1e-6)
-    assert np.allclose(net._ppc['branch'][:, BR_X].real, [0.136 * 1e-3, z_tk.imag], rtol=0, atol=1e-6)
-
-    ppci = net.ppci
-    tap = ppci["branch"][:, TAP].real
-    ikss1 = ppci["bus"][:, IKSS1] * np.exp(1j * np.deg2rad(ppci["bus"][:, PHI_IKSS1_DEGREE]))
-
-    v_1 = ikss1[2] * z_l / 0.4 * np.sqrt(3)  # kA * Ohm / V_base -> p.u.
-    np.abs(v_1)
-    np.angle(v_1, deg=True)
-
-    v_0 = v_1 + ikss1[2] * z_tk / 0.4 * np.sqrt(3) * 0.4 / 0.41
-    np.abs(v_0)
-    np.angle(v_0, deg=True)
-
-    v_0_ref = v_1 + ikss1[2] * z_tk / 0.4 * np.sqrt(3)
-
-    Yf = ppci["internal"]["Yf"]
-    Yt = ppci["internal"]["Yt"]
-    V_diff = np.ones_like(net.bus.index.values, dtype=np.complex128)
-    V_diff[0] = v_0
-    V_diff[1] = v_1
-    V_diff[2] = 0
-    i_f = Yf.dot(V_diff) / ppci["internal"]["baseI"][ppci["branch"][:, F_BUS].real.astype(np.int64)]
-    i_t = Yt.dot(V_diff) / ppci["internal"]["baseI"][ppci["branch"][:, T_BUS].real.astype(np.int64)]
-    abs(i_f)
-    abs(i_t)
+    # z_l = 0.416 * 1e-3 + 1j * 0.136 * 1e-3  # Ohm
+    #
+    # v_1 = ikss1[2] * z_l / 0.4 * np.sqrt(3)  # kA * Ohm / V_base -> p.u.
+    # np.abs(v_1)
+    # np.angle(v_1, deg=True)
+    #
+    # v_0 = v_1 + ikss1[2] * z_tk / 0.4 * np.sqrt(3) * 0.4 / 0.41
+    # np.abs(v_0)
+    # np.angle(v_0, deg=True)
+    #
+    # v_0_ref = v_1 + ikss1[2] * z_tk / 0.4 * np.sqrt(3)
+    #
+    # Yf = ppci["internal"]["Yf"]
+    # Yt = ppci["internal"]["Yt"]
+    # V_diff = np.ones_like(net.bus.index.values, dtype=np.complex128)
+    # V_diff[0] = v_0
+    # V_diff[1] = v_1
+    # V_diff[2] = 0
+    # i_f = Yf.dot(V_diff) / ppci["internal"]["baseI"][ppci["branch"][:, F_BUS].real.astype(np.int64)]
+    # i_t = Yt.dot(V_diff) / ppci["internal"]["baseI"][ppci["branch"][:, T_BUS].real.astype(np.int64)]
+    # abs(i_f)
+    # abs(i_t)
 
 
 @pytest.mark.parametrize("inverse_y", (True, False), ids=("Inverse Y", "LU factorization"))

pandapower/test/shortcircuit/test_iec60909_4.py

@@ -349,7 +349,10 @@ def test_iec_60909_4_3ph_min():
     net.ext_grid["rx_min"] = net.ext_grid["rx_max"]
     calc_sc(net, fault="3ph", case="min", ip=True, tk_s=0.1, kappa_method="C")
 
-    ikss_min = [5.0501, 12.2915, 10.3292, 9.4708, 11.8604, 28.3052, 18.6148, 10.9005, 44.5098, 67.9578]
+    ikss_min = [5.0257, 12.0645, 10.2108, 9.3820, 11.6761,
+                27.7655, 18.3930, 10.9024, 44.4310, 67.8216]
+
+    # ikss_min = [5.0501, 12.2915, 10.3292, 9.4708, 11.8604, 28.3052, 18.6148, 10.9005, 44.5098, 67.9578]
 
     assert np.allclose(net.res_bus_sc.ikss_ka.values[:10], np.array(ikss_min), atol=1e-3)