Branch into itself: core implementation

power_grid_model_c/power_grid_model/include/power_grid_model/common/exception.hpp

@@ -73,21 +73,6 @@ class ConflictVoltage : public PowerGridError {
               id1, u1, id2, u2)} {}
 };
 
-class InvalidBranch : public PowerGridError {
-  public:
-    InvalidBranch(ID branch_id, ID node_id)
-        : PowerGridError{std::format("Branch {} has the same from- and to-node {},\n This is not allowed!\n", branch_id,
-                                     node_id)} {}
-};
-
-class InvalidBranch3 : public PowerGridError {
-  public:
-    InvalidBranch3(ID branch3_id, ID node_1_id, ID node_2_id, ID node_3_id)
-        : PowerGridError{std::format(
-              "Branch3 {} is connected to the same node at least twice. Node 1/2/3: {}/{}/{},\n This is not allowed!\n",
-              branch3_id, node_1_id, node_2_id, node_3_id)} {}
-};
-
 class InvalidTransformerClock : public PowerGridError {
   public:
     InvalidTransformerClock(ID id, IntS clock)

power_grid_model_c/power_grid_model/include/power_grid_model/common/three_phase_tensor.hpp

@@ -187,6 +187,16 @@ inline ComplexValue<asymmetric_t> piecewise_complex_value(T<Derived> const& x) {
 // piecewise factory construction for complex vector
 inline ComplexValue<asymmetric_t> piecewise_complex_value(ComplexValue<asymmetric_t> const& x) { return x; }
 
+// basic operations
+using Eigen::arg;
+using Eigen::conj;
+using Eigen::imag;
+using Eigen::real;
+using std::arg;
+using std::conj;
+using std::imag;
+using std::real;
+
 // abs
 template <std::floating_point Float> inline Float cabs(Float x) { return std::abs(x); }
 template <std::floating_point Float> inline Float abs2(Float x) { return x * x; }

power_grid_model_c/power_grid_model/include/power_grid_model/component/branch.hpp

@@ -41,11 +41,7 @@ class Branch : public Base {
           from_node_{branch_input.from_node},
           to_node_{branch_input.to_node},
           from_status_{static_cast<bool>(branch_input.from_status)},
-          to_status_{static_cast<bool>(branch_input.to_status)} {
-        if (from_node_ == to_node_) {
-            throw InvalidBranch{id(), from_node_};
-        }
-    }
+          to_status_{static_cast<bool>(branch_input.to_status)} {}
 
     // getter
     constexpr ID from_node() const { return from_node_; }

power_grid_model_c/power_grid_model/include/power_grid_model/component/branch3.hpp

@@ -41,11 +41,7 @@ class Branch3 : public Base {
           node_3_{branch3_input.node_3},
           status_1_{static_cast<bool>(branch3_input.status_1)},
           status_2_{static_cast<bool>(branch3_input.status_2)},
-          status_3_{static_cast<bool>(branch3_input.status_3)} {
-        if (node_1_ == node_2_ || node_1_ == node_3_ || node_2_ == node_3_) {
-            throw InvalidBranch3{id(), node_1_, node_2_, node_3_};
-        }
-    }
+          status_3_{static_cast<bool>(branch3_input.status_3)} {}
 
     // getter
     constexpr ID node_1() const { return node_1_; }

power_grid_model_c/power_grid_model/include/power_grid_model/topology.hpp

@@ -4,6 +4,7 @@
 
 #pragma once
 
+#include "common/typing.hpp"
 #include "sparse_ordering.hpp"
 
 #include "calculation_parameters.hpp"
@@ -129,7 +130,7 @@ class Topology {
           phase_shift_(comp_topo_.n_node_total(), 0.0),
           predecessors_(
               boost::counting_iterator<GraphIdx>{0}, // Predecessors is initialized as 0, 1, 2, ..., n_node_total() - 1
-              boost::counting_iterator<GraphIdx>{(GraphIdx)comp_topo_.n_node_total()}),
+              boost::counting_iterator<GraphIdx>{narrow_cast<GraphIdx>(comp_topo_.n_node_total())}),
           node_status_(comp_topo_.n_node_total(), -1) {}
 
     // build topology
@@ -202,12 +203,12 @@ class Topology {
             auto const [i, j] = branch_node_idx;
             auto const [i_status, j_status] = branch_connected;
             // node_i - node_j
-            if (i_status != 0 && j_status != 0) {
+            if (i_status != 0 && j_status != 0 && i != j) {
                 // node_j - node_i
-                edges.emplace_back((GraphIdx)i, (GraphIdx)j);
+                edges.emplace_back(narrow_cast<GraphIdx>(i), narrow_cast<GraphIdx>(j));
                 edge_props.push_back({-phase_shift});
                 // node_i - node_j
-                edges.emplace_back((GraphIdx)j, (GraphIdx)i);
+                edges.emplace_back(narrow_cast<GraphIdx>(j), narrow_cast<GraphIdx>(i));
                 edge_props.push_back({phase_shift});
             }
         }
@@ -219,17 +220,17 @@ class Topology {
             for (Idx m = 0; m != 3; ++m) {
                 if (i_status[m] != 0) {
                     // node_internal - node_i
-                    edges.emplace_back((GraphIdx)i[m], (GraphIdx)j_internal);
+                    edges.emplace_back(narrow_cast<GraphIdx>(i[m]), narrow_cast<GraphIdx>(j_internal));
                     edge_props.push_back({-phase_shift[m]});
                     // node_i - node_internal
-                    edges.emplace_back((GraphIdx)j_internal, (GraphIdx)i[m]);
+                    edges.emplace_back(narrow_cast<GraphIdx>(j_internal), narrow_cast<GraphIdx>(i[m]));
                     edge_props.push_back({phase_shift[m]});
                 }
             }
         }
         // build graph
         global_graph_ = GlobalGraph{boost::edges_are_unsorted_multi_pass, edges.cbegin(), edges.cend(),
-                                    edge_props.cbegin(), (GraphIdx)comp_topo_.n_node_total()};
+                                    edge_props.cbegin(), narrow_cast<GraphIdx>(comp_topo_.n_node_total())};
         for_all_vertices(global_graph_, [this](boost::graph_traits<GlobalGraph>::vertex_descriptor const& v) {
             global_graph_[v].color = boost::default_color_type::white_color;
         });
@@ -256,7 +257,7 @@ class Topology {
             std::vector<std::pair<GraphIdx, GraphIdx>> back_edges;
             // start dfs search
             boost::depth_first_visit(
-                global_graph_, (GraphIdx)source_node,
+                global_graph_, narrow_cast<GraphIdx>(source_node),
                 GlobalDFSVisitor{math_solver_idx, comp_coup_.node, phase_shift_, dfs_node, predecessors_, back_edges},
                 boost::get(&GlobalVertex::color, global_graph_));
 

src/power_grid_model/_core/error_handling.py

@@ -21,8 +21,6 @@
     IDNotFound,
     IDWrongType,
     InvalidArguments,
-    InvalidBranch,
-    InvalidBranch3,
     InvalidCalculationMethod,
     InvalidMeasuredObject,
     InvalidRegulatedObject,
@@ -63,11 +61,6 @@ class _PgmCErrorCode(IntEnum):
     r"Conflicting voltage for line (-?\d+)\n voltage at from node (-?\d+) is (.*)\n"
     r" voltage at to node (-?\d+) is (.*)\n"
 )
-_INVALID_BRANCH_RE = re.compile(r"Branch (-?\d+) has the same from- and to-node (-?\d+),\n This is not allowed!\n")
-_INVALID_BRANCH3_RE = re.compile(
-    r"Branch3 (-?\d+) is connected to the same node at least twice. Node 1\/2\/3: (-?\d+)\/(-?\d+)\/(-?\d+),\n"
-    r" This is not allowed!\n"
-)
 _INVALID_TRANSFORMER_CLOCK_RE = re.compile(r"Invalid clock for transformer (-?\d+), clock (-?\d+)\n")
 _SPARSE_MATRIX_ERROR_RE = re.compile(r"Sparse matrix error")  # multiple different flavors
 _NOT_OBSERVABLE_ERROR_RE = re.compile(r"Not enough measurements available for state estimation.\n")
@@ -107,8 +100,6 @@ class _PgmCErrorCode(IntEnum):
     _MISSING_CASE_FOR_ENUM_RE: MissingCaseForEnumError,
     _INVALID_ARGUMENTS_RE: InvalidArguments,
     _CONFLICT_VOLTAGE_RE: ConflictVoltage,
-    _INVALID_BRANCH_RE: InvalidBranch,
-    _INVALID_BRANCH3_RE: InvalidBranch3,
     _INVALID_TRANSFORMER_CLOCK_RE: InvalidTransformerClock,
     _SPARSE_MATRIX_ERROR_RE: SparseMatrixError,
     _NOT_OBSERVABLE_ERROR_RE: NotObservableError,

src/power_grid_model/_core/errors.py

@@ -6,8 +6,12 @@
 This file contains error classes for library-internal use.
 """
 
+import warnings
+
 import numpy as np
 
+_DEPRECATED_ERROR_MSG = "This error type is deprecated and may be reduced in a future release."
+
 
 class PowerGridError(RuntimeError):
     """Generic power grid error."""
@@ -37,11 +41,23 @@ class ConflictVoltage(PowerGridError):
 
 
 class InvalidBranch(PowerGridError):
-    """A branch is invalid."""
+    """A branch is invalid.
+
+    [DEPRECATED] This error is no longer relevant and may be reduced in a future release."""
+
+    def __init__(self, *args, **kwargs):
+        warnings.warn(_DEPRECATED_ERROR_MSG, DeprecationWarning)
+        super().__init__(*args, **kwargs)
 
 
 class InvalidBranch3(PowerGridError):
-    """A branch3 is invalid."""
+    """A branch3 is invalid.
+
+    [DEPRECATED] This error is no longer relevant and may be reduced in a future release."""
+
+    def __init__(self, *args, **kwargs):
+        warnings.warn(_DEPRECATED_ERROR_MSG, DeprecationWarning)
+        super().__init__(*args, **kwargs)
 
 
 class InvalidTransformerClock(PowerGridError):

tests/cpp_unit_tests/test_exceptions.cpp

@@ -105,19 +105,6 @@ TEST_CASE("Exceptions") {
                               .what()} ==
               "Conflicting voltage for line 0\n voltage at from node 0 is inf\n voltage at to node 0 is -inf\n");
     }
-    SUBCASE("InvalidBranch") {
-        CHECK(std::string{InvalidBranch{ID{0}, ID{1}}.what()} ==
-              "Branch 0 has the same from- and to-node 1,\n This is not allowed!\n");
-        CHECK(std::string{InvalidBranch{na_IntID, na_IntID}.what()} ==
-              "Branch -2147483648 has the same from- and to-node -2147483648,\n This is not allowed!\n");
-    }
-    SUBCASE("InvalidBranch3") {
-        CHECK(std::string{InvalidBranch3{ID{0}, ID{4}, ID{5}, ID{6}}.what()} ==
-              "Branch3 0 is connected to the same node at least twice. Node 1/2/3: 4/5/6,\n This is not allowed!\n");
-        CHECK(std::string{InvalidBranch3{na_IntID, na_IntID, na_IntID, na_IntID}.what()} ==
-              "Branch3 -2147483648 is connected to the same node at least twice. Node 1/2/3: "
-              "-2147483648/-2147483648/-2147483648,\n This is not allowed!\n");
-    }
     SUBCASE("InvalidTransformerClock") {
         CHECK(std::string{InvalidTransformerClock{ID{0}, IntS{1}}.what()} ==
               "Invalid clock for transformer 0, clock 1\n");

tests/cpp_unit_tests/test_line.cpp

@@ -61,6 +61,10 @@ TEST_CASE("Test line") {
     ComplexValue<asymmetric_t> const uaf{1.0};
     ComplexValue<asymmetric_t> const uat{0.9};
 
+    // regular result
+    BranchSolverOutput<symmetric_t> const branch_solver_output{
+        .s_f = 1.0 - 1.5i, .s_t = 1.5 - 1.5i, .i_f = 1.0 - 2.0i, .i_t = 2.0 - 1.0i};
+
     // Short circuit results
     DoubleComplex const if_sc{1.0, 1.0};
     DoubleComplex const it_sc{2.0, 2.0 * sqrt3};
@@ -69,7 +73,7 @@ TEST_CASE("Test line") {
 
     CHECK(line.math_model_type() == ComponentType::branch);
 
-    SUBCASE("Voltge error") { CHECK_THROWS_AS(Line(input, 50.0, 10.0e3, 50.0e3), ConflictVoltage); }
+    SUBCASE("Voltage error") { CHECK_THROWS_AS(Line(input, 50.0, 10.0e3, 50.0e3), ConflictVoltage); }
 
     SUBCASE("General") {
         CHECK(branch.from_node() == 2);
@@ -140,11 +144,6 @@ TEST_CASE("Test line") {
     }
 
     SUBCASE("Symmetric results with direct power and current output") {
-        BranchSolverOutput<symmetric_t> branch_solver_output{};
-        branch_solver_output.i_f = 1.0 - 2.0i;
-        branch_solver_output.i_t = 2.0 - 1.0i;
-        branch_solver_output.s_f = 1.0 - 1.5i;
-        branch_solver_output.s_t = 1.5 - 1.5i;
         BranchOutput<symmetric_t> const output = branch.get_output<symmetric_t>(branch_solver_output);
         CHECK(output.id == 1);
         CHECK(output.energized);
@@ -246,6 +245,56 @@ TEST_CASE("Test line") {
         CHECK(inv.from_status == expected.from_status);
         CHECK(inv.to_status == expected.to_status);
     }
+
+    SUBCASE("Lines into itself") {
+        auto line_into_itself_input = input;
+        line_into_itself_input.to_node = 2;
+        Line const line_into_itself{line_into_itself_input, 50.0, 10.0e3, 10.0e3};
+        Branch const& branch_into_itself = line_into_itself;
+
+        CHECK(branch_into_itself.from_node() == branch.from_node());
+        CHECK(branch_into_itself.to_node() == branch_into_itself.from_node());
+
+        CHECK(branch_into_itself.from_status() == branch.from_status());
+        CHECK(branch_into_itself.to_status() == branch.to_status());
+        CHECK(branch_into_itself.branch_status() == branch.branch_status());
+        CHECK(branch_into_itself.status(BranchSide::from) == branch.status(BranchSide::from));
+        CHECK(branch_into_itself.status(BranchSide::to) == branch.status(BranchSide::to));
+        CHECK(branch_into_itself.base_i_from() == branch.base_i_from());
+        CHECK(branch_into_itself.base_i_to() == branch.base_i_to());
+        CHECK(branch_into_itself.phase_shift() == branch.phase_shift());
+        CHECK(branch_into_itself.is_param_mutable() == branch.is_param_mutable());
+
+        SUBCASE("Symmetric parameters") {
+            auto const params = branch_into_itself.calc_param<symmetric_t>();
+            auto const ref_params = branch.calc_param<symmetric_t>();
+            CHECK(params.yff() == ref_params.yff());
+            CHECK(params.ytt() == ref_params.ytt());
+            CHECK(params.ytf() == ref_params.ytf());
+            CHECK(params.yft() == ref_params.yft());
+        }
+        SUBCASE("Asymmetric parameters") {
+            auto const params = branch_into_itself.calc_param<asymmetric_t>();
+            auto const ref_params = branch.calc_param<asymmetric_t>();
+            CHECK((cabs(params.yff() - ref_params.yff()) < numerical_tolerance).all());
+            CHECK((cabs(params.ytt() - ref_params.ytt()) < numerical_tolerance).all());
+            CHECK((cabs(params.ytf() - ref_params.ytf()) < numerical_tolerance).all());
+            CHECK((cabs(params.yft() - ref_params.yft()) < numerical_tolerance).all());
+        }
+        SUBCASE("Sym output") {
+            auto const branch_into_itself_solver_output =
+                BranchSolverOutput<symmetric_t>{.s_f = branch_solver_output.s_f,
+                                                .s_t = branch_solver_output.s_f, // same node, so should be s_f
+                                                .i_f = branch_solver_output.i_f,
+                                                .i_t = branch_solver_output.i_f};
+            auto const branch_into_itself_output =
+                branch_into_itself.get_output<symmetric_t>(branch_into_itself_solver_output);
+            CHECK(branch_into_itself_output.s_from == branch.get_output<symmetric_t>(branch_solver_output).s_from);
+            CHECK(branch_into_itself_output.s_to == branch_into_itself_output.s_from);
+            CHECK(branch_into_itself_output.i_from == branch.get_output<symmetric_t>(branch_solver_output).i_from);
+            CHECK(branch_into_itself_output.i_to == branch_into_itself_output.i_from);
+        }
+    }
 }
 
 } // namespace power_grid_model

tests/cpp_unit_tests/test_link.cpp

@@ -11,13 +11,11 @@
 #include <power_grid_model/calculation_parameters.hpp>
 #include <power_grid_model/common/common.hpp>
 #include <power_grid_model/common/enum.hpp>
-#include <power_grid_model/common/exception.hpp>
 #include <power_grid_model/common/three_phase_tensor.hpp>
 #include <power_grid_model/component/base.hpp>
 
 #include <doctest/doctest.h>
 
-#include <cmath>
 #include <complex>
 
 namespace power_grid_model {
@@ -50,9 +48,9 @@ TEST_CASE("Test link") {
         CHECK(branch.phase_shift() == 0.0);
     }
 
-    SUBCASE("Invalid branch") {
+    SUBCASE("Branch into itself is supported") {
         input.to_node = 2;
-        CHECK_THROWS_AS(Link(input, 10e3, 50e3), InvalidBranch);
+        CHECK_NOTHROW(Link(input, 10e3, 50e3));
     }
 
     SUBCASE("Symmetric parameters") {