Repository: e2nIEE/pandapower

Issue:
BUG: distributed_slack=True and enforce_q_lims=True cause error in pfsoln

Discussion:

Patch:
diff --git a/pandapower/estimation/results.py b/pandapower/estimation/results.py
index 87d10c91..08ffa31b 100644
--- a/pandapower/estimation/results.py
+++ b/pandapower/estimation/results.py
@@ -17,7 +17,7 @@ def _calc_power_flow(ppci, V):
     # calculate branch results (in ppc_i)
     baseMVA, bus, gen, branch, ref, pv, pq, _, _, _, ref_gens = _get_pf_variables_from_ppci(ppci)
     Ybus, Yf, Yt = ppci['internal']['Ybus'], ppci['internal']['Yf'], ppci['internal']['Yt']
-    ppci['bus'], ppci['gen'], ppci['branch'] =\
+    ppci['bus'], ppci['gen'], ppci['branch'] = \
         pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, ref_gens)
 
     # calculate bus power injections
diff --git a/pandapower/pf/pfsoln_numba.py b/pandapower/pf/pfsoln_numba.py
index 48569700..4c77811e 100644
--- a/pandapower/pf/pfsoln_numba.py
+++ b/pandapower/pf/pfsoln_numba.py
@@ -11,7 +11,7 @@
 """Updates bus, gen, branch data structures to match power flow soln.
 """
 
-from numpy import conj, zeros, complex128, abs, float64, sqrt, real
+from numpy import conj, zeros, complex128, abs, float64, sqrt, real, isin, arange
 from numpy import finfo, c_, flatnonzero as find, setdiff1d, r_
 
 from pandapower.pypower.idx_brch import F_BUS, T_BUS, PF, PT, QF, QT
@@ -27,7 +27,7 @@ except ImportError:
 EPS = finfo(float).eps
 
 
-def pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, ref_gens, Ibus=None):
+def pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, ref_gens, Ibus=None, limited_gens=None):
     """Updates bus, gen, branch data structures to match power flow soln.
 
     @author: Ray Zimmerman (PSERC Cornell)
@@ -42,13 +42,17 @@ def pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, ref_gens, Ibus=None)
 
     # compute total injected bus powers
     Ibus = zeros(len(V)) if Ibus is None else Ibus
-    Sbus = V[gbus] * conj(Ybus[gbus, :] * V - Ibus[gbus])
-    Sbus_xw = V[xbus] * conj(Ybus[xbus, :] * V - Ibus[xbus])
+    Sbus = V * conj(Ybus * V - Ibus)
 
     _update_v(bus, V)
     # update gen results
-    _update_q(baseMVA, bus, gen, gbus, Sbus, on)
-    _update_p(baseMVA, bus, gen, ref, gbus, on, r_[Sbus, Sbus_xw], ref_gens)
+    _update_q(baseMVA, bus, gen, gbus, Sbus[gbus], on)
+
+    if limited_gens is not None and len(limited_gens) > 0:
+        on = find((gen[:, GEN_STATUS] > 0) | isin(arange(len(gen)), limited_gens))
+        gbus = gen[on, GEN_BUS].astype(int)
+
+    _update_p(baseMVA, bus, gen, ref, gbus, Sbus, ref_gens)
 
     # ----- update/compute branch power flows -----
     branch = _update_branch_flows(Yf, Yt, V, baseMVA, branch)
@@ -56,7 +60,7 @@ def pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, ref_gens, Ibus=None)
     return bus, gen, branch
 
 
-def pf_solution_single_slack(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, ref_gens, Ibus=None):
+def pf_solution_single_slack(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, ref_gens, Ibus=None, limited_gens=None):
     """
     faster version of pfsoln for a grid with a single slack bus
 
diff --git a/pandapower/pf/run_newton_raphson_pf.py b/pandapower/pf/run_newton_raphson_pf.py
index d69a09b3..a839fc92 100644
--- a/pandapower/pf/run_newton_raphson_pf.py
+++ b/pandapower/pf/run_newton_raphson_pf.py
@@ -63,7 +63,7 @@ def _run_newton_raphson_pf(ppci, options):
     return ppci
 
 
-def ppci_to_pfsoln(ppci, options):
+def ppci_to_pfsoln(ppci, options, limited_gens=None):
     internal = ppci["internal"]
     if options["only_v_results"]:
         # time series relevant hack which ONLY saves V from ppci
