Repository: e2nIEE/pandapower

Issue:
Multiple ext_grids in a single network
## Assumptions
- An ext_grid in pandapower is synonymous to an ideal voltage source in circuit analysis. 

- Since a circuit can have multiple voltage sources, and be perfectly well behaved - a pandapower model should be able to have multiple ext_grids. 

The key idea being that ext_grid is a way of specifying voltage at a bus - just like a load is a way of specifying power output at a bus.


## Basic example
```py
net = pp.create_empty_network()
pp.create_buses(net, nr_buses=2, vn_kv=1)
pp.create_line_from_parameters(net, from_bus=0, to_bus=1, length_km=1, 
                               r_ohm_per_km=1, x_ohm_per_km=1, 
                               c_nf_per_km=0, max_i_ka=1)
pp.create_ext_grid(net, bus=0, vm_pu=1)
pp.create_ext_grid(net, bus=1, vm_pu=1) # Adding this ext_grid breaks the powerflow
pp.runpp(net)
```



## Expected result
Since the voltage on both sides of the line is the same - there is no voltage drop, and hence also no power transfer/current.

---

Setting the voltage of the second ext_grid to 0.5 would mean a large amount of power is transferred from the first ext_grid to the second one through the line.

## Actual result
A value-error is raised:

<details><summary>Error logs</summary>

```py
ValueError                                Traceback (most recent call last)
<ipython-input-28-b5de297f7233> in <module>
      4 pp.create_ext_grid(net, bus=0, vm_pu=1)
      5 pp.create_ext_grid(net, bus=1, vm_pu=1)
----> 6 pp.runpp(net)

/opt/conda/lib/python3.8/site-packages/pandapower/run.py in runpp(net, algorithm, calculate_voltage_angles, init, max_iteration, tolerance_mva, trafo_model, trafo_loading, enforce_q_lims, check_connectivity, voltage_depend_loads, consider_line_temperature, run_control, **kwargs)
    223         _check_bus_index_and_print_warning_if_high(net)
    224         _check_gen_index_and_print_warning_if_high(net)
--> 225         _powerflow(net, **kwargs)
    226 
    227 

/opt/conda/lib/python3.8/site-packages/pandapower/powerflow.py in _powerflow(net, **kwargs)
     72 
     73     # ----- run the powerflow -----
---> 74     result = _run_pf_algorithm(ppci, net["_options"], **kwargs)
     75     # read the results (=ppci with results) to net
     76     _ppci_to_net(result, net)

/opt/conda/lib/python3.8/site-packages/pandapower/powerflow.py in _run_pf_algorithm(ppci, options, **kwargs)
    149             result = _run_bfswpf(ppci, options, **kwargs)[0]
    150         elif algorithm in ['nr', 'iwamoto_nr']:
--> 151             result = _run_newton_raphson_pf(ppci, options)
    152         elif algorithm in ['fdbx', 'fdxb', 'gs']:  # algorithms existing within pypower
    153             result = _runpf_pypower(ppci, options, **kwargs)[0]

/opt/conda/lib/python3.8/site-packages/pandapower/pf/run_newton_raphson_pf.py in _run_newton_raphson_pf(ppci, options)
     49         ppci, success, iterations, bus, gen, branch = _run_ac_pf_with_qlims_enforced(ppci, options)
     50     else:
---> 51         ppci, success, iterations = _run_ac_pf_without_qlims_enforced(ppci, options)
     52         # update data matrices with solution store in ppci
     53         bus, gen, branch = ppci_to_pfsoln(ppci, options)

/opt/conda/lib/python3.8/site-packages/pandapower/pf/run_newton_raphson_pf.py in _run_ac_pf_without_qlims_enforced(ppci, options)
    129     # run the newton power flow
    130     newton = newton_ls if lightsim2grid_available and options["lightsim2grid"] else newtonpf
--> 131     V, success, iterations, J, Vm_it, Va_it = newton(Ybus, Sbus, V0, pv, pq, ppci, options)
    132 
    133     # keep "internal" variables in  memory / net["_ppc"]["internal"] -> needed for recycle.

/opt/conda/lib/python3.8/site-packages/pandapower/pypower/newtonpf.py in newtonpf(Ybus, Sbus, V0, pv, pq, ppci, options)
     86     # evaluate F(x0)
     87     F = _evaluate_Fx(Ybus, V, Sbus, pv, pq)
---> 88     converged = _check_for_convergence(F, tol)
     89 
     90     Ybus = Ybus.tocsr()

/opt/conda/lib/python3.8/site-packages/pandapower/pypower/newtonpf.py in _check_for_convergence(F, tol)
    139 def _check_for_convergence(F, tol):
    140     # calc infinity norm
--> 141     return linalg.norm(F, Inf) < tol

<__array_function__ internals> in norm(*args, **kwargs)

/opt/conda/lib/python3.8/site-packages/numpy/linalg/linalg.py in norm(x, ord, axis, keepdims)
   2547     if len(axis) == 1:
   2548         if ord == Inf:
-> 2549             return abs(x).max(axis=axis, keepdims=keepdims)
   2550         elif ord == -Inf:
   2551             return abs(x).min(axis=axis, keepdims=keepdims)

/opt/conda/lib/python3.8/site-packages/numpy/core/_methods.py in _amax(a, axis, out, keepdims, initial, where)
     37 def _amax(a, axis=None, out=None, keepdims=False,
     38           initial=_NoValue, where=True):
---> 39     return umr_maximum(a, axis, None, out, keepdims, initial, where)
     40 
     41 def _amin(a, axis=None, out=None, keepdims=False,

ValueError: zero-size array to reduction operation maximum which has no identity
```

