Merge branch 'feature/current-sensor-measurements-observability' into feature/current-sensor-math-solver

Author: mgovers

power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp

@@ -218,43 +218,75 @@ template <std::ranges::view RandVarsView>
                           UniformComplexRandVar<typename std::ranges::range_value_t<RandVarsView>::sym>> ||
              std::same_as<std::ranges::range_value_t<RandVarsView>,
                           IndependentComplexRandVar<typename std::ranges::range_value_t<RandVarsView>::sym>>
-constexpr auto accumulate(RandVarsView rand_vars) {
+constexpr auto combine(RandVarsView rand_vars) {
     using RandVarType = std::ranges::range_value_t<RandVarsView>;
     using ValueType = decltype(RandVarType::value);
     using VarianceType = decltype(RandVarType::variance);
 
     VarianceType accumulated_inverse_variance{};
-    ValueType accumulated_value{};
+    ValueType weighted_accumulated_value{};
 
-    std::ranges::for_each(rand_vars, [&accumulated_inverse_variance, &accumulated_value](auto const& measurement) {
-        accumulated_inverse_variance += VarianceType{1.0} / measurement.variance;
-        accumulated_value += measurement.value / measurement.variance;
-    });
+    std::ranges::for_each(rand_vars,
+                          [&accumulated_inverse_variance, &weighted_accumulated_value](auto const& measurement) {
+                              accumulated_inverse_variance += VarianceType{1.0} / measurement.variance;
+                              weighted_accumulated_value += measurement.value / measurement.variance;
+                          });
 
-    if (is_normal(accumulated_inverse_variance)) {
-        return RandVarType{.value = accumulated_value / accumulated_inverse_variance,
-                           .variance = VarianceType{1.0} / accumulated_inverse_variance};
+    if (!is_normal(accumulated_inverse_variance)) {
+        return RandVarType{.value = weighted_accumulated_value,
+                           .variance = VarianceType{std::numeric_limits<double>::infinity()}};
     }
-    return RandVarType{.value = accumulated_value, .variance = VarianceType{std::numeric_limits<double>::infinity()}};
+    return RandVarType{.value = weighted_accumulated_value / accumulated_inverse_variance,
+                       .variance = VarianceType{1.0} / accumulated_inverse_variance};
 }
 
 template <std::ranges::view RandVarsView>
     requires std::same_as<std::ranges::range_value_t<RandVarsView>,
                           DecomposedComplexRandVar<typename std::ranges::range_value_t<RandVarsView>::sym>>
-constexpr auto accumulate(RandVarsView rand_vars) {
+constexpr auto combine(RandVarsView rand_vars) {
     using sym = std::ranges::range_value_t<RandVarsView>::sym;
 
     DecomposedComplexRandVar<sym> result{
         .real_component =
-            accumulate(rand_vars | std::views::transform([](auto const& x) -> auto& { return x.real_component; })),
+            combine(rand_vars | std::views::transform([](auto const& x) -> auto& { return x.real_component; })),
         .imag_component =
-            accumulate(rand_vars | std::views::transform([](auto const& x) -> auto& { return x.imag_component; }))};
+            combine(rand_vars | std::views::transform([](auto const& x) -> auto& { return x.imag_component; }))};
 
     if (!(is_normal(result.real_component.variance) && is_normal(result.imag_component.variance))) {
         result.real_component.variance = RealValue<sym>{std::numeric_limits<double>::infinity()};
         result.imag_component.variance = RealValue<sym>{std::numeric_limits<double>::infinity()};
     }
     return result;
 }
+
+namespace detail {
+template <symmetry_tag sym> inline RealValue<sym> cabs_or_real(ComplexValue<sym> const& value) {
+    if (is_nan(imag(value))) {
+        return real(value); // only keep real part
+    }
+    return cabs(value); // get abs of the value
+}
+} // namespace detail
+
+template <std::ranges::view RandVarsView>
+    requires std::same_as<std::ranges::range_value_t<RandVarsView>,
+                          UniformComplexRandVar<typename std::ranges::range_value_t<RandVarsView>::sym>>
+constexpr auto combine_magnitude(RandVarsView rand_vars) {
+    using sym = std::ranges::range_value_t<RandVarsView>::sym;
+
+    auto const weighted_average_magnitude_measurement =
+        statistics::combine(rand_vars | std::views::transform([](auto const& measurement) {
+                                return UniformRealRandVar<sym>{.value = detail::cabs_or_real<sym>(measurement.value),
+                                                               .variance = measurement.variance};
+                            }));
+    return UniformComplexRandVar<sym>{.value =
+                                          [&weighted_average_magnitude_measurement]() {
+                                              ComplexValue<sym> abs_value =
+                                                  piecewise_complex_value<sym>(DoubleComplex{0.0, nan});
+                                              abs_value += weighted_average_magnitude_measurement.value;
+                                              return abs_value;
+                                          }(),
+                                      .variance = weighted_average_magnitude_measurement.variance};
+}
 } // namespace statistics
 } // namespace power_grid_model

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp

@@ -15,16 +15,6 @@ Collect all measured Values
 #include <memory>
 
 namespace power_grid_model::math_solver {
-
-namespace detail {
-template <symmetry_tag sym> inline RealValue<sym> cabs_or_real(ComplexValue<sym> const& value) {
-    if (is_nan(imag(value))) {
-        return real(value); // only keep real part
-    }
-    return cabs(value); // get abs of the value
-}
-} // namespace detail
-
 // processed measurement struct
 // combined all measurement of the same quantity
 // accumulate for bus injection measurement
@@ -71,6 +61,7 @@ template <symmetry_tag sym> class MeasuredValues {
 
     constexpr bool has_angle() const { return n_voltage_angle_measurements_ > 0; }
     constexpr bool has_voltage_measurements() const { return n_voltage_measurements_ > 0; }
+    constexpr bool has_global_angle_current() const { return n_global_angle_current_measurements_ > 0; }
 
     constexpr bool has_voltage(Idx bus) const { return idx_voltage_[bus] >= 0; }
     constexpr bool has_angle_measurement(Idx bus) const { return !is_nan(imag(voltage(bus))); }
@@ -215,6 +206,7 @@ template <symmetry_tag sym> class MeasuredValues {
 
     Idx n_voltage_measurements_{};
     Idx n_voltage_angle_measurements_{};
+    Idx n_global_angle_current_measurements_{};
 
     // average angle shift of voltages with angle measurement
     // default is zero is no voltage has angle measurement
@@ -406,15 +398,15 @@ template <symmetry_tag sym> class MeasuredValues {
         connected at that side. For each branch the checker checks if the from and to side are connected by checking if
         branch_bus_idx = disconnected.
 
-        If the branch_bus_idx = disconnected, idx_branch_(to|from)_(power|current)_ is set to disconnected.
+        If the branch_bus_idx = disconnected, idx_branch_(to/from)_(power/current)_ is set to disconnected.
         If the side is connected, but there are no measurements in this branch side
-        idx_branch_(to|from)_(power|current)_ is set to disconnected.
-        Else, idx_branch_(to|from)_(power|current)_ is set to the index of the aggregated data in
+        idx_branch_(to/from)_(power/current)_ is set to disconnected.
+        Else, idx_branch_(to/from)_(power/current)_ is set to the index of the aggregated data in
         power/current_main_value_.
 
         All measurement values for a single side of a branch are combined in a weighted average, which is appended to
         power/current_main_value_. The values in power/current_main_value_ can be found using
-        idx_branch_(to|from)_(power|current)_.
+        idx_branch_(to/from)_(power/current)_.
         */
         MathModelTopology const& topo = math_topology();
         static constexpr auto branch_from_checker = [](BranchIdx x) { return x[0] != -1; };
@@ -437,12 +429,10 @@ template <symmetry_tag sym> class MeasuredValues {
                 process_one_object(branch, topo.current_sensors_per_branch_to, topo.branch_bus_idx,
                                    input.measured_branch_to_current, current_main_value_, branch_to_checker);
 
-            if ((!has_angle()) && std::ranges::any_of(current_main_value_, [](auto const& measurement) {
+            n_global_angle_current_measurements_ =
+                std::ranges::count_if(current_main_value_, [](auto const& measurement) {
                     return measurement.angle_measurement_type == AngleMeasurementType::global_angle;
-                })) {
-                throw ConflictingAngleMeasurementType{
-                    "Global angle current measurements require a voltage angle measurement in the connected subgrid."};
-            }
+                });
         }
     }
 
@@ -455,28 +445,16 @@ template <symmetry_tag sym> class MeasuredValues {
                                                             IdxRange const& sensors) {
         auto complex_measurements = sensors | std::views::transform([&data](Idx pos) -> auto& { return data[pos]; });
         if constexpr (only_magnitude) {
-            auto const combined_magnitude_measurement = statistics::accumulate(
-                complex_measurements | std::views::transform([](auto const& measurement) {
-                    return UniformRealRandVar<sym>{.value = detail::cabs_or_real<sym>(measurement.value),
-                                                   .variance = measurement.variance};
-                }));
-            return UniformComplexRandVar<sym>{.value =
-                                                  [&combined_magnitude_measurement]() {
-                                                      ComplexValue<sym> abs_value =
-                                                          piecewise_complex_value<sym>(DoubleComplex{0.0, nan});
-                                                      abs_value += combined_magnitude_measurement.value;
-                                                      return abs_value;
-                                                  }(),
-                                              .variance = combined_magnitude_measurement.variance};
+            return statistics::combine_magnitude(complex_measurements);
         } else {
-            return statistics::accumulate(complex_measurements);
+            return statistics::combine(complex_measurements);
         }
     }
     template <bool only_magnitude = false>
         requires(!only_magnitude)
     static PowerSensorCalcParam<sym> combine_measurements(std::vector<PowerSensorCalcParam<sym>> const& data,
                                                           IdxRange const& sensors) {
-        return statistics::accumulate(sensors | std::views::transform([&data](Idx pos) -> auto& { return data[pos]; }));
+        return statistics::combine(sensors | std::views::transform([&data](Idx pos) -> auto& { return data[pos]; }));
     }
     template <bool only_magnitude = false>
         requires(!only_magnitude)
@@ -485,15 +463,15 @@ template <symmetry_tag sym> class MeasuredValues {
         auto const params = sensors | std::views::transform([&data](Idx pos) -> auto& { return data[pos]; });
         auto const angle_measurement_type = sensors.empty() ? AngleMeasurementType::local_angle // fallback
                                                             : params.front().angle_measurement_type;
-        if (std::ranges::any_of(params, [angle_measurement_type](auto const& params) {
-                return params.angle_measurement_type != angle_measurement_type;
+        if (std::ranges::any_of(params, [angle_measurement_type](auto const& param) {
+                return param.angle_measurement_type != angle_measurement_type;
             })) {
             throw ConflictingAngleMeasurementType{
                 "Cannot mix local and global angle current measurements on the same terminal."};
         }
 
         return {.angle_measurement_type = angle_measurement_type,
-                .measurement = statistics::accumulate(
+                .measurement = statistics::combine(
                     params | std::views::transform([](auto const& param) -> auto& { return param.measurement; }))};
     }
 

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

@@ -214,6 +214,8 @@ template <symmetry_tag sym_type> class NewtonRaphsonSESolver {
     typename SparseLUSolver<NRSEGainBlock<sym>, NRSERhs<sym>, NRSEUnknown<sym>>::BlockPermArray perm_;
 
     void initialize_unknown(ComplexValueVector<sym>& initial_u, MeasuredValues<sym> const& measured_values) {
+        using statistics::detail::cabs_or_real;
+
         reset_unknown();
         RealValue<sym> const mean_angle_shift = measured_values.mean_angle_shift();
         for (Idx bus = 0; bus != n_bus_; ++bus) {
@@ -224,7 +226,7 @@ template <symmetry_tag sym_type> class NewtonRaphsonSESolver {
                 if (measured_values.has_angle_measurement(bus)) {
                     estimated_result.theta() = arg(measured_values.voltage(bus));
                 }
-                estimated_result.v() = detail::cabs_or_real<sym>(measured_values.voltage(bus));
+                estimated_result.v() = cabs_or_real<sym>(measured_values.voltage(bus));
             }
             initial_u[bus] = estimated_result.v() * exp(1.0i * estimated_result.theta());
         }
@@ -532,14 +534,16 @@ template <symmetry_tag sym_type> class NewtonRaphsonSESolver {
     /// @param bus bus with voltage measurement
     void process_voltage_measurements(NRSEGainBlock<sym>& block, NRSERhs<sym>& rhs_block,
                                       MeasuredValues<sym> const& measured_values, Idx const& bus) {
+        using statistics::detail::cabs_or_real;
+
         if (!measured_values.has_voltage(bus)) {
             return;
         }
 
         // G += 1.0 / variance
         // for 3x3 tensor, fill diagonal
         auto const w_v = RealTensor<sym>{1.0 / measured_values.voltage_var(bus)};
-        auto const abs_measured_v = detail::cabs_or_real<sym>(measured_values.voltage(bus));
+        auto const abs_measured_v = cabs_or_real<sym>(measured_values.voltage(bus));
         auto const delta_v = abs_measured_v - x_[bus].v();
 
         auto const virtual_angle_measurement_bus = measured_values.has_voltage(math_topo_->slack_bus)

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

@@ -155,6 +155,11 @@ inline ObservabilityResult necessary_observability_check(MeasuredValues<sym> con
         throw NotObservableError{};
     }
 
+    if (measured_values.has_global_angle_current() && n_voltage_phasor_sensor == 0) {
+        throw NotObservableError{
+            "Global angle current sensors require at least one voltage angle measurement as a reference point.\n"};
+    }
+
     // 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) {
@@ -163,8 +168,9 @@ inline ObservabilityResult necessary_observability_check(MeasuredValues<sym> con
         if (Idx const n_flow_sensor_new =
                 std::reduce(flow_sensors.cbegin(), flow_sensors.cend(), Idx{}, std::plus<Idx>{});
             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"};
+            throw NotObservableError{
+                "The number of power/current sensors appears sufficient, but they are not independent "
+                "enough. The system is still not observable.\n"};
         }
         result.is_sufficiently_observable = true;
     }