Current sensor component functionality

power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp

@@ -108,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/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,26 @@ template <symmetry_tag current_sensor_symmetry_> class CurrentSensor : public Ge
     }
 
   private:
+    double base_current_{};
+    double base_current_inv_{};
     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_{};
 
-    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(mgovers) when filling the functions below take in mind that i_angle_sigma is optional
 
-    // TODO when filling the functions below take in mind that i_angle_sigma is optional
+    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>(); }
 
-    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;
+    template <symmetry_tag sym_calc> CurrentSensorCalcParam<sym_calc> calc_decomposed_param() const {
+        auto const i_polar = PolarComplexRDV<current_sensor_symmetry>(
+            {i_measured_, i_sigma_ * i_sigma_}, {i_angle_measured_, i_angle_sigma_ * i_angle_sigma_});
+        auto const i_decomposed = DecomposedComplexRDV<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 +151,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/test_current_sensor.cpp

@@ -6,6 +6,8 @@
 
 #include <doctest/doctest.h>
 
+TEST_SUITE_BEGIN("test_current_sensor");
+
 namespace power_grid_model {
 namespace {
 auto const r_nan = RealValue<asymmetric_t>{nan};
@@ -30,7 +32,7 @@ TEST_CASE("Test current sensor") {
         for (auto const terminal_type :
              {MeasuredTerminalType::branch_from, MeasuredTerminalType::branch_to, MeasuredTerminalType::branch3_1,
               MeasuredTerminalType::branch3_2, MeasuredTerminalType::branch3_3}) {
-            CAPTURE(terminal_type);
+            // CAPTURE(terminal_type); // TODO check why this inerferes with normal output
 
             CurrentSensorInput<symmetric_t> sym_current_sensor_input{};
             sym_current_sensor_input.id = 0;
@@ -40,10 +42,14 @@ TEST_CASE("Test current sensor") {
             sym_current_sensor_input.i_sigma = 1.0;
             sym_current_sensor_input.i_measured = 1.0 * 1e3;
             sym_current_sensor_input.i_angle_measured = 0.0;
-            sym_current_sensor_input.i_angle_sigma = nan;
+            sym_current_sensor_input.i_angle_sigma = 0.2;
 
             double const u_rated = 10.0e3;
             double const base_current = base_power_3p / u_rated / sqrt3;
+            double const i_pu = 1.0e3 / base_current;
+            double const i_sigma_pu = 1.0 / base_current;
+            double const i_variance_pu = pow(i_sigma_pu, 2);
+            double const i_angle_variance_pu = pow(0.2, 2);
 
             ComplexValue<symmetric_t> const i_sym = (1.0 * 1e3 + 1i * 0.0) / base_current;
             ComplexValue<asymmetric_t> const i_asym = i_sym * RealValue<asymmetric_t>{1.0};
@@ -60,26 +66,28 @@ TEST_CASE("Test current sensor") {
 
             // Check symmetric sensor output for symmetric parameters
             CHECK(sym_sensor_param.angle_measurement_type == AngleMeasurementType::local);
-            CHECK(sym_sensor_param.i_real_variance == doctest::Approx(0.0));
-            CHECK(sym_sensor_param.i_imag_variance == doctest::Approx(0.0));
-            CHECK(real(sym_sensor_param.value) == doctest::Approx(0.0));
+            CHECK(sym_sensor_param.i_real_variance == doctest::Approx(i_variance_pu));
+            CHECK(sym_sensor_param.i_imag_variance == doctest::Approx(i_angle_variance_pu * i_pu * i_pu));
+            CHECK(real(sym_sensor_param.value) == doctest::Approx(i_pu));
             CHECK(imag(sym_sensor_param.value) == doctest::Approx(0.0));
 
-            CHECK(is_nan(sym_sensor_output.id));
-            CHECK(is_nan(sym_sensor_output.energized));
-            CHECK(is_nan(sym_sensor_output.i_residual));
-            CHECK(is_nan(sym_sensor_output.i_angle_residual));
+            CHECK(sym_sensor_output.id == 0);
+            CHECK(sym_sensor_output.energized == 1);
+            CHECK(sym_sensor_output.i_residual == doctest::Approx(0.0));
+            CHECK(sym_sensor_output.i_angle_residual == doctest::Approx(0.0));
 
             // Check symmetric sensor output for asymmetric parameters
-            CHECK(asym_sensor_param.i_real_variance == doctest::Approx(0.0));
-            CHECK(asym_sensor_param.i_imag_variance == doctest::Approx(0.0));
-            CHECK(real(asym_sensor_param.value[0]) == doctest::Approx(0.0));
-            CHECK(imag(asym_sensor_param.value[1]) == doctest::Approx(0.0));
-
-            CHECK(is_nan(sym_sensor_output_asym_param.id));
-            CHECK(is_nan(sym_sensor_output_asym_param.energized));
-            CHECK(is_nan(sym_sensor_output_asym_param.i_residual[0]));
-            CHECK(is_nan(sym_sensor_output_asym_param.i_angle_residual[1]));
+            CHECK(asym_sensor_param.i_real_variance[0] == doctest::Approx(i_variance_pu));
+            CHECK(asym_sensor_param.i_imag_variance[1] ==
+                  doctest::Approx(i_variance_pu * sin(deg_240) * sin(deg_240) +
+                                  i_angle_variance_pu * i_pu * i_pu * cos(deg_240) * cos(deg_240)));
+            CHECK(real(asym_sensor_param.value[0]) == doctest::Approx(i_pu));
+            CHECK(imag(asym_sensor_param.value[1]) == doctest::Approx(i_pu * sin(deg_240)));
+
+            CHECK(sym_sensor_output_asym_param.id == 0);
+            CHECK(sym_sensor_output_asym_param.energized == 1);
+            CHECK(sym_sensor_output_asym_param.i_residual[0] == doctest::Approx(0.0));
+            CHECK(sym_sensor_output_asym_param.i_angle_residual[1] == doctest::Approx(0.0));
 
             CHECK(sym_current_sensor.get_terminal_type() == terminal_type);
 
@@ -94,6 +102,52 @@ TEST_CASE("Test current sensor") {
                                 InvalidMeasuredTerminalType);
             }
         }
+        SUBCASE("Symmetric calculation parameters") {
+            double const u_rated = 10.0e3;
+            double const base_current = base_power_3p / u_rated / sqrt3;
+
+            CurrentSensor<symmetric_t> sym_current_sensor{
+                {1, 1, MeasuredTerminalType::branch3_1, AngleMeasurementType::local}, u_rated};
+
+            SUBCASE("No phase shift") {
+                sym_current_sensor.update(
+                    {.id = 1, .i_sigma = 1.0, .i_angle_sigma = 0.2, .i_measured = 1.0, .i_angle_measured = 0.0});
+                auto const sym_param = sym_current_sensor.calc_param<symmetric_t>();
+
+                CHECK(sym_param.angle_measurement_type == AngleMeasurementType::local);
+                CHECK(sym_param.i_real_variance == doctest::Approx(pow(1.0 / base_current, 2)));
+                CHECK(sym_param.i_imag_variance == doctest::Approx(pow(0.2 / base_current, 2)));
+                CHECK(real(sym_param.value) == doctest::Approx(1.0 / base_current));
+                CHECK(imag(sym_param.value) == doctest::Approx(0.0 / base_current));
+            }
+
+            SUBCASE("Perpendicular phase shift") {
+                sym_current_sensor.update(
+                    {.id = 1, .i_sigma = 1.0, .i_angle_sigma = 0.2, .i_measured = 1.0, .i_angle_measured = pi / 2});
+                auto const sym_param = sym_current_sensor.calc_param<symmetric_t>();
+
+                CHECK(sym_param.angle_measurement_type == AngleMeasurementType::local);
+                CHECK(sym_param.i_real_variance == doctest::Approx(pow(0.2 / base_current, 2)));
+                CHECK(sym_param.i_imag_variance == doctest::Approx(pow(1.0 / base_current, 2)));
+                CHECK(real(sym_param.value) == doctest::Approx(0.0 / base_current));
+                CHECK(imag(sym_param.value) == doctest::Approx(1.0 / base_current));
+            }
+
+            SUBCASE("45deg phase shift") {
+                using std::numbers::sqrt2;
+                constexpr auto inv_sqrt2 = sqrt2 / 2;
+
+                sym_current_sensor.update(
+                    {.id = 1, .i_sigma = 1.0, .i_angle_sigma = 0.2, .i_measured = 1.0, .i_angle_measured = pi / 4});
+                auto const sym_param = sym_current_sensor.calc_param<symmetric_t>();
+
+                CHECK(sym_param.angle_measurement_type == AngleMeasurementType::local);
+                CHECK(sym_param.i_real_variance == doctest::Approx(1.04 / 2.0 / (base_current * base_current)));
+                CHECK(sym_param.i_imag_variance == doctest::Approx(sym_param.i_real_variance));
+                CHECK(real(sym_param.value) == doctest::Approx(inv_sqrt2 / base_current));
+                CHECK(imag(sym_param.value) == doctest::Approx(real(sym_param.value)));
+            }
+        }
     }
     SUBCASE("Update inverse - sym") {
         constexpr auto i_measured = 1.0;
@@ -234,3 +288,5 @@ TEST_CASE("Test current sensor") {
 }
 
 } // namespace power_grid_model
+
+TEST_SUITE_END();