@@ -93,7 +93,7 @@ def ppci_to_pfsoln(ppci, options):
             internal["Ybus"], internal["Yf"], internal["Yt"] = makeYbus(internal["baseMVA"], internal["bus"], internal["branch"])
 
         result_pfsoln = pfsoln(internal["baseMVA"], internal["bus"], internal["gen"], internal["branch"], internal["Ybus"],
-                      internal["Yf"], internal["Yt"], internal["V"], ref, ref_gens)
+                      internal["Yf"], internal["Yt"], internal["V"], ref, ref_gens, limited_gens=limited_gens)
         return result_pfsoln
 
 def _get_Y_bus(ppci, options, makeYbus, baseMVA, bus, branch):
@@ -174,6 +174,7 @@ def _run_ac_pf_without_qlims_enforced(ppci, options):
 def _run_ac_pf_with_qlims_enforced(ppci, options):
     baseMVA, bus, gen, branch, ref, pv, pq, on, _, V0, ref_gens = _get_pf_variables_from_ppci(ppci)
     bus_backup_p_q = bus[:, [PD, QD]].copy()
+    gen_backup_p = gen[:, PG].copy()
 
     qlim = options["enforce_q_lims"]
     limited = []  # list of indices of gens @ Q lims
@@ -181,7 +182,9 @@ def _run_ac_pf_with_qlims_enforced(ppci, options):
 
     while True:
         ppci, success, iterations = _run_ac_pf_without_qlims_enforced(ppci, options)
-        bus, gen, branch = ppci_to_pfsoln(ppci, options)
+        gen[:, PG] = gen_backup_p
+        bus[:, PD] = bus_backup_p_q[:, 0]
+        bus, gen, branch = ppci_to_pfsoln(ppci, options, limited)
 
         # find gens with violated Q constraints
         gen_status = gen[:, GEN_STATUS] > 0
@@ -210,11 +213,6 @@ def _run_ac_pf_with_qlims_enforced(ppci, options):
 
             # convert to PQ bus
             gen[mx, QG] = fixedQg[mx]  # set Qg to binding
-            for i in range(len(mx)):  # [one at a time, since they may be at same bus]
-                gen[mx[i], GEN_STATUS] = 0  # temporarily turn off gen,
-                bi = gen[mx[i], GEN_BUS].astype(int)  # adjust load accordingly,
-                bus[bi, [PD, QD]] = (bus[bi, [PD, QD]] - gen[mx[i], [PG, QG]])
-
             #            if len(ref) > 1 and any(bus[gen[mx, GEN_BUS].astype(int), BUS_TYPE] == REF):
             #                raise ValueError('Sorry, pandapower cannot enforce Q '
             #                                 'limits for slack buses in systems '
@@ -227,6 +225,11 @@ def _run_ac_pf_with_qlims_enforced(ppci, options):
             ref, pv, pq = bustypes(bus, gen)
 
             limited = r_[limited, mx].astype(int)
+
+            for i in range(len(limited)):  # [one at a time, since they may be at same bus]
+                gen[limited[i], GEN_STATUS] = 0  # temporarily turn off gen,
+                bi = gen[limited[i], GEN_BUS].astype(int)  # adjust load accordingly,
+                bus[bi, [PD, QD]] = (bus[bi, [PD, QD]] - gen[limited[i], [PG, QG]])
         else:
             break  # no more generator Q limits violated
 
@@ -234,10 +237,6 @@ def _run_ac_pf_with_qlims_enforced(ppci, options):
         # restore injections from limited gens [those at Q limits]
         bus[setdiff1d(changed_gens, ref), BUS_TYPE] = PV  # & set bus type back to PV
         gen[limited, QG] = fixedQg[limited]  # restore Qg value,
-        for i in range(len(limited)):  # [one at a time, since they may be at same bus]
-            bi = gen[limited[i], GEN_BUS].astype(int)  # re-adjust load,
-            gen[limited[i], PG] = bus_backup_p_q[bi, 0] - bus[bi, PD]
-            bus[bi, [PD, QD]] = bus_backup_p_q[bi, :]
-            gen[limited[i], GEN_STATUS] = 1  # and turn gen back on
-
+        gen[limited, GEN_STATUS] = 1  # turn gens back on
+        bus[:, [PD, QD]] = bus_backup_p_q
     return ppci, success, iterations, bus, gen, branch
