Improve potentials propagation (#205)

Related to #200, I tried to get the maximum information from the voltage
sources to determine good potentials (good meaning potentials close to
the final ones).
And if it is not enough to determine all the potentials of the buses,
come back to something simple $[V_m, V_m e^{-i2\pi/3}, V_m e^{i2\pi/3},
0]$, as suggested by @alihamdan.

roseau/load_flow/network.py

@@ -1281,46 +1281,21 @@ def _reset_inputs(self) -> None:
     def _propagate_potentials(self) -> None:
         """Set the bus potentials that have not been initialized yet."""
         uninitialized = False
+        all_phases = set()
         for bus in self.buses.values():
             if not bus._initialized:
                 uninitialized = True
+                all_phases |= set(bus.phases)
 
         if uninitialized:
-            max_voltages = 0.0
-            voltage_source = None
-            potentials = None
-            for source in self.sources.values():
-                # if there are multiple voltage sources, start from the higher one
-                source_voltages = source._voltages
-                source_voltages_avg = np.average(np.abs(source_voltages))
-                if source_voltages_avg > max_voltages:
-                    max_voltages = source_voltages_avg
-                    voltage_source = source
-                    if "n" in source.phases:
-                        # Assume Vn = 0
-                        potentials = np.append(source_voltages, 0.0)
-                    elif len(source.phases) == 2:
-                        # Assume V1 + V2 = 0
-                        u = source_voltages[0]
-                        potentials = np.array([u / 2, -u / 2])
-                    else:
-                        assert len(source.phases) == 3
-                        # Assume Va + Vb + Vc = 0
-                        u_ab = source_voltages[0]
-                        u_bc = source_voltages[1]
-                        v_b = (u_bc - u_ab) / 3
-                        v_c = v_b - u_bc
-                        v_a = v_b + u_ab
-                        potentials = np.array([v_a, v_b, v_c, 0.0])
-
-            elements = [(voltage_source, potentials)]
+            potentials, starting_source = self._get_potentials(all_phases)
+            elements = [(starting_source, potentials)]
             visited = set()
             while elements:
                 element, potentials = elements.pop(-1)
                 visited.add(element)
                 if isinstance(element, Bus) and not element._initialized:
-                    bus_n = element._n
-                    element.potentials = potentials[0:bus_n]
+                    element.potentials = np.array([potentials[p] for p in element.phases], dtype=np.complex128)
                     element._initialized_by_the_user = False  # only used for serialization
                 for e in element._connected_elements:
                     if e not in visited and isinstance(e, (AbstractBranch, Bus)):
@@ -1329,13 +1304,52 @@ def _propagate_potentials(self) -> None:
                             phase_displacement = element.parameters.phase_displacement
                             if phase_displacement is None:
                                 phase_displacement = 0
-                            if element.parameters.type == "center" and "n" not in element.bus1.phases:
-                                # "n" is mandatory in the bus2 but not in the bus1
-                                potentials = np.append(potentials, 0.0)
-                            elements.append((e, potentials * k * np.exp(phase_displacement * -1j * np.pi / 6.0)))
+                            new_potentials = potentials.copy()
+                            for key, p in new_potentials.items():
+                                new_potentials[key] = p * k * np.exp(phase_displacement * -1j * np.pi / 6.0)
+                            elements.append((e, new_potentials))
                         else:
                             elements.append((e, potentials))
 
+    def _get_potentials(self, all_phases: set[str]) -> tuple[dict[str, complex], VoltageSource]:
+        """Compute initial potentials from the voltages sources of the network, get also the starting source"""
+        starting_source = None
+        potentials = {"n": 0.0}
+        # if there are multiple voltage sources, start from the higher one (the last one in the sorted below)
+        for source in sorted(self.sources.values(), key=lambda x: np.average(np.abs(x._voltages))):
+            source_voltages = source._voltages
+            starting_source = source
+            if "n" in source.phases:
+                # Assume Vn = 0
+                for phase, voltage in zip(source.phases[:-1], source_voltages, strict=True):
+                    potentials[phase] = voltage
+            elif len(source.phases) == 2:
+                # Assume V1 + V2 = 0
+                u = source_voltages[0]
+                potentials[source.phases[0]] = u / 2
+                potentials[source.phases[1]] = -u / 2
+            else:
+                assert source.phases == "abc"
+                # Assume Va + Vb + Vc = 0
+                u_ab = source_voltages[0]
+                u_bc = source_voltages[1]
+                v_b = (u_bc - u_ab) / 3
+                v_c = v_b - u_bc
+                v_a = v_b + u_ab
+                potentials["a"] = v_a
+                potentials["b"] = v_b
+                potentials["c"] = v_c
+
+        if len(potentials) < len(all_phases):
+            # We failed to determine all the potentials (the sources are strange), fallback to something simple
+            v = np.average(np.abs(starting_source._voltages))
+            potentials["a"] = v
+            potentials["b"] = v * np.exp(-2j * np.pi / 3)
+            potentials["c"] = v * np.exp(2j * np.pi / 3)
+            potentials["n"] = 0.0
+
+        return potentials, starting_source
+
     @staticmethod
     def _check_ref(elements: Iterable[Element]) -> None:
         """Check the number of potential references to avoid having a singular jacobian matrix."""

roseau/load_flow/tests/test_electrical_network.py

@@ -1623,6 +1623,50 @@ def _reset_inputs():
     assert not reset_inputs_called
 
 
+def test_propagate_potentials():
+    # Delta source
+    source_bus = Bus("source_bus", phases="abc")
+    _ = VoltageSource(id="source", bus=source_bus, voltages=20e3 * np.array([np.exp(1j * np.pi / 6), -1j, 0.0]))
+    _ = PotentialRef("pref", element=source_bus)
+    load_bus = Bus("load_bus", phases="abc")
+    _ = Switch("switch", bus1=source_bus, bus2=load_bus)
+
+    assert not load_bus._initialized
+    assert not source_bus._initialized
+    _ = ElectricalNetwork.from_element(source_bus)
+    assert load_bus._initialized
+    assert source_bus._initialized
+    un = 20e3 / np.sqrt(3)
+    expected_potentials = un * np.array([1, np.exp(-2j * np.pi / 3), np.exp(2j * np.pi / 3)])
+    assert np.allclose(load_bus.potentials.m, expected_potentials)
+    assert np.allclose(source_bus.potentials.m, expected_potentials)
+
+    # Multiple sources
+    source_bus = Bus("source_bus", phases="abcn")
+    _ = VoltageSource(id="VSa", bus=source_bus, voltages=[100], phases="an")
+    _ = VoltageSource(id="VSbc", bus=source_bus, voltages=[200, 300], phases="bcn")
+    _ = PotentialRef("pref", element=source_bus)
+    load_bus = Bus("load_bus", phases="abcn")
+    _ = Switch("switch", bus1=source_bus, bus2=load_bus)
+
+    assert not load_bus._initialized
+    _ = ElectricalNetwork.from_element(source_bus)
+    assert load_bus._initialized
+    assert np.allclose(load_bus.potentials.m, [100, 200, 300, 0])
+
+    # Do not define a source for all phases
+    source_bus = Bus("source_bus", phases="abcn")
+    _ = VoltageSource(id="VSa", bus=source_bus, voltages=[100], phases="an")
+    _ = PotentialRef("pref", element=source_bus)
+    load_bus = Bus("load_bus", phases="abcn")
+    _ = Switch("switch", bus1=source_bus, bus2=load_bus)
+
+    assert not load_bus._initialized
+    _ = ElectricalNetwork.from_element(source_bus)
+    assert load_bus._initialized
+    assert np.allclose(load_bus.potentials.m, 100 * np.array([1, np.exp(-2j * np.pi / 3), np.exp(2j * np.pi / 3), 0]))
+
+
 def test_short_circuits():
     vn = 400 / np.sqrt(3)
     voltages = [vn, vn * np.exp(-2 / 3 * np.pi * 1j), vn * np.exp(2 / 3 * np.pi * 1j)]