Merge pull request #947 from PowerGridModel/feature/current-sensor-measurements-observability

Current sensor: statistics + observability

Author: mgovers

power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp

@@ -182,6 +182,10 @@ struct MathModelTopology {
 
     Idx n_bus_power_sensor() const { return power_sensors_per_bus.element_size(); }
 
+    Idx n_branch_from_current_sensor() const { return current_sensors_per_branch_from.element_size(); }
+
+    Idx n_branch_to_current_sensor() const { return current_sensors_per_branch_to.element_size(); }
+
     Idx n_transformer_tap_regulator() const { return tap_regulators_per_branch.element_size(); }
 };
 
@@ -235,6 +239,8 @@ template <symmetry_tag sym_type> struct StateEstimationInput {
     std::vector<PowerSensorCalcParam<sym>> measured_branch_from_power;
     std::vector<PowerSensorCalcParam<sym>> measured_branch_to_power;
     std::vector<PowerSensorCalcParam<sym>> measured_bus_injection;
+    std::vector<CurrentSensorCalcParam<sym>> measured_branch_from_current;
+    std::vector<CurrentSensorCalcParam<sym>> measured_branch_to_current;
 };
 
 struct ShortCircuitInput {

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

@@ -39,13 +39,13 @@ class InvalidArguments : public PowerGridError {
     };
 
     template <std::same_as<TypeValuePair>... Options>
-    InvalidArguments(std::string const& method, std::string const& arguments) {
+    InvalidArguments(std::string_view method, std::string_view arguments) {
         append_msg(std::format("{} is not implemented for {}!\n", method, arguments));
     }
 
     template <class... Options>
         requires(std::same_as<std::remove_cvref_t<Options>, TypeValuePair> && ...)
-    InvalidArguments(std::string const& method, Options const&... options)
+    InvalidArguments(std::string_view method, Options const&... options)
         : InvalidArguments{method, "the following combination of options"} {
         (append_msg(std::format(" {}: {}\n", options.name, options.value)), ...);
     }
@@ -54,7 +54,7 @@ class InvalidArguments : public PowerGridError {
 class MissingCaseForEnumError : public InvalidArguments {
   public:
     template <typename T>
-    MissingCaseForEnumError(std::string const& method, T const& value)
+    MissingCaseForEnumError(std::string_view method, T const& value)
         : InvalidArguments{method,
                            std::format("{} #{}", typeid(T).name(), detail::to_string(static_cast<IntS>(value)))} {}
 };
@@ -97,7 +97,7 @@ class InvalidTransformerClock : public PowerGridError {
 
 class SparseMatrixError : public PowerGridError {
   public:
-    SparseMatrixError(Idx err, std::string const& msg = "") {
+    SparseMatrixError(Idx err, std::string_view msg = "") {
         append_msg(
             std::format("Sparse matrix error with error code #{} (possibly singular)\n", detail::to_string(err)));
         if (!msg.empty()) {
@@ -117,7 +117,7 @@ class SparseMatrixError : public PowerGridError {
 
 class NotObservableError : public PowerGridError {
   public:
-    NotObservableError(std::string const& msg = "") {
+    NotObservableError(std::string_view msg = "") {
         append_msg("Not enough measurements available for state estimation.\n");
         if (!msg.empty()) {
             append_msg(std::format("{}\n", msg));
@@ -137,7 +137,7 @@ class IterationDiverge : public PowerGridError {
 
 class MaxIterationReached : public IterationDiverge {
   public:
-    MaxIterationReached(std::string const& msg = "") {
+    MaxIterationReached(std::string_view msg = "") {
         append_msg(std::format("Maximum number of iterations reached{}\n", msg));
     }
 };
@@ -161,25 +161,25 @@ class Idx2DNotFound : public PowerGridError {
 
 class InvalidMeasuredObject : public PowerGridError {
   public:
-    InvalidMeasuredObject(std::string const& object, std::string const& sensor) {
+    InvalidMeasuredObject(std::string_view object, std::string_view sensor) {
         append_msg(std::format("{} measurement is not supported for object of type {}", sensor, object));
     }
 };
 
 class InvalidMeasuredTerminalType : public PowerGridError {
   public:
-    InvalidMeasuredTerminalType(MeasuredTerminalType const terminal_type, std::string const& sensor) {
+    InvalidMeasuredTerminalType(MeasuredTerminalType const terminal_type, std::string_view sensor) {
         append_msg(std::format("{} measurement is not supported for object of type {}", sensor,
                                detail::to_string(static_cast<IntS>(terminal_type))));
     }
 };
 
 class InvalidRegulatedObject : public PowerGridError {
   public:
-    InvalidRegulatedObject(std::string const& object, std::string const& regulator) {
+    InvalidRegulatedObject(std::string_view object, std::string_view regulator) {
         append_msg(std::format("{} regulator is not supported for object of type {}", regulator, object));
     }
-    InvalidRegulatedObject(ID id, std::string const& regulator) {
+    InvalidRegulatedObject(ID id, std::string_view regulator) {
         append_msg(
             std::format("{} regulator is not supported for object with ID {}", regulator, detail::to_string(id)));
     }
@@ -203,7 +203,7 @@ class AutomaticTapCalculationError : public PowerGridError {
 
 class AutomaticTapInputError : public PowerGridError {
   public:
-    AutomaticTapInputError(std::string const& msg) {
+    AutomaticTapInputError(std::string_view msg) {
         append_msg(std::format("Automatic tap changer has invalid configuration. {}", msg));
     }
 };
@@ -215,14 +215,21 @@ class IDWrongType : public PowerGridError {
     }
 };
 
+class ConflictingAngleMeasurementType : public PowerGridError {
+  public:
+    ConflictingAngleMeasurementType(std::string_view msg) {
+        append_msg(std::format("Conflicting angle measurement type. {}", msg));
+    }
+};
+
 class CalculationError : public PowerGridError {
   public:
-    explicit CalculationError(std::string const& msg) { append_msg(msg); }
+    explicit CalculationError(std::string_view msg) { append_msg(msg); }
 };
 
 class BatchCalculationError : public CalculationError {
   public:
-    BatchCalculationError(std::string const& msg, IdxVector failed_scenarios, std::vector<std::string> err_msgs)
+    BatchCalculationError(std::string_view msg, IdxVector failed_scenarios, std::vector<std::string> err_msgs)
         : CalculationError(msg), failed_scenarios_{std::move(failed_scenarios)}, err_msgs_(std::move(err_msgs)) {}
 
     IdxVector const& failed_scenarios() const { return failed_scenarios_; }
@@ -270,12 +277,12 @@ class InvalidShortCircuitPhaseOrType : public PowerGridError {
 
 class SerializationError : public PowerGridError {
   public:
-    explicit SerializationError(std::string const& msg) { append_msg(msg); }
+    explicit SerializationError(std::string_view msg) { append_msg(msg); }
 };
 
 class DatasetError : public PowerGridError {
   public:
-    explicit DatasetError(std::string const& msg) { append_msg(std::format("Dataset error: {}", msg)); }
+    explicit DatasetError(std::string_view msg) { append_msg(std::format("Dataset error: {}", msg)); }
 };
 
 class ExperimentalFeature : public InvalidArguments {
@@ -289,7 +296,7 @@ class NotImplementedError : public PowerGridError {
 
 class UnreachableHit : public PowerGridError {
   public:
-    UnreachableHit(std::string const& method, std::string const& reason_for_assumption) {
+    UnreachableHit(std::string_view method, std::string_view reason_for_assumption) {
         append_msg(std::format("Unreachable code hit when executing {}.\n The following assumption for unreachability "
                                "was not met: {}.\n This may be a bug in the library\n",
                                method, reason_for_assumption));
@@ -299,7 +306,7 @@ class UnreachableHit : public PowerGridError {
 class TapSearchStrategyIncompatibleError : public InvalidArguments {
   public:
     template <typename T1, typename T2>
-    TapSearchStrategyIncompatibleError(std::string const& method, T1 const& value1, T2 const& value2)
+    TapSearchStrategyIncompatibleError(std::string_view method, T1 const& value1, T2 const& value2)
         : InvalidArguments{method, std::format("{} #{} and {} #{}", typeid(T1).name(),
                                                detail::to_string(static_cast<IntS>(value1)), typeid(T2).name(),
                                                detail::to_string(static_cast<IntS>(value2)))} {}

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

@@ -7,6 +7,8 @@
 #include "common.hpp"
 #include "three_phase_tensor.hpp"
 
+#include <ranges>
+
 /**
  * @file statistics.hpp
  * @brief This file contains various structures and functions for handling statistical representations of
@@ -204,4 +206,87 @@ template <symmetry_tag sym_type> struct PolarComplexRandVar {
                             9.0}};
     }
 };
+
+namespace statistics {
+// combine multiple random variables of one quantity using Kalman filter
+template <std::ranges::view RandVarsView>
+    requires std::same_as<std::ranges::range_value_t<RandVarsView>,
+                          UniformRealRandVar<typename std::ranges::range_value_t<RandVarsView>::sym>> ||
+             std::same_as<std::ranges::range_value_t<RandVarsView>,
+                          IndependentRealRandVar<typename std::ranges::range_value_t<RandVarsView>::sym>> ||
+             std::same_as<std::ranges::range_value_t<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 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 weighted_accumulated_value{};
+
+    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 = weighted_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 combine(RandVarsView rand_vars) {
+    using sym = std::ranges::range_value_t<RandVarsView>::sym;
+
+    DecomposedComplexRandVar<sym> result{
+        .real_component =
+            combine(rand_vars | std::views::transform([](auto const& x) -> auto& { return x.real_component; })),
+        .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/main_model_impl.hpp

@@ -1171,6 +1171,8 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
             se_input[i].measured_branch_from_power.resize(state.math_topology[i]->n_branch_from_power_sensor());
             se_input[i].measured_branch_to_power.resize(state.math_topology[i]->n_branch_to_power_sensor());
             se_input[i].measured_bus_injection.resize(state.math_topology[i]->n_bus_power_sensor());
+            se_input[i].measured_branch_from_current.resize(state.math_topology[i]->n_branch_from_current_sensor());
+            se_input[i].measured_branch_to_current.resize(state.math_topology[i]->n_branch_to_current_sensor());
         }
 
         prepare_input_status<sym, &StateEstimationInput<sym>::shunt_status, Shunt>(state, state.topo_comp_coup->shunt,
@@ -1217,6 +1219,22 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
             state, state.topo_comp_coup->power_sensor, se_input,
             [&state](Idx i) { return state.comp_topo->power_sensor_terminal_type[i] == MeasuredTerminalType::node; });
 
+        prepare_input<StateEstimationInput<sym>, CurrentSensorCalcParam<sym>,
+                      &StateEstimationInput<sym>::measured_branch_from_current, GenericCurrentSensor>(
+            state, state.topo_comp_coup->current_sensor, se_input, [&state](Idx i) {
+                using enum MeasuredTerminalType;
+                return state.comp_topo->current_sensor_terminal_type[i] == branch_from ||
+                       // all branch3 sensors are at from side in the mathematical model
+                       state.comp_topo->current_sensor_terminal_type[i] == branch3_1 ||
+                       state.comp_topo->current_sensor_terminal_type[i] == branch3_2 ||
+                       state.comp_topo->current_sensor_terminal_type[i] == branch3_3;
+            });
+        prepare_input<StateEstimationInput<sym>, CurrentSensorCalcParam<sym>,
+                      &StateEstimationInput<sym>::measured_branch_to_current, GenericCurrentSensor>(
+            state, state.topo_comp_coup->current_sensor, se_input, [&state](Idx i) {
+                return state.comp_topo->current_sensor_terminal_type[i] == MeasuredTerminalType::branch_to;
+            });
+
         return se_input;
     }
 

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

@@ -142,7 +142,8 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
 
   private:
     // array selection function pointer
-    static constexpr std::array has_branch_{&MeasuredValues<sym>::has_branch_from, &MeasuredValues<sym>::has_branch_to};
+    static constexpr std::array has_branch_power_{&MeasuredValues<sym>::has_branch_from_power,
+                                                  &MeasuredValues<sym>::has_branch_to_power};
     static constexpr std::array branch_power_{&MeasuredValues<sym>::branch_from_power,
                                               &MeasuredValues<sym>::branch_to_power};
 
@@ -210,7 +211,7 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
                         // measured at from-side: 0, to-side: 1
                         for (IntS const measured_side : std::array<IntS, 2>{0, 1}) {
                             // has measurement
-                            if (std::invoke(has_branch_[measured_side], measured_value, obj)) {
+                            if (std::invoke(has_branch_power_[measured_side], measured_value, obj)) {
                                 // G += Y{side, b0}^H * (variance^-1) * Y{side, b1}
                                 auto const& power = std::invoke(branch_power_[measured_side], measured_value, obj);
                                 block.g() +=
@@ -301,7 +302,7 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
                     // measured at from-side: 0, to-side: 1
                     for (IntS const measured_side : std::array<IntS, 2>{0, 1}) {
                         // has measurement
-                        if (std::invoke(has_branch_[measured_side], measured_value, obj)) {
+                        if (std::invoke(has_branch_power_[measured_side], measured_value, obj)) {
                             PowerSensorCalcParam<sym> const& m =
                                 std::invoke(branch_power_[measured_side], measured_value, obj);
                             // the current needs to be calculated with the voltage of the measured bus side