diff --git a/pandapower/pypower/pfsoln.py b/pandapower/pypower/pfsoln.py
index 04496cf7..bcf1e40d 100644
--- a/pandapower/pypower/pfsoln.py
+++ b/pandapower/pypower/pfsoln.py
@@ -12,7 +12,7 @@
 """
 
 from numpy import asarray, angle, pi, conj, zeros, ones, finfo, c_, ix_, real, flatnonzero as find, \
-    setdiff1d, intersect1d, r_
+    setdiff1d, intersect1d, r_, isin, arange
 from scipy.sparse import csr_matrix
 
 from pandapower.pypower.idx_brch import F_BUS, T_BUS, BR_STATUS, PF, PT, QF, QT
@@ -22,7 +22,7 @@ from pandapower.pypower.idx_gen import GEN_BUS, GEN_STATUS, PG, QG, QMIN, QMAX,
 EPS = finfo(float).eps
 
 
-def pfsoln(baseMVA, bus0, gen0, branch0, Ybus, Yf, Yt, V, ref, ref_gens, Ibus=None):
+def pfsoln(baseMVA, bus0, gen0, branch0, Ybus, Yf, Yt, V, ref, ref_gens, Ibus=None, limited_gens=None):
     """Updates bus, gen, branch data structures to match power flow soln.
 
     @author: Ray Zimmerman (PSERC Cornell)
