Fixed 3ph load/sgen mapping. e2nIEE#96

pp2ppc_zero lacked connectivity checks, which resulted in wrong net._is_elements entries. Now, all sequences are handled by pd2ppc using the sequence argument. All tests pass except for test_transformer_3ph_diff_kvabase, which has to be investigated.
panza238 · Jul 20, 2018 · c6663c8 · c6663c8
1 parent 5a50b93
commit c6663c8
Show file tree

Hide file tree

Showing 7 changed files with 845 additions and 591 deletions.
diff --git a/pandapower/auxiliary.py b/pandapower/auxiliary.py
@@ -32,6 +32,7 @@
 import pandas as pd
 import scipy as sp
 import six
+from numpy import complex128
 
 from pandapower.idx_brch import F_BUS, T_BUS, BR_STATUS
 from pandapower.idx_bus import BUS_I, BUS_TYPE, NONE, PD, QD, VM, VA, REF
@@ -696,8 +697,20 @@ def I_from_V(Y, V):
 # =============================================================================
 # Calculating Power
 # =============================================================================
-def S_from_VI(V, I):
+def S_from_VI_elementwise(V, I):
     return np.multiply(V, I.conjugate())
 
-def I_from_SV(S, V):
+def I_from_SV_elementwise(S, V):
     return np.conjugate(np.divide(S, V, out=np.zeros_like(S), where=V!=0)) # Return zero if div by zero
+
+def SVabc_from_SV012(S012, V012, n_res=None, idx=None):
+    if n_res==None:
+        n_res = S012.shape[1]
+    if idx==None:
+        idx = np.ones(n_res, dtype="bool")
+    I012 = np.matrix(np.zeros((3, n_res), dtype=complex))
+    I012[:, idx] = I_from_SV_elementwise(S012[:, idx], V012[:, idx])
+    Vabc = sequence_to_phase(V012[:, idx])
+    Iabc = sequence_to_phase(I012[:, idx])
+    Sabc = S_from_VI_elementwise(Vabc, Iabc)
+    return Sabc, Vabc
diff --git a/pandapower/build_bus.py b/pandapower/build_bus.py
@@ -9,6 +9,7 @@
 
 import numpy as np
 import pandas as pd
+from numpy import complex128
 
 from pandapower.auxiliary import _sum_by_group
 from pandapower.idx_bus import BUS_I, BASE_KV, PD, QD, GS, BS, VMAX, VMIN, BUS_TYPE, NONE, VM, VA, CID, CZD, bus_cols
@@ -276,7 +277,7 @@ def _build_bus_ppc(net, ppc):
     net["_pd2ppc_lookups"]["bus"] = bus_lookup
 
 
-def _calc_pq_elements_and_add_on_ppc(net, ppc):
+def _calc_pq_elements_and_add_on_ppc(net, ppc, sequence=None):
     # init values
     b, p, q = np.array([], dtype=int), np.array([]), np.array([])
 
@@ -287,8 +288,11 @@ def _calc_pq_elements_and_add_on_ppc(net, ppc):
     # distinguish calculation modes
     mode = net["_options"]["mode"]
 
+    # Determine sequence
+    pos_sequence = sequence==None or sequence==1
+
     # if mode == powerflow...
-    if mode == "pf" or mode == "pf_3ph":
+    if (mode == "pf" or mode == "pf_3ph") and pos_sequence:
         l = net["load"]
         if len(l) > 0:
             voltage_depend_loads = net["_options"]["voltage_depend_loads"]
@@ -390,6 +394,7 @@ def _calc_pq_elements_and_add_on_ppc(net, ppc):
         b, vp, vq = _sum_by_group(b, p, q)
         ppc["bus"][b, PD] = vp
         ppc["bus"][b, QD] = vq
+        # Todo: Actually, P and Q have to be divided by 3 because Sabc=3*S012 (we are writing pos. seq. values here!)
 
 
 def _calc_shunts_and_add_on_ppc(net, ppc):

