Merge pull request #896 from PowerGridModel/feature/statistics-module

Measurement statistics module

Author: Jerry-Jinfeng-Guo

power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp

@@ -7,6 +7,7 @@
 #include "common/common.hpp"
 #include "common/enum.hpp"
 #include "common/grouped_index_vector.hpp"
+#include "common/statistics.hpp"
 #include "common/three_phase_tensor.hpp"
 
 namespace power_grid_model {
@@ -78,21 +79,10 @@ template <symmetry_tag sym_type> struct ApplianceShortCircuitSolverOutput {
     ComplexValue<sym> i{};
 };
 
-// Complex measured value of a sensor in p.u. with a uniform variance across all phases and axes of the complex plane
-// (circularly symmetric)
-template <symmetry_tag sym_type> struct UniformComplexRandomVariable {
-    using sym = sym_type;
-
-    static constexpr bool symmetric{is_symmetric_v<sym>};
-
-    ComplexValue<sym> value{};
-    double variance{}; // variance (sigma^2) of the error range, in p.u.
-};
-
 // voltage sensor calculation parameters for state estimation
 // The value is the complex voltage
 // If the imaginary part is NaN, it means the angle calculation is not correct
-template <symmetry_tag sym> using VoltageSensorCalcParam = UniformComplexRandomVariable<sym>;
+template <symmetry_tag sym> using VoltageSensorCalcParam = UniformComplexRandVar<sym>;
 
 // power sensor calculation parameters for state estimation
 // The value is the complex power
@@ -118,8 +108,8 @@ template <symmetry_tag sym_type> struct CurrentSensorCalcParam {
 
     AngleMeasurementType angle_measurement_type{};
     ComplexValue<sym> value{};
-    double i_real_variance{}; // variance (sigma^2) of the error range of real part of the current, in p.u.
-    double i_imag_variance{}; // variance (sigma^2) of the error range of imaginary part of the current, in p.u.
+    RealValue<sym> i_real_variance{}; // variance (sigma^2) of the error range of real part of the current, in p.u.
+    RealValue<sym> i_imag_variance{}; // variance (sigma^2) of the error range of imaginary part of the current, in p.u.
 };
 
 template <typename T>

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

@@ -0,0 +1,207 @@
+// SPDX-FileCopyrightText: Contributors to the Power Grid Model project <powergridmodel@lfenergy.org>
+//
+// SPDX-License-Identifier: MPL-2.0
+
+#pragma once
+
+#include "common.hpp"
+#include "three_phase_tensor.hpp"
+
+/**
+ * @file statistics.hpp
+ * @brief This file contains various structures and functions for handling statistical representations of
+ * random variables (RandVar) used in the Power Grid Model, like in the State Estimation algorithms to
+ * handle measurements.
+ *
+ * The structures provided in this file are used to represent measured values of sensors
+ * with different types of variances. These structures support both symmetric and asymmetric representations and
+ * provide conversion operators to transform between these representations.
+ *
+ * A random variable in PGM can have following characteristics:
+ *  - Uniform: Single total variance for all phases
+ *  - Independent: all phases are independent from each other
+ *  - Scalar: Named as `Real` in this file, a scalar value `RealValue`, eg. real axis: RealValue (* 1), imaginary axis:
+ * RealValue (* i).
+ *  - Complex: a complex value with real and imaginary parts.
+ *
+ * Based on these, we use combined variables in Decomposed/Polar forms:
+ * - Decomposed: treat random variables individually as in cartesian co-ordinates with separated variances for both real
+ * and imaginary part.
+ * - Polar: random variables are in polar co-ordinates, with magnitudes and angles.
+ *
+ */
+
+namespace power_grid_model {
+template <symmetry_tag sym_type> struct UniformRealRandVar {
+    using sym = sym_type;
+
+    RealValue<sym> value{};
+    double variance{}; // variance (sigma^2) of the error range
+
+    explicit operator UniformRealRandVar<asymmetric_t>() const
+        requires(is_symmetric_v<sym>)
+    {
+        return {.value = RealValue<asymmetric_t>{std::piecewise_construct, value}, .variance = variance};
+    }
+    explicit operator UniformRealRandVar<symmetric_t>() const
+        requires(is_asymmetric_v<sym>)
+    {
+        return {.value = mean_val(value), .variance = variance / 3.0};
+    }
+};
+
+template <symmetry_tag sym_type> struct IndependentRealRandVar {
+    using sym = sym_type;
+
+    RealValue<sym> value{};
+    RealValue<sym> variance{}; // variance (sigma^2) of the error range
+
+    explicit operator UniformRealRandVar<symmetric_t>() const {
+        constexpr auto scale = is_asymmetric_v<sym> ? 3.0 : 1.0;
+        return {.value = mean_val(value), .variance = mean_val(variance) / scale};
+    }
+    explicit operator UniformRealRandVar<asymmetric_t>() const {
+        return {.value = value, .variance = mean_val(variance)};
+    }
+    explicit operator IndependentRealRandVar<asymmetric_t>() const
+        requires(is_symmetric_v<sym>)
+    {
+        return {.value = RealValue<asymmetric_t>{std::piecewise_construct, value},
+                .variance = RealValue<asymmetric_t>{std::piecewise_construct, variance}};
+    }
+    explicit operator IndependentRealRandVar<symmetric_t>() const
+        requires(is_asymmetric_v<sym>)
+    {
+        return {.value = mean_val(value), .variance = mean_val(variance) / 3.0};
+    }
+};
+
+// Complex measured value of a sensor with a uniform variance across all phases and axes of the complex plane
+// (rotationally symmetric)
+template <symmetry_tag sym_type> struct UniformComplexRandVar {
+    using sym = sym_type;
+
+    ComplexValue<sym> value{};
+    double variance{}; // variance (sigma^2) of the error range
+};
+
+inline UniformComplexRandVar<symmetric_t> pos_seq(UniformComplexRandVar<asymmetric_t> const& var) {
+    return {.value = pos_seq(var.value), .variance = var.variance / 3.0};
+}
+inline UniformComplexRandVar<asymmetric_t> three_phase(UniformComplexRandVar<symmetric_t> const& var) {
+    return {.value = ComplexValue<asymmetric_t>{var.value}, .variance = var.variance};
+}
+
+// Complex measured value of a sensor with separate variances per phase (but rotationally symmetric in the
+// complex plane)
+template <symmetry_tag sym_type> struct IndependentComplexRandVar {
+    using sym = sym_type;
+
+    ComplexValue<sym> value{};
+    RealValue<sym> variance{}; // variance (sigma^2) of the error range
+
+    explicit operator UniformComplexRandVar<sym>() const {
+        return UniformComplexRandVar<sym>{.value = value, .variance = sum_val(variance)};
+    }
+};
+
+// Complex measured value of a sensor modeled as separate real and imaginary components with independent
+// variances (rotationally symmetric)
+template <symmetry_tag sym_type> struct DecomposedComplexRandVar {
+    using sym = sym_type;
+
+    IndependentRealRandVar<sym> real_component;
+    IndependentRealRandVar<sym> imag_component;
+
+    ComplexValue<sym> value() const { return {real_component.value, imag_component.value}; }
+
+    explicit operator UniformComplexRandVar<sym>() const {
+        return static_cast<UniformComplexRandVar<sym>>(static_cast<IndependentComplexRandVar<sym>>(*this));
+    }
+    explicit operator IndependentComplexRandVar<sym>() const {
+        return IndependentComplexRandVar<sym>{.value = value(),
+                                              .variance = real_component.variance + imag_component.variance};
+    }
+};
+
+// Complex measured value of a sensor in polar coordinates (magnitude and angle)
+// (rotationally symmetric)
+template <symmetry_tag sym_type> struct PolarComplexRandVar {
+    using sym = sym_type;
+
+    UniformRealRandVar<sym> magnitude;
+    UniformRealRandVar<sym> angle;
+
+    ComplexValue<sym> value() const { return magnitude.value * exp(1.0i * angle.value); }
+
+    explicit operator UniformComplexRandVar<sym>() const {
+        return static_cast<UniformComplexRandVar<sym>>(static_cast<IndependentComplexRandVar<sym>>(*this));
+    }
+    explicit operator IndependentComplexRandVar<sym>() const {
+        return IndependentComplexRandVar<sym>{
+            .value = value(), .variance = magnitude.variance + magnitude.value * magnitude.value * angle.variance};
+    }
+
+    // For sym to sym conversion:
+    // Var(I_Re) ≈ Var(I) * cos^2(θ) + Var(θ) * I^2 * sin^2(θ)
+    // Var(I_Im) ≈ Var(I) * sin^2(θ) + Var(θ) * I^2 * cos^2(θ)
+    // For asym to asym conversion:
+    // Var(I_Re,p) ≈ Var(I_p) * cos^2(θ_p) + Var(θ_p) * I_p^2 * sin^2(θ_p)
+    // Var(I_Im,p) ≈ Var(I_p) * sin^2(θ_p) + Var(θ_p) * I_p^2 * cos^2(θ_p)
+    explicit operator DecomposedComplexRandVar<sym>() const {
+        auto const cos_theta = cos(angle.value);
+        auto const sin_theta = sin(angle.value);
+        auto const real_component = magnitude.value * cos_theta;
+        auto const imag_component = magnitude.value * sin_theta;
+        return DecomposedComplexRandVar<sym>{
+            .real_component = {.value = real_component,
+                               .variance = magnitude.variance * cos_theta * cos_theta +
+                                           imag_component * imag_component * angle.variance},
+            .imag_component = {.value = imag_component,
+                               .variance = magnitude.variance * sin_theta * sin_theta +
+                                           real_component * real_component * angle.variance}};
+    }
+
+    // Var(I_Re,p) ≈ Var(I) * cos^2(θ - 2πp/3) + Var(θ) * I^2 * sin^2(θ - 2πp/3)
+    // Var(I_Im,p) ≈ Var(I) * sin^2(θ - 2πp/3) + Var(θ) * I^2 * cos^2(θ - 2πp/3)
+    explicit operator DecomposedComplexRandVar<asymmetric_t>() const
+        requires(is_symmetric_v<sym>)
+    {
+        ComplexValue<asymmetric_t> const unit_complex{exp(1.0i * angle.value)};
+        ComplexValue<asymmetric_t> const complex = magnitude.value * unit_complex;
+        return DecomposedComplexRandVar<asymmetric_t>{
+            .real_component = {.value = real(complex),
+                               .variance = magnitude.variance * real(unit_complex) * real(unit_complex) +
+                                           imag(complex) * imag(complex) * angle.variance},
+            .imag_component = {.value = imag(complex),
+                               .variance = magnitude.variance * imag(unit_complex) * imag(unit_complex) +
+                                           real(complex) * real(complex) * angle.variance}};
+    }
+
+    // Var(I_Re) ≈ (1 / 9) * sum_p(Var(I_p) * cos^2(theta_p + 2 * pi * p / 3) + Var(theta_p) * I_p^2 * sin^2(theta_p + 2
+    // * pi * p / 3))
+    // Var(I_Im) ≈ (1 / 9) * sum_p(Var(I_p) * sin^2(theta_p + 2 * pi * p / 3) + Var(theta_p) * I_p^2 *
+    // cos^2(theta_p + 2 * pi * p / 3))
+    explicit operator DecomposedComplexRandVar<symmetric_t>() const
+        requires(is_asymmetric_v<sym>)
+    {
+        ComplexValue<asymmetric_t> const unit_complex{exp(1.0i * angle.value)};
+        ComplexValue<asymmetric_t> const unit_pos_seq_per_phase{unit_complex(0), a * unit_complex(1),
+                                                                a2 * unit_complex(2)};
+        DoubleComplex const pos_seq_value = pos_seq(magnitude.value * unit_complex);
+        return DecomposedComplexRandVar<symmetric_t>{
+            .real_component = {.value = real(pos_seq_value),
+                               .variance = sum_val(magnitude.variance * real(unit_pos_seq_per_phase) *
+                                                       real(unit_pos_seq_per_phase) +
+                                                   imag(unit_pos_seq_per_phase) * imag(unit_pos_seq_per_phase) *
+                                                       magnitude.value * magnitude.value * angle.variance) /
+                                           9.0},
+            .imag_component = {
+                .value = imag(pos_seq_value),
+                .variance = sum_val(magnitude.variance * imag(unit_pos_seq_per_phase) * imag(unit_pos_seq_per_phase) +
+                                    real(unit_pos_seq_per_phase) * real(unit_pos_seq_per_phase) * magnitude.value *
+                                        magnitude.value * angle.variance) /
+                            9.0}};
+    }
+};
+} // namespace power_grid_model

power_grid_model_c/power_grid_model/include/power_grid_model/component/current_sensor.hpp

@@ -10,6 +10,7 @@
 #include "../common/common.hpp"
 #include "../common/enum.hpp"
 #include "../common/exception.hpp"
+#include "../common/statistics.hpp"
 
 namespace power_grid_model {
 
@@ -48,6 +49,10 @@ class GenericCurrentSensor : public Sensor {
         }
     }
 
+  protected:
+    MeasuredTerminalType terminal_type() const { return terminal_type_; }
+    AngleMeasurementType angle_measurement_type() const { return angle_measurement_type_; }
+
   private:
     MeasuredTerminalType terminal_type_;
     AngleMeasurementType angle_measurement_type_;
@@ -73,12 +78,12 @@ template <symmetry_tag current_sensor_symmetry_> class CurrentSensor : public Ge
 
     explicit CurrentSensor(CurrentSensorInput<current_sensor_symmetry> const& current_sensor_input, double u_rated)
         : GenericCurrentSensor{current_sensor_input},
+          base_current_{base_power_3p * inv_sqrt3 / u_rated},
+          base_current_inv_{1.0 / base_current_},
           i_angle_measured_{current_sensor_input.i_angle_measured},
           i_angle_sigma_{current_sensor_input.i_angle_sigma},
-          base_current_{base_power_3p * inv_sqrt3 / u_rated},
-          base_current_inv_{1.0 / base_current_} {
-        set_current(current_sensor_input);
-
+          i_sigma_{current_sensor_input.i_sigma * base_current_inv_},
+          i_measured_{current_sensor_input.i_measured * base_current_inv_} {
         switch (current_sensor_input.measured_terminal_type) {
             using enum MeasuredTerminalType;
         case branch_from:
@@ -93,14 +98,13 @@ template <symmetry_tag current_sensor_symmetry_> class CurrentSensor : public Ge
     };
 
     UpdateChange update(CurrentSensorUpdate<current_sensor_symmetry> const& update_data) {
-        if (!is_nan(update_data.i_sigma)) {
-            i_sigma_ = update_data.i_sigma * base_current_inv_;
-        }
-        if (!is_nan(update_data.i_angle_sigma)) {
-            i_angle_sigma_ = update_data.i_angle_sigma;
-        }
+        assert(update_data.id == this->id() || is_nan(update_data.id));
+
+        update_real_value<symmetric_t>(update_data.i_sigma, i_sigma_, base_current_inv_);
+        update_real_value<symmetric_t>(update_data.i_angle_sigma, i_angle_sigma_, 1.0);
         update_real_value<current_sensor_symmetry>(update_data.i_measured, i_measured_, base_current_inv_);
         update_real_value<current_sensor_symmetry>(update_data.i_angle_measured, i_angle_measured_, 1.0);
+
         return {false, false};
     }
 
@@ -117,29 +121,24 @@ template <symmetry_tag current_sensor_symmetry_> class CurrentSensor : public Ge
     }
 
   private:
-    RealValue<current_sensor_symmetry> i_measured_{};
-    RealValue<current_sensor_symmetry> i_angle_measured_{};
-    double i_sigma_{};
-    double i_angle_sigma_{};
     double base_current_{};
     double base_current_inv_{};
+    RealValue<current_sensor_symmetry> i_angle_measured_{};
+    double i_angle_sigma_{};
+    double i_sigma_{};
+    RealValue<current_sensor_symmetry> i_measured_{};
 
-    void set_current(CurrentSensorInput<current_sensor_symmetry> const& input) {
-        i_sigma_ = input.i_sigma * base_current_inv_;
-        i_measured_ = input.i_measured * base_current_inv_;
-    }
-
-    // TODO when filling the functions below take in mind that i_angle_sigma is optional
-
-    CurrentSensorCalcParam<symmetric_t> sym_calc_param() const final {
-        CurrentSensorCalcParam<symmetric_t> calc_param{};
-        // TODO
-        return calc_param;
-    }
-    CurrentSensorCalcParam<asymmetric_t> asym_calc_param() const final {
-        CurrentSensorCalcParam<asymmetric_t> calc_param{};
-        // TODO
-        return calc_param;
+    CurrentSensorCalcParam<symmetric_t> sym_calc_param() const final { return calc_decomposed_param<symmetric_t>(); }
+    CurrentSensorCalcParam<asymmetric_t> asym_calc_param() const final { return calc_decomposed_param<asymmetric_t>(); }
+
+    template <symmetry_tag sym_calc> CurrentSensorCalcParam<sym_calc> calc_decomposed_param() const {
+        auto const i_polar = PolarComplexRandVar<current_sensor_symmetry>(
+            {i_measured_, i_sigma_ * i_sigma_}, {i_angle_measured_, i_angle_sigma_ * i_angle_sigma_});
+        auto const i_decomposed = DecomposedComplexRandVar<sym_calc>(i_polar);
+        return CurrentSensorCalcParam<sym_calc>{.angle_measurement_type = angle_measurement_type(),
+                                                .value = i_decomposed.value(),
+                                                .i_real_variance = i_decomposed.real_component.variance,
+                                                .i_imag_variance = i_decomposed.imag_component.variance};
     }
     CurrentSensorOutput<symmetric_t> get_sym_output(ComplexValue<symmetric_t> const& i) const final {
         return get_generic_output<symmetric_t>(i);
@@ -150,8 +149,11 @@ template <symmetry_tag current_sensor_symmetry_> class CurrentSensor : public Ge
     template <symmetry_tag sym_calc>
     CurrentSensorOutput<sym_calc> get_generic_output(ComplexValue<sym_calc> const& i) const {
         CurrentSensorOutput<sym_calc> output{};
-        // TODO
-        (void)i; // Suppress unused variable warning
+        output.id = id();
+        ComplexValue<sym_calc> const i_residual{process_mean_val<sym_calc>(i_measured_ - i) * base_current_};
+        output.energized = 1; // current sensor is always energized
+        output.i_residual = cabs(i_residual);
+        output.i_angle_residual = arg(i_residual);
         return output;
     }
 };

power_grid_model_c/power_grid_model/include/power_grid_model/component/power_sensor.hpp

@@ -87,11 +87,11 @@ template <symmetry_tag power_sensor_symmetry_> class PowerSensor : public Generi
     };
 
     UpdateChange update(PowerSensorUpdate<power_sensor_symmetry> const& update_data) {
+        assert(update_data.id == this->id() || is_nan(update_data.id));
+
         set_power(update_data.p_measured, update_data.q_measured);
 
-        if (!is_nan(update_data.power_sigma)) {
-            apparent_power_sigma_ = update_data.power_sigma * inv_base_power;
-        }
+        update_real_value<symmetric_t>(update_data.power_sigma, apparent_power_sigma_, inv_base_power);
         update_real_value<power_sensor_symmetry>(update_data.p_sigma, p_sigma_, inv_base_power);
         update_real_value<power_sensor_symmetry>(update_data.q_sigma, q_sigma_, inv_base_power);
 

power_grid_model_c/power_grid_model_c/src/handle.hpp

@@ -28,8 +28,8 @@ struct PGM_Handle {
 };
 
 template <class Exception = std::exception, class Functor>
-auto call_with_catch(PGM_Handle* handle, Functor func, PGM_Idx error_code,
-                     std::string_view extra_msg = {}) -> std::invoke_result_t<Functor> {
+inline auto call_with_catch(PGM_Handle* handle, Functor func, PGM_Idx error_code,
+                            std::string_view extra_msg = {}) -> std::invoke_result_t<Functor> {
     if (handle) {
         PGM_clear_error(handle);
     }

tests/cpp_unit_tests/CMakeLists.txt

@@ -12,6 +12,7 @@ set(PROJECT_SOURCES
     "test_component_output.cpp"
     "test_component_update.cpp"
     "test_three_phase_tensor.cpp"
+    "test_statistics.cpp"
     "test_node.cpp"
     "test_line.cpp"
     "test_generic_branch.cpp"