@@ -38,17 +38,18 @@ def pfsoln(baseMVA, bus0, gen0, branch0, Ybus, Yf, Yt, V, ref, ref_gens, Ibus=No
     on = find(gen[:, GEN_STATUS] > 0)  # which generators are on?
     gbus = gen[on, GEN_BUS].astype(int)  # what buses are they at?
 
-    # xward: add ref buses that are not at the generators
-    xbus = setdiff1d(ref, gbus)
-
     # compute total injected bus powers
     Ibus = zeros(len(V)) if Ibus is None else Ibus
-    Sbus = V[gbus] * conj(Ybus[gbus, :] * V - Ibus[gbus])
-    Sbus_xw = V[xbus] * conj(Ybus[xbus, :] * V - Ibus[xbus])
+    Sbus = V * conj(Ybus * V - Ibus)
 
     _update_v(bus, V)
-    _update_q(baseMVA, bus, gen, gbus, Sbus, on)
-    _update_p(baseMVA, bus, gen, ref, gbus, on, r_[Sbus, Sbus_xw], ref_gens)
+    _update_q(baseMVA, bus, gen, gbus, Sbus[gbus], on)
+
+    if limited_gens is not None and len(limited_gens) > 0:
+        on = find((gen[:, GEN_STATUS] > 0) | isin(arange(len(gen)), limited_gens))
+        gbus = gen[on, GEN_BUS].astype(int)
+
+    _update_p(baseMVA, bus, gen, ref, gbus, Sbus, ref_gens)
 
     # ----- update/compute branch power flows -----
     out = find(branch[:, BR_STATUS] == 0)  # out-of-service branches
@@ -72,7 +73,25 @@ def _update_v(bus, V):
     bus[:, VA] = angle(V) * 180. / pi
 
 
-def _update_p(baseMVA, bus, gen, ref, gbus, on, Sbus, ref_gens):
+def _split_p_for_gens_at_same_bus(gen, p_bus, gens_at_bus, ref_gens):
+    if len(gens_at_bus) > 1:  # more than one generator at this ref bus
+        # subtract off what is generated by other gens at this bus
+        ext_grids = intersect1d(gens_at_bus, ref_gens)
+        pv_gens = setdiff1d(gens_at_bus, ext_grids)
+        p_ext_grids = p_bus - sum(gen[pv_gens, PG])
+        slack_weights = gen[ext_grids, SL_FAC]
+        sum_slack_weights = sum(slack_weights)
+        if sum_slack_weights > 0:
+            # distribute bus slack power according to the distributed slack contribution factors
+            gen[ext_grids, PG] = gen[ext_grids, PG] + (
+                        p_ext_grids - sum(gen[ext_grids, PG])) * slack_weights / sum_slack_weights
+        else:
+            gen[ext_grids, PG] = p_ext_grids / len(ext_grids)
+    else:
+        gen[gens_at_bus[0], PG] = p_bus
+
+
+def _update_p(baseMVA, bus, gen, ref, gbus, Sbus, ref_gens):
     # update Pg for slack bus(es)
     # inj P + local Pd
     for slack_bus in ref:
@@ -81,21 +100,9 @@ def _update_p(baseMVA, bus, gen, ref, gbus, on, Sbus, ref_gens):
         if len(gens_at_bus) == 0:
             bus[slack_bus, PD] = -Sbus[slack_bus].real * baseMVA  # negative because it is a PQ bus
             continue
-        p_bus = Sbus[gens_at_bus[0]].real * baseMVA + bus[slack_bus, PD]
-        if len(gens_at_bus) > 1:  # more than one generator at this ref bus
-            # subtract off what is generated by other gens at this bus
-            ext_grids = intersect1d(gens_at_bus, ref_gens)
-            pv_gens = setdiff1d(gens_at_bus, ext_grids)
-            p_ext_grids = p_bus - sum(gen[pv_gens, PG])
-            slack_weights = gen[ext_grids, SL_FAC]
-            sum_slack_weights = sum(slack_weights)
-            if sum_slack_weights > 0:
-                # distribute bus slack power according to the distributed slack contribution factors
-                gen[ext_grids, PG] = gen[ext_grids, PG] + (p_ext_grids - sum(gen[ext_grids, PG])) * slack_weights / sum_slack_weights
-            else:
-                gen[ext_grids, PG] = p_ext_grids / len(ext_grids)
-        else:
-            gen[on[gens_at_bus[0]], PG] = p_bus
+        # p_bus = Sbus[gens_at_bus[0]].real * baseMVA + bus[slack_bus, PD]
+        p_bus = Sbus[slack_bus].real * baseMVA + bus[slack_bus, PD]
+        _split_p_for_gens_at_same_bus(gen, p_bus, gens_at_bus, ref_gens)
 
 
 def _update_q(baseMVA, bus, gen, gbus, Sbus, on):


Test patch:
diff --git a/pandapower/test/loadflow/test_dist_slack.py b/pandapower/test/loadflow/test_dist_slack.py
index 70c7fed3..a67522c2 100644
--- a/pandapower/test/loadflow/test_dist_slack.py
+++ b/pandapower/test/loadflow/test_dist_slack.py
@@ -517,6 +517,26 @@ def test_dist_slack_with_enforce_q_lims():
     net.ext_grid['slack_weight'] = 1 / 3
     net.gen['slack_weight'] = 1 / 3
 
+    net.gen.at[0, 'max_q_mvar'] = 10
+    net.gen.at[1, 'max_q_mvar'] = 0.5
+    pp.runpp(net, distributed_slack=True, enforce_q_lims=True, numba=False)
+    assert net._options["distributed_slack"] and net._options["enforce_q_lims"]
+    assert np.allclose(net.res_gen.q_mvar, net.gen.max_q_mvar, rtol=0, atol=1e-6)
+
+    assert_results_correct(net, tol=1e-6)
+
+
+def test_dist_slack_with_enforce_q_lims_duplicate_gens():
+    net = pp.networks.case9()
+    pp.create_gen(net, net.ext_grid.bus.at[0], 1, slack=False, max_q_mvar=0.07)
+    pp.create_gen(net, net.ext_grid.bus.at[0], 1, slack=False, max_q_mvar=0.01)
+    pp.create_gen(net, net.ext_grid.bus.at[0], 2, slack=False, max_q_mvar=0.2)
+    pp.create_gen(net, net.ext_grid.bus.at[0], 2, slack=False, max_q_mvar=0.1)
+    pp.create_gen(net, net.ext_grid.bus.at[0], 0, slack=False, max_q_mvar=0.05)
+    pp.create_gen(net, net.ext_grid.bus.at[0], 0, slack=False, max_q_mvar=0.25)
+    net.gen['slack_weight'] = 1 / 8
+    net.ext_grid['slack_weight'] = 1 / 8
+
     net.gen.at[0, 'max_q_mvar'] = 10
     net.gen.at[1, 'max_q_mvar'] = 0.5
     pp.runpp(net, distributed_slack=True, enforce_q_lims=True)