diff --git a/pandapower/pd2ppc.py b/pandapower/pd2ppc.py
@@ -5,12 +5,13 @@
 
 
 import copy
+import math
 
 import numpy as np
 
 from pandapower.idx_area import PRICE_REF_BUS
-from pandapower.idx_brch import F_BUS, T_BUS, BR_STATUS
-from pandapower.idx_bus import NONE, BUS_I, BUS_TYPE
+from pandapower.idx_brch import F_BUS, T_BUS, BR_STATUS, branch_cols, TAP, SHIFT, BR_R, BR_X, BR_B
+from pandapower.idx_bus import NONE, BUS_I, BUS_TYPE, BASE_KV, GS, BS
 from pandapower.idx_gen import GEN_BUS, GEN_STATUS
 
 try:
@@ -21,7 +22,8 @@
 
 import pandapower.auxiliary as aux
 from pandapower.build_branch import _build_branch_ppc, _switch_branches, _branches_with_oos_buses, \
-    _update_trafo_trafo3w_ppc
+    _update_trafo_trafo3w_ppc, _initialize_branch_lookup, _calc_tap_from_dataframe, _calc_nominal_ratio_from_dataframe, \
+    _transformer_correction_factor
 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
 from pandapower.build_gen import _build_gen_ppc, _update_gen_ppc
@@ -63,13 +65,13 @@ def _pd2ppc(net, sequence=None):
         **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)
+    net["_is_elements"] = aux._select_is_elements_numba(net, sequence=sequence)
 
     # get options
     mode = net["_options"]["mode"]
     check_connectivity = net["_options"]["check_connectivity"]
 
-    ppc = _init_ppc(net, sequence)
+    ppc = _init_ppc(net, sequence=sequence)
 
     if mode == "opf":
         # additional fields in ppc
@@ -81,14 +83,19 @@ def _pd2ppc(net, sequence=None):
     _build_bus_ppc(net, ppc)
     # generate ppc['gen'] and fills ppc['bus'] with generator values (PV, REF nodes)
     _build_gen_ppc(net, ppc)
+    if sequence==0:
+    # Add external grid impedance to the zero system for 3ph and sc calculations
+        _add_ext_grid_sc_impedance_zero(net, ppc)
     # generate ppc['branch'] and directly generates branch values
-    _build_branch_ppc(net, ppc)
+        _build_branch_ppc_zero(net, ppc)
+    else:
+        _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_motor_impedances_ppc(net, ppc)
     else:
-        _calc_pq_elements_and_add_on_ppc(net, ppc)
+        _calc_pq_elements_and_add_on_ppc(net, ppc, sequence=sequence)
         # adds P and Q for shunts, wards and xwards (to PQ nodes)
         _calc_shunts_and_add_on_ppc(net, ppc)
 
@@ -288,7 +295,7 @@ def _update_ppc(net, sequence=None):
     ppc = net["_ppc"] if sequence is None else net["_ppc%s"% sequence]
     ppci = copy.deepcopy(ppc)
     # adds P and Q for loads / sgens in ppc['bus'] (PQ nodes)
-    _calc_pq_elements_and_add_on_ppc(net, ppc)
+    _calc_pq_elements_and_add_on_ppc(net, ppc, sequence=sequence)
     # adds P and Q for shunts, wards and xwards (to PQ nodes)
     _calc_shunts_and_add_on_ppc(net, ppc)
     # updates values for gen
@@ -309,3 +316,236 @@ def _update_ppc(net, sequence=None):
     ppci["gen"] = ppc["gen"][gs]
 
     return ppc, ppci