</details>

Discussion:

Patch:
diff --git a/pandapower/powerflow.py b/pandapower/powerflow.py
index 4309008e..cace0b8e 100644
--- a/pandapower/powerflow.py
+++ b/pandapower/powerflow.py
@@ -14,6 +14,7 @@ from pandapower.pf.run_bfswpf import _run_bfswpf
 from pandapower.pf.run_dc_pf import _run_dc_pf
 from pandapower.pf.run_newton_raphson_pf import _run_newton_raphson_pf
 from pandapower.pf.runpf_pypower import _runpf_pypower
+from pandapower.pypower.bustypes import bustypes
 from pandapower.pypower.idx_bus import VM
 from pandapower.pypower.makeYbus import makeYbus as makeYbus_pypower
 from pandapower.pypower.pfsoln import pfsoln as pfsoln_pypower
@@ -142,9 +143,10 @@ def _run_pf_algorithm(ppci, options, **kwargs):
     ac = options["ac"]
 
     if ac:
+        _, pv, pq = bustypes(ppci["bus"], ppci["gen"])
         # ----- run the powerflow -----
-        if ppci["branch"].shape[0] == 0:
-            result = _pf_without_branches(ppci, options)
+        if pq.shape[0] == 0 and pv.shape[0] == 0:
+            result = _bypass_pf_and_set_results(ppci, options)
         elif algorithm == 'bfsw':  # forward/backward sweep power flow algorithm
             result = _run_bfswpf(ppci, options, **kwargs)[0]
         elif algorithm in ['nr', 'iwamoto_nr']:
@@ -182,7 +184,7 @@ def _ppci_to_net(result, net):
     _clean_up(net)
 
 
-def _pf_without_branches(ppci, options):
+def _bypass_pf_and_set_results(ppci, options):
     Ybus, Yf, Yt = makeYbus_pypower(ppci["baseMVA"], ppci["bus"], ppci["branch"])
     baseMVA, bus, gen, branch, ref, _, pq, _, _, V0, ref_gens = _get_pf_variables_from_ppci(ppci)
     V = ppci["bus"][:, VM]


Test patch:
diff --git a/pandapower/test/loadflow/test_runpp.py b/pandapower/test/loadflow/test_runpp.py
index e0f3bbb5..a84920ba 100644
--- a/pandapower/test/loadflow/test_runpp.py
+++ b/pandapower/test/loadflow/test_runpp.py
@@ -1036,6 +1036,46 @@ def test_dc_with_ext_grid_at_one_bus():
     assert np.allclose(net.res_ext_grid.p_mw.values, [0, 0])
 
 
+def test_no_branches():
+    net = pp.create_empty_network()
+    pp.create_buses(net, 3, 110)
+    pp.create_ext_grid(net, 0)
+    pp.create_sgen(net, 1, 10)
+    pp.create_load(net, 2, 10)
+    pp.runpp(net)
+    assert net.res_ext_grid.p_mw.at[0] == 0.
+    assert net.res_ext_grid.q_mvar.at[0] == 0.
+    assert net.res_bus.vm_pu.at[0] == 1.
+    assert net.res_bus.va_degree.at[0] == 0.
+    assert np.all(pd.isnull(net.res_bus.loc[[1, 2], 'vm_pu']))
+
+
+def test_only_ref_buses():
+    net = pp.create_empty_network()
+    pp.create_buses(net, nr_buses=2, vn_kv=1)
+    pp.create_line_from_parameters(net, from_bus=0, to_bus=1, length_km=1,
+                                   r_ohm_per_km=1, x_ohm_per_km=1,
+                                   c_nf_per_km=0, max_i_ka=1)
+    pp.create_ext_grid(net, bus=0, vm_pu=1)
+    pp.create_ext_grid(net, bus=1, vm_pu=1)
+    pp.runpp(net)
+    assert np.all(net.res_bus.vm_pu == 1.)
+    assert np.all(net.res_bus.va_degree == 0.)
+    assert net.res_line.loading_percent.at[0] == 0.
+    assert np.all(net.res_ext_grid.p_mw == 0.)
+    assert np.all(net.res_ext_grid.q_mvar == 0.)
+
+    net.ext_grid.vm_pu.at[1] = 0.5
+    pp.runpp(net)
+    assert np.allclose(net.res_ext_grid.p_mw.values, np.array([0.25, -0.125]), rtol=0, atol=1e-12)
+    assert np.allclose(net.res_ext_grid.q_mvar.values, np.array([0.25, -0.125]), rtol=0, atol=1e-12)
+    assert abs(net.res_line.p_from_mw.at[0] - 0.25) < 1e-12
+    assert abs(net.res_line.q_from_mvar.at[0] - 0.25) < 1e-12
+    assert abs(net.res_line.p_to_mw.at[0] + 0.125) < 1e-12
+    assert abs(net.res_line.q_to_mvar.at[0] + 0.125) < 1e-12
+    assert abs(net.res_line.i_ka.at[0] - 0.20412415) < 1e-6
+
+
 def test_init_results_without_results():
     # should switch to "auto" mode and not fail
     net = example_multivoltage()