From 6ae6f6d4abd402f6a5f6f00e1e2cdf8a0f681d29 Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Tue, 29 Apr 2025 14:21:15 +0200 Subject: [PATCH 1/2] voltage sensor angle residual mod 2pi Signed-off-by: Martijn Govers --- docs/user_manual/components.md | 4 +- .../component/voltage_sensor.hpp | 4 +- tests/cpp_unit_tests/test_voltage_sensor.cpp | 41 +++++++++++++++++++ 3 files changed, 46 insertions(+), 3 deletions(-) diff --git a/docs/user_manual/components.md b/docs/user_manual/components.md index 10418ce73..b23829afa 100644 --- a/docs/user_manual/components.md +++ b/docs/user_manual/components.md @@ -627,10 +627,12 @@ A sensor only has output for state estimation. For other calculation types, sens $$ \begin{eqnarray} & u_{\text{residual}} = u_{\text{measured}} - u_{\text{state}} \\ - & \theta_{\text{residual}} = \theta_{\text{measured}} - \theta_{\text{state}} + & \theta_{\text{residual}} = \theta_{\text{measured}} - \theta_{\text{state}} \pmod{2 \pi} \end{eqnarray} $$ +The $\pmod{2\pi}$ is handled such that $-\pi \lt \theta_{\text{angle},\text{residual}} \leq \pi$. + ### Generic Power Sensor * type name: `generic_power_sensor` diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/component/voltage_sensor.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/component/voltage_sensor.hpp index 276bec23e..8f21f2a88 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/component/voltage_sensor.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/component/voltage_sensor.hpp @@ -151,7 +151,7 @@ template class VoltageSensor : public GenericVoltageSensor { } else { value.u_residual = (real(u1_measured) - cabs(u)) * u_rated_; } - value.u_angle_residual = arg(u1_measured) - arg(u); + value.u_angle_residual = arg(ComplexValue{exp(1.0i * (arg(u1_measured) - arg(u)))}); return value; } @@ -160,7 +160,7 @@ template class VoltageSensor : public GenericVoltageSensor { value.id = id(); value.energized = 1; value.u_residual = (u_measured_ - cabs(u)) * u_rated_ / sqrt3; - value.u_angle_residual = u_angle_measured_ - arg(u); + value.u_angle_residual = arg(ComplexValue{exp(1.0i * (u_angle_measured_ - arg(u)))}); return value; } }; diff --git a/tests/cpp_unit_tests/test_voltage_sensor.cpp b/tests/cpp_unit_tests/test_voltage_sensor.cpp index 5c926a8f1..8bfee4807 100644 --- a/tests/cpp_unit_tests/test_voltage_sensor.cpp +++ b/tests/cpp_unit_tests/test_voltage_sensor.cpp @@ -363,6 +363,47 @@ TEST_CASE("Test voltage sensor") { CHECK(sym_voltage_sensor_asym_output.u_angle_residual[2] == doctest::Approx(-0.2)); } + SUBCASE("Angle = ± pi ∓ 0.1") { + RealValue const u_measured{10.1e3}; + RealValue const u_angle_measured{pi - 0.1}; + double const u_sigma = 1.0; + double const u_rated = 10.0e3; + + VoltageSensorInput voltage_sensor_input{}; + voltage_sensor_input.id = 0; + voltage_sensor_input.measured_object = 1; + voltage_sensor_input.u_sigma = u_sigma; + voltage_sensor_input.u_measured = u_measured; + voltage_sensor_input.u_angle_measured = u_angle_measured; + + VoltageSensor const voltage_sensor{voltage_sensor_input, u_rated}; + + ComplexValue const u_calc_sym{1.02 * exp(1i * (-pi + 0.1))}; + VoltageSensorOutput sym_voltage_sensor_sym_output = + voltage_sensor.get_output(u_calc_sym); + + ComplexValue const u_calc_asym{1.02 * exp(1i * (-pi + 0.1)), 1.03 * exp(1i * (-pi + 0.2)), + 1.04 * exp(1i * (-pi + 0.3))}; + VoltageSensorOutput sym_voltage_sensor_asym_output = + voltage_sensor.get_output(u_calc_asym); + + // Check sym output + CHECK(sym_voltage_sensor_sym_output.id == 0); + CHECK(sym_voltage_sensor_sym_output.energized == 1); + CHECK(sym_voltage_sensor_sym_output.u_residual == doctest::Approx(-100.0)); + CHECK(sym_voltage_sensor_sym_output.u_angle_residual == doctest::Approx(-0.2).epsilon(1e-12)); + + // Check asym output + CHECK(sym_voltage_sensor_asym_output.id == 0); + CHECK(sym_voltage_sensor_asym_output.energized == 1); + CHECK(sym_voltage_sensor_asym_output.u_residual[0] == doctest::Approx(-100.0 / sqrt3)); + CHECK(sym_voltage_sensor_asym_output.u_residual[1] == doctest::Approx(-200.0 / sqrt3)); + CHECK(sym_voltage_sensor_asym_output.u_residual[2] == doctest::Approx(-300.0 / sqrt3)); + CHECK(sym_voltage_sensor_asym_output.u_angle_residual[0] == doctest::Approx(-0.2).epsilon(1e-12)); + CHECK(sym_voltage_sensor_asym_output.u_angle_residual[1] == doctest::Approx(-0.3)); + CHECK(sym_voltage_sensor_asym_output.u_angle_residual[2] == doctest::Approx(-0.4)); + } + SUBCASE("Angle = nan") { RealValue const u_measured{10.1e3}; RealValue const u_angle_measured{nan}; From ab9cff1a3e452e81cdb19a0f704a9ed8ab5900d4 Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Tue, 29 Apr 2025 16:28:35 +0200 Subject: [PATCH 2/2] move logic to common phase_mod_2pi functionality Signed-off-by: Martijn Govers --- .../common/three_phase_tensor.hpp | 8 +++++ .../component/current_sensor.hpp | 4 +-- .../component/voltage_sensor.hpp | 4 +-- .../test_three_phase_tensor.cpp | 32 +++++++++++++++++++ 4 files changed, 43 insertions(+), 5 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/common/three_phase_tensor.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/common/three_phase_tensor.hpp index 0d6ff9a79..252e214f2 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/common/three_phase_tensor.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/common/three_phase_tensor.hpp @@ -176,6 +176,14 @@ inline ComplexValue phase_shift(ComplexValue const& return {phase_shift(m(0)), phase_shift(m(1)), phase_shift(m(2))}; } +// arg(e^(i * phase)) = phase (mod 2pi). By convention restrict to [-pi, pi]. +inline auto phase_mod_2pi(double phase) { + return RealValue{arg(ComplexValue{exp(1.0i * phase)})}; +} +inline auto phase_mod_2pi(RealValue const& phase) { + return RealValue{arg(ComplexValue{exp(1.0i * phase)})}; +} + // calculate kron product of two vector inline double vector_outer_product(double x, double y) { return x * y; } inline DoubleComplex vector_outer_product(DoubleComplex x, DoubleComplex y) { return x * y; } diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/component/current_sensor.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/component/current_sensor.hpp index 5074492ef..0cdf37bab 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/component/current_sensor.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/component/current_sensor.hpp @@ -174,9 +174,7 @@ template class CurrentSensor : public Ge } }(); output.i_residual = (cabs(i_measured_complex) - cabs(i_output)) * base_current_; - // arg(e^i * u_angle) = u_angle in (-pi, pi] - auto const unbounded_i_angle_residual = arg(i_measured_complex) - arg(i_output); - output.i_angle_residual = arg(ComplexValue{exp(1.0i * unbounded_i_angle_residual)}); + output.i_angle_residual = phase_mod_2pi(arg(i_measured_complex) - arg(i_output)); return output; } }; diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/component/voltage_sensor.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/component/voltage_sensor.hpp index 8f21f2a88..6eb264ace 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/component/voltage_sensor.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/component/voltage_sensor.hpp @@ -151,7 +151,7 @@ template class VoltageSensor : public GenericVoltageSensor { } else { value.u_residual = (real(u1_measured) - cabs(u)) * u_rated_; } - value.u_angle_residual = arg(ComplexValue{exp(1.0i * (arg(u1_measured) - arg(u)))}); + value.u_angle_residual = phase_mod_2pi(arg(u1_measured) - arg(u)); return value; } @@ -160,7 +160,7 @@ template class VoltageSensor : public GenericVoltageSensor { value.id = id(); value.energized = 1; value.u_residual = (u_measured_ - cabs(u)) * u_rated_ / sqrt3; - value.u_angle_residual = arg(ComplexValue{exp(1.0i * (u_angle_measured_ - arg(u)))}); + value.u_angle_residual = phase_mod_2pi(u_angle_measured_ - arg(u)); return value; } }; diff --git a/tests/cpp_unit_tests/test_three_phase_tensor.cpp b/tests/cpp_unit_tests/test_three_phase_tensor.cpp index 42d805020..48167f4ca 100644 --- a/tests/cpp_unit_tests/test_three_phase_tensor.cpp +++ b/tests/cpp_unit_tests/test_three_phase_tensor.cpp @@ -206,6 +206,38 @@ TEST_CASE("Three phase tensor") { CHECK(is_nan(vec(1))); CHECK(vec(2) == 6.0); } + SUBCASE("phase_mod_2pi") { + auto const check = [](double value) { + CAPTURE(value); + CHECK_GE(value, -pi); + CHECK_LE(value, pi); + if (value != pi && value != -pi) { + CHECK(phase_mod_2pi(value) == doctest::Approx(value)); + } + }; + auto const check_asym = [&check](RealValue const& value) { + for (Idx i : {0, 1, 2}) { + CAPTURE(i); + check(value(i)); + } + }; + check(phase_mod_2pi(0.0)); + check(phase_mod_2pi(2.0 * pi)); + check(phase_mod_2pi(2.0 * pi + 1.0)); + check(phase_mod_2pi(-1.0)); + check(phase_mod_2pi(-pi)); + check(phase_mod_2pi(pi)); + check(phase_mod_2pi(-3.0 * pi)); + check(phase_mod_2pi(3.0 * pi)); + check(phase_mod_2pi(pi * (1.0 + std::numeric_limits::epsilon()))); + check(phase_mod_2pi(pi * (1.0 - std::numeric_limits::epsilon()))); + check(phase_mod_2pi(-pi * (1.0 + std::numeric_limits::epsilon()))); + check(phase_mod_2pi(-pi * (1.0 - std::numeric_limits::epsilon()))); + + check_asym(phase_mod_2pi(RealValue{0.0, 2.0 * pi, 2.0 * pi + 1.0})); + check_asym(phase_mod_2pi(RealValue{0.0, 2.0 * pi, 2.0 * pi + 1.0})); + check_asym(phase_mod_2pi(RealValue{0.0, 2.0 * pi, 2.0 * pi + 1.0})); + } } } // namespace power_grid_model \ No newline at end of file