+
+
+def _build_branch_ppc_zero(net, ppc):
+    """
+    Takes the empty ppc network and fills it with the zero imepdance branch values. The branch
+    datatype will be np.complex 128 afterwards.
+
+    .. note:: The order of branches in the ppc is:
+            1. Lines
+            2. Transformers
+
+    **INPUT**:
+        **net** -The pandapower format network
+
+        **ppc** - The PYPOWER format network to fill in values
+
+    """
+    length = _initialize_branch_lookup(net)
+    lookup = net._pd2ppc_lookups["branch"]
+    mode = net._options["mode"]
+    ppc["branch"] = np.zeros(shape=(length, branch_cols), dtype=np.complex128)
+    if mode == "sc":
+        from pandapower.shortcircuit.idx_brch import branch_cols_sc
+        branch_sc = np.empty(shape=(length, branch_cols_sc), dtype=float)
+        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 NotImplemented("Three winding transformers are not implemented for unbalanced calculations")
+
+
+def _add_trafo_sc_impedance_zero(net, ppc, trafo_df=None):
+    if trafo_df is None:
+        trafo_df = net["trafo"]
+    branch_lookup = net["_pd2ppc_lookups"]["branch"]
+    if not "trafo" in branch_lookup:
+        return
+    bus_lookup = net["_pd2ppc_lookups"]["bus"]
+    mode = net["_options"]["mode"]
+    f, t = branch_lookup["trafo"]
+    trafo_df["_ppc_idx"] = range(f, t)
+    bus_lookup = net["_pd2ppc_lookups"]["bus"]
+    buses_all, gs_all, bs_all = np.array([], dtype=int), np.array([]), np.array([])
+    for vector_group, trafos in trafo_df.groupby("vector_group"):
+        ppc_idx = trafos["_ppc_idx"].values.astype(int)
+        ppc["branch"][ppc_idx, BR_STATUS] = 0
+
+        if vector_group in ["Yy", "Yd", "Dy", "Dd"]:
+            continue
+
+        vsc_percent = trafos["vsc_percent"].values.astype(float)
+        vscr_percent = trafos["vscr_percent"].values.astype(float)
+        trafo_kva = trafos["sn_kva"].values.astype(float)
+        vsc0_percent = trafos["vsc0_percent"].values.astype(float)
+        vscr0_percent = trafos["vscr0_percent"].values.astype(float)
+        lv_buses = trafos["lv_bus"].values.astype(int)
+        hv_buses = trafos["hv_bus"].values.astype(int)
+        lv_buses_ppc = bus_lookup[lv_buses]
+        hv_buses_ppc = bus_lookup[hv_buses]
+        mag0_percent = trafos.mag0_percent.values.astype(float)
+        mag0_rx = trafos["mag0_rx"].values.astype(float)
+        si0_hv_partial = trafos.si0_hv_partial.values.astype(float)
+        parallel = trafos.parallel.values.astype(float)
+        in_service = trafos["in_service"].astype(int)
+
+        ppc["branch"][ppc_idx, F_BUS] = hv_buses_ppc
+        ppc["branch"][ppc_idx, T_BUS] = lv_buses_ppc
+
+        vn_trafo_hv, vn_trafo_lv, shift = _calc_tap_from_dataframe(net, trafos)
+        #        if mode == 'pf3ph':
+        #            vn_trafo_hv = vn_trafo_hv/np.sqrt(3)
+        #            vn_trafo_lv = vn_trafo_lv/np.sqrt(3)
+        vn_lv = ppc["bus"][lv_buses_ppc, BASE_KV]
+        ratio = _calc_nominal_ratio_from_dataframe(ppc, trafos, vn_trafo_hv, vn_trafo_lv,
+                                                   bus_lookup)
+        ppc["branch"][ppc_idx, TAP] = ratio
+        ppc["branch"][ppc_idx, SHIFT] = shift
+
+        # zero seq. transformer impedance
+        tap_lv = np.square(vn_trafo_lv / vn_lv)  # adjust for low voltage side voltage converter
+        if mode == 'pf_3ph':
+            tap_lv = np.square(vn_trafo_lv / (vn_lv / np.sqrt(3)))
+        z_sc = vsc0_percent / 100. / trafo_kva * tap_lv * net.sn_kva
+        r_sc = vscr0_percent / 100. / trafo_kva * tap_lv * net.sn_kva
+
+        z_sc = z_sc.astype(float)
+        r_sc = r_sc.astype(float)
+        x_sc = np.sign(z_sc) * np.sqrt(z_sc ** 2 - r_sc ** 2)
+        z0_k = (r_sc + x_sc * 1j) / parallel
+        if mode == "sc":
+            from pandapower.shortcircuit.idx_bus import C_MAX
+            cmax = ppc["bus"][lv_buses_ppc, C_MAX]
+            kt = _transformer_correction_factor(vsc_percent, vscr_percent, trafo_kva, cmax)
+            z0_k *= kt
+        y0_k = 1 / z0_k
+        # zero sequence transformer magnetising impedance
+        z_m = ((vsc0_percent * mag0_percent) / 100.) * (net.sn_kva / trafo_kva) * tap_lv
+        x_m = z_m / np.sqrt(mag0_rx ** 2 + 1)
+        r_m = x_m * mag0_rx
+        r0_trafo_mag = r_m / parallel
+        x0_trafo_mag = x_m / parallel
+        z0_mag = r0_trafo_mag + x0_trafo_mag * 1j
+
+        if vector_group == "Dyn":
+            buses_all = np.hstack([buses_all, lv_buses_ppc])
+            gs_all = np.hstack([gs_all, y0_k.real * in_service * int(ppc["baseMVA"])])
+            bs_all = np.hstack([bs_all, y0_k.imag * in_service * int(ppc["baseMVA"])])
+
+        elif vector_group == "YNd":
+            buses_all = np.hstack([buses_all, hv_buses_ppc])
+            gs_all = np.hstack([gs_all, y0_k.real * in_service * int(ppc["baseMVA"])])
+            bs_all = np.hstack([bs_all, y0_k.imag * in_service * int(ppc["baseMVA"])])
+
+        elif vector_group == "Yyn":
+            buses_all = np.hstack([buses_all, lv_buses_ppc])
+            y = 1 / (z0_mag + z0_k).astype(complex) * int(ppc["baseMVA"])
+            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
+            # convert the t model to pi model
+            z1 = si0_hv_partial * z0_k
+            z2 = (1 - si0_hv_partial) * z0_k
+            z3 = z0_mag
+
+            z_temp = z1 * z2 + z2 * z3 + z1 * z3
+            za = z_temp / z2
+            zb = z_temp / z1
+            zc = z_temp / z3
+
+            ppc["branch"][ppc_idx, BR_R] = zc.real
+            ppc["branch"][ppc_idx, BR_X] = zc.imag
+            y = 2 / za
+            ppc["branch"][ppc_idx, BR_B] = y.imag - y.real * 1j
+            # add a shunt element parallel to zb if the leakage impedance distribution is unequal
+            # TODO: this only necessary if si0_hv_partial!=0.5 --> test
+            zs = (za * zb) / (za - zb)
+            ys = 1 / zs.astype(complex)
+            buses_all = np.hstack([buses_all, lv_buses_ppc])
+            gs_all = np.hstack([gs_all, ys.real * in_service * int(ppc["baseMVA"])])
+            bs_all = np.hstack([bs_all, ys.imag * in_service * int(ppc["baseMVA"])])
+        elif vector_group == "YNy":
+            buses_all = np.hstack([buses_all, hv_buses_ppc])
+            y = 1 / (z0_mag + z0_k).astype(complex) * int(ppc["baseMVA"])
+            gs_all = np.hstack([gs_all, y.real * in_service])
+            bs_all = np.hstack([bs_all, y.imag * in_service])
+        elif vector_group[-1].isdigit():
+            raise ValueError(
+                "Unknown transformer vector group %s - please specify vector group without phase shift number. Phase shift can be specified in net.trafo.shift_degree" % vector_group)
+        else:
+            raise ValueError("Transformer vector group %s is unknown / not implemented" % vector_group)
+
+    buses, gs, bs = aux._sum_by_group(buses_all, gs_all, bs_all)
+    ppc["bus"][buses, GS] += gs
+    ppc["bus"][buses, BS] += bs
+    del net.trafo["_ppc_idx"]
+
+
+def _add_ext_grid_sc_impedance_zero(net, ppc):
+    mode = net["_options"]["mode"]
+
+    if mode == "sc":
+        from pandapower.shortcircuit.idx_bus import C_MAX, C_MIN
+        case = net._options["case"]
+    else:
+        case = "max"
+    bus_lookup = net["_pd2ppc_lookups"]["bus"]
+    eg = net["ext_grid"][net._is_elements["ext_grid"]]
+    if len(eg) == 0:
+        return
+    eg_buses = eg.bus.values
+    eg_buses_ppc = bus_lookup[eg_buses]
+
+    if mode == "sc":
+        c = ppc["bus"][eg_buses_ppc, C_MAX] if case == "max" else ppc["bus"][eg_buses_ppc, C_MIN]
+    elif mode == 'pf_3ph':
+        c = 1.1  # Todo: Where does that value come from?
+    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")
+    s_sc = eg["s_sc_%s_mva" % case].values
+    if not "rx_%s" % case in eg:
+        raise ValueError("short circuit R/X rate rx_%s needs to be specified for external grid" %
+                         case)
+    rx = eg["rx_%s" % case].values
+    z_grid = c / s_sc
+    if mode == 'pf_3ph':
+        z_grid = c / (s_sc / 3)  # 3 phase power divided to get 1 ph power
+    x_grid = z_grid / np.sqrt(rx ** 2 + 1)
+    r_grid = rx * x_grid
+    eg["r"] = r_grid
+    eg["x"] = x_grid
+
+    # ext_grid zero sequence impedance
+    if case == "max":
+        x0_grid = net.ext_grid["x0x_%s" % case] * x_grid
+        r0_grid = net.ext_grid["r0x0_%s" % case] * x0_grid
+    elif case == "min":
+        x0_grid = net.ext_grid["x0x_%s" % case] * x_grid
+        r0_grid = net.ext_grid["r0x0_%s" % case] * x0_grid
+    y0_grid = 1 / (r0_grid + x0_grid * 1j)
+    buses, gs, bs = aux._sum_by_group(eg_buses_ppc, y0_grid.real, y0_grid.imag)
+    ppc["bus"][buses, GS] = gs
+    ppc["bus"][buses, BS] = bs
+
+
+def _add_line_sc_impedance_zero(net, ppc):
+    branch_lookup = net["_pd2ppc_lookups"]["branch"]
+    mode = net["_options"]["mode"]
+    if not "line" in branch_lookup:
+        return
+    line = net["line"]
+    bus_lookup = net["_pd2ppc_lookups"]["bus"]
+    length = line["length_km"].values
+    parallel = line["parallel"].values
+
+    fb = bus_lookup[line["from_bus"].values]
+    tb = bus_lookup[line["to_bus"].values]
+    baseR = np.square(ppc["bus"][fb, BASE_KV]) / ppc["baseMVA"]
+    if mode == 'pf_3ph':
+        baseR = np.square(ppc["bus"][fb, BASE_KV] / np.sqrt(3)) / ppc["baseMVA"]
+    f, t = branch_lookup["line"]
+    # line zero sequence impedance
+    ppc["branch"][f:t, F_BUS] = fb
+    ppc["branch"][f:t, T_BUS] = tb
+    ppc["branch"][f:t, BR_R] = line["r0_ohm_per_km"].values * length / baseR / parallel
+    ppc["branch"][f:t, BR_X] = line["x0_ohm_per_km"].values * length / baseR / parallel
+    ppc["branch"][f:t, BR_B] = (
+                2 * net["f_hz"] * math.pi * line["c0_nf_per_km"].values * 1e-9 * baseR * length * parallel)
+    ppc["branch"][f:t, BR_STATUS] = line["in_service"].astype(int)