Merge branch 'main' into feature/document-high-level-design

Author: mgovers

.github/workflows/build-test-release.yml

@@ -208,7 +208,7 @@ jobs:
           xcode-version: latest-stable
 
       - name: Build wheels
-        uses: pypa/cibuildwheel@v2.23.2
+        uses: pypa/cibuildwheel@v2.23.3
         # GitHub Actions specific build parameters
         env:
           # pass GitHub runner info into Linux container

.github/workflows/ci.yml

@@ -61,3 +61,48 @@ jobs:
 
     steps:
       - run: echo "ci passed"
+
+  publish:
+    name: Publish to PyPI
+    runs-on: ubuntu-latest
+    permissions:
+      contents: write
+      id-token: write  # Required for Trusted Publishing
+    needs: build-test-release
+    if: (github.event_name == 'workflow_dispatch') || github.event_name == 'push'
+
+    steps:
+      - name: Download assets from GitHub release
+        uses: robinraju/release-downloader@v1
+        with:
+          repository: ${{ github.repository }}
+          # download the latest release
+          latest: true
+          # don't download pre-releases
+          preRelease: false
+          fileName: "*"
+          # don't download GitHub-generated source tar and zip files
+          tarBall: false
+          zipBall: false
+          # create a directory to store the downloaded assets
+          out-file-path: assets-to-publish
+          # don't extract downloaded files
+          extract: false
+
+      - name: List downloaded assets
+        run: ls -la assets-to-publish
+      
+      - name: Upload assets to PyPI
+        uses: pypa/gh-action-pypi-publish@release/v1
+        with:
+          # To test, use the TestPyPI:
+          # repository-url: https://test.pypi.org/legacy/
+          # You must also create an account and project on TestPyPI,
+          # as well as set the trusted-publisher in the project settings:
+          # https://docs.pypi.org/trusted-publishers/adding-a-publisher/
+          # To publish to the official PyPI repository, just keep
+          # repository-url commented out.
+          packages-dir: assets-to-publish
+          skip-existing: true
+          print-hash: true
+          verbose: true

.github/workflows/publish-pypi.yml

@@ -721,6 +721,117 @@ $$
    \end{eqnarray}
 $$
 
