Merge pull request #858 from PowerGridModel/feature/full-observabilit…

…y-radial

Run full observability check for radial grid without phasor measurement

docs/user_manual/calculations.md

@@ -57,10 +57,10 @@ Output:
 - Power flow through branches
 - Deviation between measurement values and estimated state
 
-In order to perform a state estimation, the system should be observable. If the system is not observable, the calculation will 
-raise a sparse matrix error; the matrix reprensentation of the system is too sparse to solve and particular so when it is
-singular. In short, meeting the requirement of observability indicates that the system is either an overdetermined system (when the number of measurements is larger than the number of
-unknowns) or a balanced system (when the number of measurements is equal to the number of unknowns). For each node, there are two unknowns, `u` and `u_angle`, so the following conditions should be met:
+In order to perform a state estimation, the system should be observable. If the system is not observable, the calculation will raise either a `NotObservableError` or 
+a `SparseMatrixError`.
+In short, meeting the requirement of observability indicates that the system is either an overdetermined system (when the number of measurements is larger than the number of unknowns.
+For each node, there are two unknowns, `u` and `u_angle`. Due to the relative nature of `u_angle` (relevant only in systems with at least two nodes), in total the following conditions should be met:
 
 $$
     \begin{eqnarray}
@@ -72,7 +72,7 @@ Where
 
 $$
     \begin{eqnarray}
-        n_{unknowns}    & = & 2 & \cdot & n_{nodes}
+        n_{unknowns}    & = & 2 & \cdot & n_{nodes} - 1
     \end{eqnarray}
 $$
 
@@ -96,7 +96,7 @@ In observable systems this helps better outputting correct results. On the other
 
 ##### Necessary observability condition
 
-Based on the requirements of observability mentioned above, user needs to satisfy  at least the following conditions for state estimation calculation in power-grid-model.
+Based on the requirements of observability mentioned above, users need to satisfy at least the following conditions for state estimation calculation in `power-grid-model`.
 
 - `n_voltage_sensor >= 1`
 - If no voltage phasor sensors are available, then the following conditions should be satisfied:  `n_unique_power_sensor >= n_bus - 1`. Otherwise: `n_unique_power_sensor + n_voltage_sensor_with_phasor >= n_bus`
@@ -108,6 +108,16 @@ Based on the requirements of observability mentioned above, user needs to satisf
 - Any sensor on a `Branch` for all branches: Parallel branches with either side of measurements count as one.
 - All `Branch3` sensors.
 
+##### Sufficient observability condition
+
+The condition check above only checks the necessary condition for observability. When the measurements are not independent enough, the system may still be unobservable even if the necessary condition is met.
+It is rather complicated to do a full sufficient and necessary observability check in generic cases. However, `power-grid-model` performs the sufficient condition check when the following conditions are met:
+
+1. The system is a radial network.
+2. The system does not have voltage phasor measurements.
+
+In this case, the validation of the independent measurements is rather straightforward. If the system is not observable, the calculation will raise a `NotObservableError` instead of `SparseMatrixError`.
+
 #### Short circuit calculations
 
 Short circuit calculation is carried out to analyze the worst case scenario when a fault has occurred.

power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp

@@ -141,6 +141,7 @@ using Branch3Idx = std::array<Idx, 3>;
 
 struct MathModelTopology {
     Idx slack_bus{};
+    bool is_radial{};
     std::vector<double> phase_shift;
     std::vector<BranchIdx> branch_bus_idx;
     std::vector<BranchIdx> fill_in;

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

@@ -106,7 +106,9 @@ class SparseMatrixError : public PowerGridError {
     SparseMatrixError() {
         append_msg("Sparse matrix error, possibly singular matrix!\n" +
                    std::string("If you get this error from state estimation, ") +
-                   "it usually means the system is not fully observable, i.e. not enough measurements.");
+                   "it might mean the system is not fully observable, i.e. not enough measurements.\n" +
+                   "It might also mean that you are running into a corner case where PGM cannot resolve yet." +
+                   "See https://github.com/PowerGridModel/power-grid-model/issues/853.");
     }
 };
 

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp

@@ -95,7 +95,7 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
         // preprocess measured value
         sub_timer = Timer(calculation_info, 2221, "Pre-process measured value");
         MeasuredValues<sym> const measured_values{y_bus.shared_topology(), input};
-        necessary_observability_check(measured_values, y_bus.shared_topology());
+        necessary_observability_check(measured_values, y_bus.math_topology(), y_bus.y_bus_structure());
 
         // prepare matrix, including pre-factorization
         sub_timer = Timer(calculation_info, 2222, "Prepare matrix, including pre-factorization");

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp

@@ -159,7 +159,7 @@ template <symmetry_tag sym_type> class NewtonRaphsonSESolver {
         // preprocess measured value
         sub_timer = Timer(calculation_info, 2221, "Pre-process measured value");
         MeasuredValues<sym> const measured_values{y_bus.shared_topology(), input};
-        necessary_observability_check(measured_values, y_bus.shared_topology());
+        necessary_observability_check(measured_values, y_bus.math_topology(), y_bus.y_bus_structure());
 
         // initialize voltage with initial angle
         sub_timer = Timer(calculation_info, 2223, "Initialize voltages");

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp

@@ -5,44 +5,14 @@
 #pragma once
 
 #include "measured_values.hpp"
+#include "y_bus.hpp"
 
 #include "../common/exception.hpp"
 
 namespace power_grid_model::math_solver {
 
 namespace detail {
 
-template <symmetry_tag sym>
-Idx count_branch_sensors(const std::vector<BranchIdx>& branch_bus_idx, const Idx n_bus,
-                         const MeasuredValues<sym>& measured_values) {
-    Idx n_branch_sensor{};
-    std::vector<bool> measured_nodes(n_bus, false);
-    for (Idx branch = 0; branch != static_cast<Idx>(branch_bus_idx.size()); ++branch) {
-        auto const& [node_from, node_to] = branch_bus_idx[branch];
-        if (node_from == -1 || node_to == -1) {
-            continue;
-        }
-        if ((measured_values.has_branch_from(branch) || measured_values.has_branch_to(branch)) &&
-            !(measured_nodes[node_from] && measured_nodes[node_to])) {
-            n_branch_sensor++;
-            measured_nodes[node_from] = true;
-            measured_nodes[node_to] = true;
-        }
-    }
-    return n_branch_sensor;
-}
-
-template <symmetry_tag sym>
-Idx count_bus_injection_sensors(const Idx n_bus, const MeasuredValues<sym>& measured_values) {
-    Idx n_injection_sensor{};
-    for (Idx bus = 0; bus != n_bus; ++bus) {
-        if (measured_values.has_bus_injection(bus)) {
-            n_injection_sensor++;
-        }
-    }
-    return n_injection_sensor;
-}
-
 template <symmetry_tag sym>
 std::tuple<Idx, Idx> count_voltage_sensors(const Idx n_bus, const MeasuredValues<sym>& measured_values) {
     Idx n_voltage_sensor{};
@@ -58,29 +28,114 @@ std::tuple<Idx, Idx> count_voltage_sensors(const Idx n_bus, const MeasuredValues
     return std::make_tuple(n_voltage_sensor, n_voltage_phasor_sensor);
 }
 
-} // namespace detail
+// count flow sensors into integer matrix with ybus structure
+// lower triangle part is always zero
+// for diagonal part, it will be one if there is bus injection
+// for upper triangle part, it will be one if there is branch flow sensor and the branch is fully connected
 template <symmetry_tag sym>
-inline void necessary_observability_check(MeasuredValues<sym> const& measured_values,
-                                          std::shared_ptr<MathModelTopology const> const& topo) {
+std::vector<int8_t> count_flow_sensors(MeasuredValues<sym> const& measured_values, MathModelTopology const& topo,
+                                       YBusStructure const& y_bus_structure) {
+    Idx const n_bus{topo.n_bus()};
+    std::vector<int8_t> flow_sensors(y_bus_structure.row_indptr.back(), 0); // initialize all to zero
+    for (Idx row = 0; row != n_bus; ++row) {
+        // lower triangle is ignored and kept as zero
+        // diagonal for bus injection measurement
+        if (measured_values.has_bus_injection(row)) {
+            flow_sensors[y_bus_structure.bus_entry[row]] = 1;
+        }
+        // upper triangle for branch flow measurement
+        for (Idx ybus_index = y_bus_structure.bus_entry[row] + 1; ybus_index != y_bus_structure.row_indptr[row + 1];
+             ++ybus_index) {
+            for (Idx element_index = y_bus_structure.y_bus_entry_indptr[ybus_index];
+                 element_index != y_bus_structure.y_bus_entry_indptr[ybus_index + 1]; ++element_index) {
+                YBusElement const& element = y_bus_structure.y_bus_element[element_index];
+                // shunt should not be considered
+                // if the branch is fully connected and measured, we consider it as a valid flow sensor
+                // we only need one flow sensor, so the loop will break
+                if (element.element_type != YBusElementType::shunt) {
+                    Idx const branch = element.idx;
+                    if ((measured_values.has_branch_from(branch) || measured_values.has_branch_to(branch)) &&
+                        topo.branch_bus_idx[branch][0] != -1 && topo.branch_bus_idx[branch][1] != -1) {
+                        flow_sensors[ybus_index] = 1;
+                        break;
+                    }
+                }
+            }
+        }
+    }
+    return flow_sensors;
+}
+
+// re-organize the flow sensor for radial grid without phasor measurement
+// this mutate the flow sensor vector to try to assign injection sensor to branch sensor
+// all the branch should be measured if the system is observable
+// this is a sufficient condition check
+// if the grid is not radial, the behavior is undefined.
+inline void assign_injection_sensor_radial(YBusStructure const& y_bus_structure, std::vector<int8_t>& flow_sensors) {
+    Idx const n_bus = std::ssize(y_bus_structure.row_indptr) - 1;
+    // loop the row without the last bus
+    for (Idx row = 0; row != n_bus - 1; ++row) {
+        Idx const current_bus = row;
+        Idx const bus_entry_current = y_bus_structure.bus_entry[current_bus];
+        Idx const branch_entry_upstream = bus_entry_current + 1;
+        // there should be only one upstream branch in the upper diagonal
+        // so the next of branch_entry_upstream is already the next row
+        // because the grid is radial
+        assert(y_bus_structure.row_indptr[current_bus + 1] == branch_entry_upstream + 1);
+        Idx const upstream_bus = y_bus_structure.col_indices[branch_entry_upstream];
+        Idx const bus_entry_upstream = y_bus_structure.bus_entry[upstream_bus];
+        // if the upstream branch is not measured
+        if (flow_sensors[branch_entry_upstream] == 0) {
+            if (flow_sensors[bus_entry_current] == 1) {
+                // try to steal from current bus
+                std::swap(flow_sensors[branch_entry_upstream], flow_sensors[bus_entry_current]);
+            } else if (flow_sensors[bus_entry_upstream] == 1) {
+                // if not possible, steal from upstream bus
+                std::swap(flow_sensors[branch_entry_upstream], flow_sensors[bus_entry_upstream]);
+            }
+        }
+        // remove the current bus injection regardless of the original state
+        flow_sensors[bus_entry_current] = 0;
+    }
+    // set last bus injection to zero
+    flow_sensors[y_bus_structure.bus_entry[n_bus - 1]] = 0;
+}
+
+} // namespace detail
 
-    Idx const n_bus{topo->n_bus()};
-    std::vector<BranchIdx> const& branch_bus_idx{topo->branch_bus_idx};
+template <symmetry_tag sym>
+inline void necessary_observability_check(MeasuredValues<sym> const& measured_values, MathModelTopology const& topo,
+                                          YBusStructure const& y_bus_structure) {
+    Idx const n_bus{topo.n_bus()};
 
     auto const [n_voltage_sensor, n_voltage_phasor_sensor] = detail::count_voltage_sensors(n_bus, measured_values);
     if (n_voltage_sensor < 1) {
-        throw NotObservableError{"no voltage sensor found"};
+        throw NotObservableError{"No voltage sensor found!\n"};
     }
 
-    Idx const n_injection_sensor = detail::count_bus_injection_sensors(n_bus, measured_values);
-    Idx const n_branch_sensor = detail::count_branch_sensors(branch_bus_idx, n_bus, measured_values);
-    Idx const n_power_sensor = n_injection_sensor + n_branch_sensor;
+    std::vector<int8_t> flow_sensors = detail::count_flow_sensors(measured_values, topo, y_bus_structure);
+    // count flow sensors, note we manually specify the intial value type to avoid overflow
+    Idx const n_flow_sensor = std::reduce(flow_sensors.cbegin(), flow_sensors.cend(), Idx{}, std::plus<Idx>{});
 
-    if (n_voltage_phasor_sensor == 0 && n_power_sensor < n_bus - 1) {
+    // check nessessary condition for observability
+    if (n_voltage_phasor_sensor == 0 && n_flow_sensor < n_bus - 1) {
         throw NotObservableError{};
     }
-    if (n_voltage_phasor_sensor > 0 && n_power_sensor + n_voltage_phasor_sensor < n_bus) {
+    if (n_voltage_phasor_sensor > 0 && n_flow_sensor + n_voltage_phasor_sensor < n_bus) {
         throw NotObservableError{};
     }
+
+    // for radial grid without phasor measurement, try to assign injection sensor to branch sensor
+    // we can then check sufficient condition for observability
+    if (topo.is_radial && n_voltage_phasor_sensor == 0) {
+        detail::assign_injection_sensor_radial(y_bus_structure, flow_sensors);
+        // count flow sensors again
+        Idx const n_flow_sensor_new = std::reduce(flow_sensors.cbegin(), flow_sensors.cend(), Idx{}, std::plus<Idx>{});
+        if (n_flow_sensor_new < n_bus - 1) {
+            throw NotObservableError{"The number of power sensors appears sufficient, but they are not independent "
+                                     "enough. The system is still not observable.\n"};
+        }
+    }
 }
 
 } // namespace power_grid_model::math_solver

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/y_bus.hpp

@@ -294,6 +294,7 @@ template <symmetry_tag sym> class YBus {
     }
 
     // getter
+    YBusStructure const& y_bus_structure() const { return *y_bus_struct_; }
     Idx size() const { return static_cast<Idx>(bus_entry().size()); }
     Idx nnz() const { return row_indptr().back(); }
     Idx nnz_lu() const { return row_indptr_lu().back(); }

power_grid_model_c/power_grid_model/include/power_grid_model/topology.hpp

@@ -241,10 +241,12 @@ class Topology {
                 // no cycle, the graph is pure tree structure
                 // just reverse the node
                 std::reverse(dfs_node.begin(), dfs_node.end());
+                math_topo_single.is_radial = true;
             } else {
                 // with cycles, meshed graph
                 // use minimum degree
                 math_topo_single.fill_in = reorder_node(dfs_node, back_edges);
+                math_topo_single.is_radial = false;
             }
             // initialize phase shift
             math_topo_single.phase_shift.resize(dfs_node.size());

tests/cpp_unit_tests/test_observability.cpp

@@ -18,19 +18,21 @@ void check_not_observable(MathModelTopology const& topo, MathModelParam<symmetri
     YBus<symmetric_t> const y_bus{topo_ptr, param_ptr};
     math_solver::MeasuredValues<symmetric_t> const measured_values{y_bus.shared_topology(), se_input};
 
-    CHECK_THROWS_AS(math_solver::necessary_observability_check(measured_values, y_bus.shared_topology()),
-                    NotObservableError);
+    CHECK_THROWS_AS(
+        math_solver::necessary_observability_check(measured_values, y_bus.math_topology(), y_bus.y_bus_structure()),
+        NotObservableError);
 }
 } // namespace
 
 TEST_CASE("Necessary observability check") {
     /*
-                  /-branch_0-\
+                  /-branch_1-\
             bus_2              bus_1 --branch_0-- bus_0 -- source
-                  \-branch_1-/
+                  \-branch_2-/
     */
     MathModelTopology topo;
     topo.slack_bus = 0;
+    topo.is_radial = true;
     topo.phase_shift = {0.0, 0.0, 0.0};
     topo.branch_bus_idx = {{0, 1}, {1, 2}, {1, 2}};
     topo.sources_per_bus = {from_sparse, {0, 1, 1, 1}};
@@ -41,7 +43,7 @@ TEST_CASE("Necessary observability check") {
     topo.power_sensors_per_source = {from_sparse, {0, 0}};
     topo.power_sensors_per_load_gen = {from_sparse, {0, 0, 0}};
     topo.power_sensors_per_shunt = {from_sparse, {0}};
-    topo.power_sensors_per_branch_from = {from_sparse, {0, 0, 1, 1}};
+    topo.power_sensors_per_branch_from = {from_sparse, {0, 1, 1, 1}};
     topo.power_sensors_per_branch_to = {from_sparse, {0, 0, 0, 0}};
     topo.voltage_sensors_per_bus = {from_sparse, {0, 1, 1, 1}};
 
@@ -61,7 +63,8 @@ TEST_CASE("Necessary observability check") {
         auto topo_ptr = std::make_shared<MathModelTopology const>(topo);
         YBus<symmetric_t> const y_bus{topo_ptr, param_ptr};
         math_solver::MeasuredValues<symmetric_t> const measured_values{y_bus.shared_topology(), se_input};
-        CHECK_NOTHROW(math_solver::necessary_observability_check(measured_values, y_bus.shared_topology()));
+        CHECK_NOTHROW(math_solver::necessary_observability_check(measured_values, y_bus.math_topology(),
+                                                                 y_bus.y_bus_structure()));
     }
 
     SUBCASE("No voltage sensor") {
@@ -86,11 +89,21 @@ TEST_CASE("Necessary observability check") {
 
         SUBCASE("Parallel branch get counted as one sensor") {
             // Add sensor on branch 3 to side. Hence 2 parallel sensors
+            topo.power_sensors_per_branch_from = {from_sparse, {0, 0, 1, 1}};
             topo.power_sensors_per_branch_to = {from_sparse, {0, 0, 0, 1}};
             se_input.measured_branch_to_power = {{1.0, 1.0, 1.0}};
             check_not_observable(topo, param, se_input);
         }
     }
+    SUBCASE("Not independent") {
+        // set branch sensor to bus_1 <-branch_1-> bus_2
+        // it is not independent with injection sensor of bus_2
+        topo.power_sensors_per_branch_from = {from_sparse, {0, 0, 1, 1}};
+        // set non phasor measurement
+        se_input.measured_voltage = {{{1.0, nan}, 1.0}};
+        // this will throw NotObservableError
+        check_not_observable(topo, param, se_input);
+    }
 }
 
 } // namespace power_grid_model

tests/data/state_estimation/not-independent-sensor/input.json

@@ -0,0 +1,31 @@
+{
+  "version": "1.0",
+  "type": "input",
+  "is_batch": false,
+  "attributes": {},
+  "data": {
+    "node": [
+      {"id": 0, "u_rated": 10},
+      {"id": 1, "u_rated": 10},
+      {"id": 7, "u_rated": 10}
+    ],
+    "line": [
+      {"id": 2, "from_node": 0, "to_node": 1, "from_status": 1, "to_status": 1, "r1": 1.0, "x1": 0.0, "c1": 0.0, "tan1": 0.0, "i_n": 500},
+      {"id": 8, "from_node": 1, "to_node": 7, "from_status": 1, "to_status": 1, "r1": 1.0, "x1": 0.0, "c1": 0.0, "tan1": 0.0, "i_n": 500}
+    ],
+    "source": [
+      {"id": 3, "node": 0, "status": 1, "u_ref": 1.0}
+    ],
+    "sym_load": [
+      {"id": 4, "node": 1, "status": 1, "type": 0, "p_specified": 9, "q_specified": 0},
+      {"id": 10, "node": 7, "status": 1, "type": 0, "p_specified": 9, "q_specified": 0}
+    ],
+    "sym_power_sensor": [
+      {"id": 5, "measured_object": 3, "measured_terminal_type": 2, "power_sigma": 1, "p_measured": 0, "q_measured": 0},
+      {"id": 9, "measured_object": 2, "measured_terminal_type": 0, "power_sigma": 1, "p_measured": 0, "q_measured": 0}
+    ],
+    "sym_voltage_sensor": [
+      {"id": 6, "measured_object": 0, "u_sigma": 1, "u_measured": 10}
+    ]
+  }
+}

tests/data/state_estimation/not-independent-sensor/input.json.license

@@ -0,0 +1,3 @@
+SPDX-FileCopyrightText: Contributors to the Power Grid Model project <[email protected]>
+
+SPDX-License-Identifier: MPL-2.0

tests/data/state_estimation/not-independent-sensor/params.json

@@ -0,0 +1,12 @@
+{
+  "calculation_method": ["newton_raphson", "iterative_linear"],
+  "rtol": 1e-8,
+  "atol": {
+    "default": 1e-8,
+    ".+_residual": 5e-4
+  },
+  "fail": {
+    "raises": "NotObservableError",
+    "reason": "Power measurements are not independent"
+  }
+}