Allow meshed control transformers

docs/user_manual/calculations.md

@@ -628,7 +628,7 @@ Power flow calculations that take the behavior of these regulators into account
 
 We provide the control logic used for tap changing. For simplicity, we demonstrate the case where the regulator control side and the transformer tap side are at different sides.
 
-- Regulated transformers are ranked according to how close they are to {hoverxreftooltip}`sources <user_manual/components:source>` in terms of the amount of regulated transformers inbetween.
+- Regulated transformers are ranked according to how close they are to {hoverxreftooltip}`sources <user_manual/components:source>` in terms of the amount of regulated transformers inbetween. In the presence of meshed grids, transformers with conflicting ranks will be ranked the last.
   - Transformers are regulated in order according to their ranks.
 - Initialize all transformers to their starting tap position (see {hoverxreftooltip}`user_manual/calculations:Initialization and exploitation of regulated transformers`)
 - Find the optimal state using the following procedure

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

@@ -137,8 +137,8 @@ class IterationDiverge : public PowerGridError {
 
 class MaxIterationReached : public IterationDiverge {
   public:
-    MaxIterationReached(std::string_view msg = "") {
-        append_msg(std::format("Maximum number of iterations reached{}\n", msg));
+    MaxIterationReached(std::string const& msg = "") {
+        append_msg(std::format("Maximum number of iterations reached {}\n", msg));
     }
 };
 

power_grid_model_c/power_grid_model/include/power_grid_model/optimizer/tap_position_optimizer.hpp

@@ -41,7 +41,8 @@ using EdgeWeight = int64_t;
 using RankedTransformerGroups = std::vector<std::vector<Idx2D>>;
 
 constexpr auto infty = std::numeric_limits<Idx>::max();
-constexpr Idx2D unregulated_idx = {.group = -1, .pos = -1};
+constexpr auto last_rank = infty - 1;
+constexpr Idx2D unregulated_idx = {-1, -1};
 struct TrafoGraphVertex {
     bool is_source{};
 };
@@ -260,11 +261,12 @@ inline auto build_transformer_graph(State const& state) -> TransformerGraph {
     return trafo_graph;
 }
 
-inline void process_edges_dijkstra(Idx v, std::vector<EdgeWeight>& vertex_distances, TransformerGraph const& graph) {
+inline void process_edges_dijkstra(Idx search_start_vertex, std::vector<EdgeWeight>& vertex_distances,
+                                   TransformerGraph const& graph) {
     using TrafoGraphElement = std::pair<EdgeWeight, TrafoGraphIdx>;
     std::priority_queue<TrafoGraphElement, std::vector<TrafoGraphElement>, std::greater<>> pq;
-    vertex_distances[v] = 0;
-    pq.emplace(0, v);
+    vertex_distances[search_start_vertex] = 0;
+    pq.emplace(0, search_start_vertex);
 
     while (!pq.empty()) {
         auto [dist, u] = pq.top();
@@ -274,19 +276,14 @@ inline void process_edges_dijkstra(Idx v, std::vector<EdgeWeight>& vertex_distan
             continue;
         }
 
-        BGL_FORALL_EDGES(e, graph, TransformerGraph) {
+        BGL_FORALL_OUTEDGES(u, e, graph, TransformerGraph) {
             auto s = boost::source(e, graph);
             auto t = boost::target(e, graph);
             const EdgeWeight weight = graph[e].weight;
 
-            // We can not use BGL_FORALL_OUTEDGES here because we need information
-            // regardless of edge direction
-            if (u == s && vertex_distances[s] + weight < vertex_distances[t]) {
+            if (vertex_distances[s] + weight < vertex_distances[t]) {
                 vertex_distances[t] = vertex_distances[s] + weight;
                 pq.emplace(vertex_distances[t], t);
-            } else if (u == t && vertex_distances[t] + weight < vertex_distances[s]) {
-                vertex_distances[s] = vertex_distances[t] + weight;
-                pq.emplace(vertex_distances[s], s);
             }
         }
     }
@@ -328,11 +325,16 @@ inline auto get_edge_weights(TransformerGraph const& graph) -> TrafoGraphEdgePro
         // situations can happen.
         // The logic still holds in meshed grids, albeit operating a more complex graph.
         if (!is_unreachable(edge_src_rank) || !is_unreachable(edge_tgt_rank)) {
-            if (edge_src_rank != edge_tgt_rank - 1) {
-                throw AutomaticTapInputError("The control side of a transformer regulator should be relatively further "
-                                             "away from the source than the tap side.\n");
+            if ((edge_src_rank == infty) || (edge_tgt_rank == infty)) {
+                throw AutomaticTapInputError("The transformer is being controlled from non source side towards source "
+                                             "side.\n");
+            } else if (edge_src_rank != edge_tgt_rank - 1) {
+                // Control side is also controlled by a closer regulated transformer.
+                // Make this transformer have the lowest possible priority.
+                result.emplace_back(graph[e].regulated_idx, last_rank);
+            } else {
+                result.emplace_back(graph[e].regulated_idx, edge_tgt_rank);
             }
-            result.emplace_back(graph[e].regulated_idx, edge_tgt_rank);
         }
     }
 
@@ -663,7 +665,7 @@ template <symmetry_tag sym> struct NodeState {
 
 class RankIteration {
   public:
-    RankIteration(std::vector<IntS> iterations_per_rank, Idx rank_index)
+    RankIteration(std::vector<uint64_t> iterations_per_rank, Idx rank_index)
         : iterations_per_rank_{std::move(iterations_per_rank)}, rank_index_{rank_index} {}
 
     // Getters
@@ -692,7 +694,7 @@ class RankIteration {
     };
 
   private:
-    std::vector<IntS> iterations_per_rank_;
+    std::vector<uint64_t> iterations_per_rank_;
     Idx rank_index_{};
 };
 
@@ -1001,7 +1003,7 @@ class TapPositionOptimizerImpl<std::tuple<TransformerTypes...>, StateCalculator,
             strategy_ == OptimizerStrategy::global_maximum || strategy_ == OptimizerStrategy::local_maximum;
         bool tap_changed = true;
         Idx rank_index = 0;
-        RankIteration rank_iterator(std::vector<IntS>(regulator_order.size()), rank_index);
+        RankIteration rank_iterator(std::vector<uint64_t>(regulator_order.size()), rank_index);
 
         while (tap_changed) {
             tap_changed = false;
@@ -1021,8 +1023,10 @@ class TapPositionOptimizerImpl<std::tuple<TransformerTypes...>, StateCalculator,
             rank_index = rank_iterator.rank_index();
 
             if (tap_changed) {
-                if (static_cast<uint64_t>(iterations_per_rank[rank_index]) > 2 * max_tap_ranges_per_rank[rank_index]) {
-                    throw MaxIterationReached{"TapPositionOptimizer::iterate"};
+                if (iterations_per_rank[rank_index] > 2 * max_tap_ranges_per_rank[rank_index]) {
+                    throw MaxIterationReached{
+                        std::format("TapPositionOptimizer::iterate {} iterations reached: {}x2 iterations in rank {}",
+                                    iterations_per_rank[rank_index], max_tap_ranges_per_rank[rank_index], rank_index)};
                 }
                 update_state(update_data);
                 result = calculate_(state, method);

tests/cpp_unit_tests/test_tap_position_optimizer.cpp

@@ -422,7 +422,10 @@ TEST_CASE("Test Transformer ranking") {
 
         SUBCASE("Ranking complete the graph") {
             // The test grid 1 is not compatible with the updated logic for step up transformers
-            CHECK_THROWS_AS(pgm_tap::rank_transformers(state), AutomaticTapInputError);
+            pgm_tap::RankedTransformerGroups order = pgm_tap::rank_transformers(state);
+            pgm_tap::RankedTransformerGroups const ref_order{
+                {{Idx2D{3, 0}, Idx2D{3, 1}, Idx2D{4, 0}}, {Idx2D{3, 2}}, {Idx2D{3, 3}, Idx2D{3, 4}}}};
+            CHECK(order == ref_order);
         }
     }
 
@@ -507,6 +510,62 @@ TEST_CASE("Test Transformer ranking") {
              {Idx2D{.group = 3, .pos = 2}, Idx2D{.group = 3, .pos = 3}, Idx2D{.group = 3, .pos = 5}}}};
         CHECK(order == ref_order);
     }
+
+    SUBCASE("Full grid 3 - For Meshed grid with low priority ranks") {
+        // =====Test Grid=====
+        // ________[0]________
+        //  |              |
+        //  |              |
+        //  |              [1]
+        //  |              |
+        // _|______[2]_____|__
+        //          |
+        //         [3]
+        TestState state;
+        std::vector<NodeInput> nodes{{0, 10e3}, {1, 10e3}, {2, 10e3}, {3, 10e3}};
+        main_core::add_component<Node>(state, nodes.begin(), nodes.end(), 50.0);
+
+        std::vector<TransformerInput> transformers{get_transformer(11, 0, 1, BranchSide::to),
+                                                   get_transformer(12, 1, 2, BranchSide::from),
+                                                   get_transformer(13, 2, 3, BranchSide::from)};
+        main_core::add_component<Transformer>(state, transformers.begin(), transformers.end(), 50.0);
+
+        std::vector<LineInput> lines{get_line_input(21, 0, 2)};
+        main_core::add_component<Line>(state, lines.begin(), lines.end(), 50.0);
+
+        std::vector<SourceInput> sources{{31, 0, 1, 1.0, 0, nan, nan, nan}};
+        main_core::add_component<Source>(state, sources.begin(), sources.end(), 50.0);
+
+        std::vector<TransformerTapRegulatorInput> regulators{get_regulator(41, 11, ControlSide::to),
+                                                             get_regulator(42, 12, ControlSide::to),
+                                                             get_regulator(43, 13, ControlSide::to)};
+        main_core::add_component<TransformerTapRegulator>(state, regulators.begin(), regulators.end(), 50.0);
+
+        state.components.set_construction_complete();
+
+        pgm_tap::RankedTransformerGroups order = pgm_tap::rank_transformers(state);
+        pgm_tap::RankedTransformerGroups const ref_order{{{Idx2D{3, 0}}, {Idx2D{3, 2}}}, {{Idx2D{3, 1}}}};
+        CHECK(order == ref_order);
+    }
+
+    SUBCASE("Controlling from non source to source transformer") {
+        TestState state;
+        std::vector<NodeInput> nodes{{.id = 0, .u_rated = 150e3}, {.id = 1, .u_rated = 10e3}};
+        main_core::add_component<Node>(state, nodes.begin(), nodes.end(), 50.0);
+
+        std::vector<TransformerInput> transformers{get_transformer(2, 0, 1, BranchSide::from)};
+        main_core::add_component<Transformer>(state, transformers.begin(), transformers.end(), 50.0);
+
+        std::vector<SourceInput> sources{SourceInput{.id = 3, .node = 0, .status = IntS{1}, .u_ref = 1.0}};
+        main_core::add_component<Source>(state, sources.begin(), sources.end(), 50.0);
+
+        std::vector<TransformerTapRegulatorInput> regulators{get_regulator(4, 2, ControlSide::from)};
+        main_core::add_component<TransformerTapRegulator>(state, regulators.begin(), regulators.end(), 50.0);
+
+        state.components.set_construction_complete();
+
+        CHECK_THROWS_AS(pgm_tap::rank_transformers(state), AutomaticTapInputError);
+    }
 }
 
 namespace optimizer::tap_position_optimizer::test {
@@ -1603,7 +1662,7 @@ TEST_CASE("Test tap position optmizer I/O") {
 TEST_CASE("Test RankIterator") {
     std::vector<std::vector<IntS>> const regulator_order = {{0, 0, 0}, {0, 0, 0}};
     bool tap_changed{false};
-    std::vector<IntS> iterations_per_rank = {2, 4, 6};
+    std::vector<uint64_t> iterations_per_rank = {2, 4, 6};
     Idx rank_index{0};
     bool update{false};
     optimizer::tap_position_optimizer::RankIteration rank_iterator(iterations_per_rank, rank_index);

tests/unit/test_error_handling.py

@@ -323,10 +323,22 @@ def test_transformer_tap_regulator_control_side_not_closer_to_source():
     transformer_input["to_node"] = [1]
     transformer_input["from_status"] = [1]
     transformer_input["to_status"] = [1]
+    transformer_input["u1"] = [1e4]
+    transformer_input["u2"] = [4e2]
+    transformer_input["sn"] = [1e5]
+    transformer_input["uk"] = [0.1]
+    transformer_input["pk"] = [1e3]
+    transformer_input["i0"] = [1.0e-6]
+    transformer_input["p0"] = [0.1]
     transformer_input["winding_from"] = [1]
     transformer_input["winding_to"] = [1]
     transformer_input["clock"] = [0]
     transformer_input["tap_side"] = [1]
+    transformer_input["tap_pos"] = [0]
+    transformer_input["tap_nom"] = [0]
+    transformer_input["tap_min"] = [-1]
+    transformer_input["tap_max"] = [1]
+    transformer_input["tap_size"] = [100]
 
     source_input = initialize_array("input", "source", 1)
     source_input["id"] = [3]
@@ -338,6 +350,8 @@ def test_transformer_tap_regulator_control_side_not_closer_to_source():
     transformer_tap_regulator_input["id"] = [4]
     transformer_tap_regulator_input["regulated_object"] = [2]
     transformer_tap_regulator_input["status"] = [1]
+    transformer_tap_regulator_input["u_set"] = [10]
+    transformer_tap_regulator_input["u_band"] = [200]
     transformer_tap_regulator_input["control_side"] = [0]
 
     model = PowerGridModel(