+### Generic Current Sensor
+
+```{warning}
+At the time of writing, this feature is still experimental and is not yet publicly available.
+```
+
+```{warning}
+At the time of writing, state estimation with current sensors is not supported by the Newton-Raphson calculation method.
+```
+
+* type name: `generic_current_sensor`
+
+`current_sensor` is an abstract class for symmetric and asymmetric current sensor and is derived from
+{hoverxreftooltip}`user_manual/components:sensor`. It measures the magnitude and angle of the current flow of a terminal.
+The terminal is connecting the from/to end of a `branch` (except `link`) and a `node`.
+
+```{note}
+Due to the high admittance of a `link` it is chosen that a current sensor cannot be coupled to a `link`, even though a link is a `branch`
+```
+
+##### Input
+
+| name                     | data type                                                                     | unit       | description                                                                                                                               | required |  update  |                     valid values                     |
+| ------------------------ | ----------------------------------------------------------------------------- | ---------- | ----------------------------------------------------------------------------------------------------------------------------------------- | :------: | :------: | :--------------------------------------------------: |
+| `measured_terminal_type` | {py:class}`MeasuredTerminalType <power_grid_model.enum.MeasuredTerminalType>` | -          | indicate the side of the `branch`                                                                                                         | &#10004; | &#10060; | the terminal type should match the `measured_object` |
+| `angle_measurement_type` | {py:class}`AngleMeasurementType <power_grid_model.enum.AngleMeasurementType>` | -          | indicate whether the measured angle is a global angle or a local angle; (see the [electric model](#local-angle-current-sensors) below)                                                                    | &#10004; | &#10060; |                                                      |
+| `i_sigma`                | `double`                                                                      | ampere (A) | standard deviation of the current (`i`) measurement error. Usually this is the absolute measurement error range divided by 3.             | &#10004; | &#10004; |                        `> 0`                         |
+| `i_angle_sigma`          | `double`                                                                      | rad        | standard deviation of the current (`i`) phase angle measurement error. Usually this is the absolute measurement error range divided by 3. | &#10004; | &#10004; |                        `> 0`                         |
+
+#### Current Sensor Concrete Types
+
+```{warning}
+At the time of writing, this feature is still experimental and is not yet publicly available.
+```
+
+```{warning}
+At the time of writing, state estimation with current sensors is not supported by the Newton-Raphson calculation method.
+```
+
+There are two concrete types of current sensor. They share similar attributes:
+the meaning of `RealValueInput` is different, as shown in the table below.
+
+| type name             | meaning of `RealValueInput` |
+| --------------------- | --------------------------- |
+| `sym_current_sensor`  | `double`                    |
+| `asym_current_sensor` | `double[3]`                 |
+
+##### Input
+
+| name               | data type        | unit       | description                               |              required              |  update  |
+| ------------------ | ---------------- | ---------- | ----------------------------------------- | :--------------------------------: | :------: |
+| `i_measured`       | `RealValueInput` | ampere (A) | measured current (`i`) magnitude          | &#10024; only for state estimation | &#10004; |
+| `i_angle_measured` | `RealValueInput` | rad        | measured phase angle of the current (`i`; see the [electric model](#local-angle-current-sensors) below) | &#10024; only for state estimation | &#10004; |
+
+See the documentation on [state estimation calculation methods](calculations.md#state-estimation-algorithms) for details per method on how the variances are taken into account for both the global and local angle measurement types and for the individual phases.
+
+##### Steady state output
+
+```{note}
+A sensor only has output for state estimation. For other calculation types, sensor output is undefined.
+```
+
+| name               | data type         | unit       | description                                                                                 |
+| ------------------ | ----------------- | ---------- | ------------------------------------------------------------------------------------------- |
+| `i_residual`       | `RealValueOutput` | ampere (A) | residual value between measured current (`i`) and calculated current (`i`)                  |
+| `i_angle_residual` | `RealValueOutput` | rad        | residual value between measured phase angle and calculated phase angle of the current (`i`) |
+
+#### Electric Model
+
+`Generic Current Sensor` is modeled by following equations:
+
+##### Global angle current sensors
+
+Current sensors with `angle_measurement_type` equal to `AngleMeasurementType.global_angle` measure the phase of the current relative to
+some reference angle that is the same across the entire grid. This reference angle must be the same one as for [voltage phasor measurements](#generic-voltage-sensor). Because the reference point may be ambiguous in the case of current sensor measurements, the power-grid-model imposes the following requirement:
+
+```{note}
+Global angle current measurements require at least one voltage angle measurement to make sense.
+```
+
+As a sign convention, the angle is the phase shift of the current relative to the reference angle, i.e.,
+
+$$
+\underline{I} = \text{i_measured} \cdot e^{j \text{i_angle_measured}} \text{ .}
+$$
+
+##### Local angle current sensors
+
+Current sensors with `angle_measurement_type` equal to `AngleMeasurementType.local_angle` measure the phase shift between
+the voltage and the current phasor, i.e., $\text{i_angle_measured} = \text{voltage_phase} - \text{current_phase}$.
+As a result, the global current phasor depends on the local voltage phase offset and is obtained using the following formula.
+
+$$
+\underline{I} = \underline{I}_{\text{local}}^{*} \frac{\underline{U}}{|\underline{U}|} = \text{i_measured} \cdot e^{\mathrm{j} \left(\theta_{U} - \text{i_angle_measured}\right)}
+$$
+
+```{note}
+As a result, the local angle current sensors have a different sign convention from the global angle current sensors.
+```
+
+##### Residuals
+
+$$
+   \begin{eqnarray}
+        & i_{\text{residual}} = i_{\text{measured}} - i_{\text{state}} && \\
+        & i_{\text{angle},\text{residual}} = i_{\text{angle},\text{measured}} - i_{\text{angle},\text{state}} \pmod 2 \pi
+   \end{eqnarray}
+$$
+
+The $\pmod 2\pi$ is handled such that $-\pi \lt i_{\text{angle},\text{residual}} \leq \pi$.
+
 ## Fault
 
 * type name: `fault`

docs/user_manual/calculations.md

@@ -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));
     }
 };
 

docs/user_manual/components.md

@@ -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);