fix residuals for global angle

Signed-off-by: Santiago Figueroa Manrique <figueroa1395@gmail.com>

Author: figueroa1395

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

@@ -147,10 +147,15 @@ template <symmetry_tag current_sensor_symmetry_> class CurrentSensor : public Ge
     CurrentSensorOutput<sym_calc> get_generic_output(ComplexValue<sym_calc> const& i) const {
         CurrentSensorOutput<sym_calc> output{};
         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);
+        auto const i_calc_param = calc_param<sym_calc>();
+        auto const i_measured_complex = i_calc_param.measurement.value();
+        if (i_calc_param.angle_measurement_type == AngleMeasurementType::global_angle) {
+            output.i_residual = (cabs(i_measured_complex) - cabs(i)) * base_current_;
+            output.i_angle_residual = arg(i_measured_complex) - arg(i);
+        } else {
+            throw NotImplementedError();
+        }
         return output;
     }
 };

tests/cpp_unit_tests/test_current_sensor.cpp

@@ -33,26 +33,36 @@ TEST_CASE("Test current sensor") {
              {MeasuredTerminalType::branch_from, MeasuredTerminalType::branch_to, MeasuredTerminalType::branch3_1,
               MeasuredTerminalType::branch3_2, MeasuredTerminalType::branch3_3}) {
             CAPTURE(terminal_type);
+            // for (auto const angle_measurement_type : {AngleMeasurementType::global_angle,
+            //                                           AngleMeasurementType::local_angle}) {
+            // }
+            auto const angle_measurement_type = AngleMeasurementType::global_angle;
+            CAPTURE(angle_measurement_type);
 
             CurrentSensorInput<symmetric_t> sym_current_sensor_input{};
             sym_current_sensor_input.id = 0;
             sym_current_sensor_input.measured_object = 1;
             sym_current_sensor_input.measured_terminal_type = terminal_type;
-            sym_current_sensor_input.angle_measurement_type = AngleMeasurementType::local_angle;
+            sym_current_sensor_input.angle_measurement_type = angle_measurement_type;
             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_measured = pi / 4.;
             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);
+            double const i_variance_pu = i_sigma_pu * i_sigma_pu;
+            double const i_angle = pi / 4.;
+            double const i_angle_sigma_pi = 0.2;
+            double const i_angle_variance_pu = i_angle_sigma_pi * i_angle_sigma_pi;
 
-            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};
+            auto const i_sym = ComplexValue<symmetric_t>{(1e3 * cos(i_angle) + 1i * 1e3 * sin(i_angle)) / base_current};
+            auto const i_asym = ComplexValue<asymmetric_t>{
+                (1e3 * cos(i_angle) + 1i * 1e3 * sin(i_angle)) / base_current,
+                (1e3 * cos(i_angle + deg_240) + 1i * 1e3 * sin(i_angle + deg_240)) / base_current,
+                (1e3 * cos(i_angle + deg_120) + 1i * 1e3 * sin(i_angle + deg_120)) / base_current};
 
             CurrentSensor<symmetric_t> const sym_current_sensor{sym_current_sensor_input, u_rated};
 
@@ -65,60 +75,73 @@ TEST_CASE("Test current sensor") {
                 sym_current_sensor.get_output<asymmetric_t>(i_asym);
 
             // Check symmetric sensor output for symmetric parameters
-            CHECK(sym_sensor_param.angle_measurement_type == AngleMeasurementType::local_angle);
-            CHECK(sym_sensor_param.measurement.real_component.variance == doctest::Approx(i_variance_pu));
+            if constexpr (angle_measurement_type == AngleMeasurementType::global_angle) {
+                CHECK(sym_sensor_param.angle_measurement_type == AngleMeasurementType::global_angle);
+            } else {
+                CHECK(sym_sensor_param.angle_measurement_type == AngleMeasurementType::local_angle);
+            }
+            // Var(I_Re) ≈ Var(I) * cos^2(pi/4) + Var(θ) * I^2 * sin^2(pi/4)
+            CHECK(sym_sensor_param.measurement.real_component.variance ==
+                  doctest::Approx(0.5 * (i_variance_pu + i_angle_variance_pu * i_pu * i_pu)));
+            // Var(I_Im) ≈ Var(I) * sin^2(pi/4) + Var(θ) * I^2 * cos^2(pi/4)
             CHECK(sym_sensor_param.measurement.imag_component.variance ==
-                  doctest::Approx(i_angle_variance_pu * i_pu * i_pu));
-            CHECK(real(sym_sensor_param.measurement.value()) == doctest::Approx(i_pu));
-            CHECK(imag(sym_sensor_param.measurement.value()) == doctest::Approx(0.0));
+                  doctest::Approx(0.5 * (i_variance_pu + i_angle_variance_pu * i_pu * i_pu)));
+            CHECK(real(sym_sensor_param.measurement.value()) == doctest::Approx(i_pu * cos(i_angle)));
+            CHECK(imag(sym_sensor_param.measurement.value()) == doctest::Approx(i_pu * sin(i_angle)));
 
             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.measurement.real_component.variance[0] == doctest::Approx(i_variance_pu));
+            CHECK(asym_sensor_param.measurement.real_component.variance[0] ==
+                  doctest::Approx(0.5 * (i_variance_pu + i_angle_variance_pu * i_pu * i_pu)));
+            auto const shifted_i_angle = i_angle + deg_240;
             CHECK(asym_sensor_param.measurement.imag_component.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.measurement.value()[0]) == doctest::Approx(i_pu));
-            CHECK(imag(asym_sensor_param.measurement.value()[1]) == doctest::Approx(i_pu * sin(deg_240)));
+                  doctest::Approx(i_variance_pu * sin(shifted_i_angle) * sin(shifted_i_angle) +
+                                  i_angle_variance_pu * i_pu * i_pu * cos(shifted_i_angle) * cos(shifted_i_angle)));
+            CHECK(real(asym_sensor_param.measurement.value()[0]) == doctest::Approx(i_pu * cos(i_angle)));
+            CHECK(imag(asym_sensor_param.measurement.value()[1]) == doctest::Approx(i_pu * sin(shifted_i_angle)));
 
             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));
+            for (auto i = 0; i < 3; ++i) {
+                CHECK(sym_sensor_output_asym_param.i_residual[i] == doctest::Approx(0.0));
+                CHECK(sym_sensor_output_asym_param.i_angle_residual[i] == doctest::Approx(0.0));
+            }
 
             CHECK(sym_current_sensor.get_terminal_type() == terminal_type);
 
-            CHECK(sym_current_sensor.get_angle_measurement_type() == AngleMeasurementType::local_angle);
+            CHECK(sym_current_sensor.get_angle_measurement_type() == AngleMeasurementType::global_angle);
         }
         SUBCASE("Wrong measured terminal type") {
             for (auto const terminal_type :
                  {MeasuredTerminalType::source, MeasuredTerminalType::shunt, MeasuredTerminalType::load,
                   MeasuredTerminalType::generator, MeasuredTerminalType::node}) {
-                CHECK_THROWS_AS((CurrentSensor<symmetric_t>{
-                                    {1, 1, terminal_type, AngleMeasurementType::local_angle, 1.0, 1.0, 1.0, 1.0}, 1.0}),
-                                InvalidMeasuredTerminalType);
+                CHECK_THROWS_AS(
+                    (CurrentSensor<symmetric_t>{
+                        {1, 1, terminal_type, AngleMeasurementType::global_angle, 1.0, 1.0, 1.0, 1.0}, 1.0}),
+                    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{{.id = 1,
-                                                           .measured_object = 1,
-                                                           .measured_terminal_type = MeasuredTerminalType::branch3_1,
-                                                           .angle_measurement_type = AngleMeasurementType::local_angle},
-                                                          u_rated};
+            CurrentSensor<symmetric_t> sym_current_sensor{
+                {.id = 1,
+                 .measured_object = 1,
+                 .measured_terminal_type = MeasuredTerminalType::branch3_1,
+                 .angle_measurement_type = AngleMeasurementType::global_angle},
+                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_angle);
+                CHECK(sym_param.angle_measurement_type == AngleMeasurementType::global_angle);
                 CHECK(sym_param.measurement.real_component.variance == doctest::Approx(pow(1.0 / base_current, 2)));
                 CHECK(sym_param.measurement.imag_component.variance == doctest::Approx(pow(0.2 / base_current, 2)));
                 CHECK(real(sym_param.measurement.value()) == doctest::Approx(1.0 / base_current));
@@ -130,7 +153,7 @@ TEST_CASE("Test current sensor") {
                     {.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_angle);
+                CHECK(sym_param.angle_measurement_type == AngleMeasurementType::global_angle);
                 CHECK(sym_param.measurement.real_component.variance == doctest::Approx(pow(0.2 / base_current, 2)));
                 CHECK(sym_param.measurement.imag_component.variance == doctest::Approx(pow(1.0 / base_current, 2)));
                 CHECK(real(sym_param.measurement.value()) == doctest::Approx(0.0 / base_current));
@@ -145,7 +168,7 @@ TEST_CASE("Test current sensor") {
                     {.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_angle);
+                CHECK(sym_param.angle_measurement_type == AngleMeasurementType::global_angle);
                 CHECK(sym_param.measurement.real_component.variance ==
                       doctest::Approx(1.04 / 2.0 / (base_current * base_current)));
                 CHECK(sym_param.measurement.imag_component.variance ==
@@ -164,7 +187,7 @@ TEST_CASE("Test current sensor") {
         CurrentSensor<symmetric_t> const current_sensor{{.id = 1,
                                                          .measured_object = 1,
                                                          .measured_terminal_type = MeasuredTerminalType::branch3_1,
-                                                         .angle_measurement_type = AngleMeasurementType::local_angle,
+                                                         .angle_measurement_type = AngleMeasurementType::global_angle,
                                                          .i_sigma = i_sigma,
                                                          .i_angle_sigma = i_angle_sigma,
                                                          .i_measured = i_measured,
@@ -289,7 +312,7 @@ TEST_CASE("Test current sensor") {
         CurrentSensor<asymmetric_t> const current_sensor{{.id = 1,
                                                           .measured_object = 1,
                                                           .measured_terminal_type = measured_terminal_type,
-                                                          .angle_measurement_type = AngleMeasurementType::local_angle,
+                                                          .angle_measurement_type = AngleMeasurementType::global_angle,
                                                           .i_sigma = i_sigma,
                                                           .i_angle_sigma = i_angle_sigma,
                                                           .i_measured = i_measured,