From e606f4c4d935e12558ace9e4d598c01ea4fd2f3e Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Wed, 29 Jan 2025 12:31:46 +0100 Subject: [PATCH 01/12] minor refactor Signed-off-by: Martijn Govers --- .../component/current_sensor.hpp | 28 ++++++++----------- .../component/power_sensor.hpp | 6 ++-- 2 files changed, 14 insertions(+), 20 deletions(-) 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 cbfd1ac8c8..b32189e30d 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 @@ -73,12 +73,12 @@ template class CurrentSensor : public Ge explicit CurrentSensor(CurrentSensorInput 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 +93,13 @@ template class CurrentSensor : public Ge }; UpdateChange update(CurrentSensorUpdate 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(update_data.i_sigma, i_sigma_, base_current_inv_); + update_real_value(update_data.i_angle_sigma, i_angle_sigma_, 1.0); update_real_value(update_data.i_measured, i_measured_, base_current_inv_); update_real_value(update_data.i_angle_measured, i_angle_measured_, 1.0); + return {false, false}; } @@ -117,17 +116,12 @@ template class CurrentSensor : public Ge } private: + double base_current_{}; + double base_current_inv_{}; RealValue i_measured_{}; RealValue i_angle_measured_{}; double i_sigma_{}; double i_angle_sigma_{}; - double base_current_{}; - double base_current_inv_{}; - - void set_current(CurrentSensorInput 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 diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/component/power_sensor.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/component/power_sensor.hpp index 83aed9efb9..46ca10f1b7 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/component/power_sensor.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/component/power_sensor.hpp @@ -87,11 +87,11 @@ template class PowerSensor : public Generi }; UpdateChange update(PowerSensorUpdate 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(update_data.power_sigma, apparent_power_sigma_, inv_base_power); update_real_value(update_data.p_sigma, p_sigma_, inv_base_power); update_real_value(update_data.q_sigma, q_sigma_, inv_base_power); From 698bc1f3494b620bf614a0dddec24e852b6aa28d Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Wed, 29 Jan 2025 12:32:20 +0100 Subject: [PATCH 02/12] add calculation parameter tests Signed-off-by: Martijn Govers --- tests/cpp_unit_tests/test_current_sensor.cpp | 28 ++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/tests/cpp_unit_tests/test_current_sensor.cpp b/tests/cpp_unit_tests/test_current_sensor.cpp index 7b8d0b1a39..66384863cf 100644 --- a/tests/cpp_unit_tests/test_current_sensor.cpp +++ b/tests/cpp_unit_tests/test_current_sensor.cpp @@ -94,6 +94,34 @@ TEST_CASE("Test current sensor") { InvalidMeasuredTerminalType); } } + SUBCASE("Symmetric calculation parameters") { + CurrentSensor sym_current_sensor{ + {1, 1, MeasuredTerminalType::branch3_1, AngleMeasurementType::local}, 1.0}; + + 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(); + + CHECK(sym_param.angle_measurement_type == AngleMeasurementType::local); + CHECK(sym_param.i_real_variance == doctest::Approx(1.0)); + CHECK(sym_param.i_imag_variance == doctest::Approx(0.2)); + CHECK(real(sym_param.value) == doctest::Approx(1.0)); + CHECK(imag(sym_param.value) == doctest::Approx(0.0)); + } + + 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(); + + CHECK(sym_param.angle_measurement_type == AngleMeasurementType::local); + CHECK(sym_param.i_real_variance == doctest::Approx(0.2)); + CHECK(sym_param.i_imag_variance == doctest::Approx(1.0)); + CHECK(real(sym_param.value) == doctest::Approx(0.0)); + CHECK(imag(sym_param.value) == doctest::Approx(1.0)); + } + } } SUBCASE("Update inverse - sym") { constexpr auto i_measured = 1.0; From dab85bcc2aaf9bd5194572058893045143aed94a Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Tue, 4 Feb 2025 08:58:21 +0100 Subject: [PATCH 03/12] add sym calculation params for symmetric current sensor Signed-off-by: Martijn Govers --- .../calculation_parameters.hpp | 12 +-- .../power_grid_model/common/statistics.hpp | 96 +++++++++++++++++++ .../component/current_sensor.hpp | 24 ++++- .../power_grid_model_c/src/handle.hpp | 4 +- tests/cpp_unit_tests/test_current_sensor.cpp | 46 ++++++--- 5 files changed, 151 insertions(+), 31 deletions(-) create mode 100644 power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp index 8e70e3f113..ffb2fcf991 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp +++ b/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,17 +79,6 @@ template struct ApplianceShortCircuitSolverOutput { ComplexValue 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 struct UniformComplexRandomVariable { - using sym = sym_type; - - static constexpr bool symmetric{is_symmetric_v}; - - ComplexValue 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 diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp new file mode 100644 index 0000000000..3599ca82cd --- /dev/null +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp @@ -0,0 +1,96 @@ +// SPDX-FileCopyrightText: Contributors to the Power Grid Model project +// +// SPDX-License-Identifier: MPL-2.0 + +#pragma once + +#include "common.hpp" +#include "three_phase_tensor.hpp" + +namespace power_grid_model { +template struct UniformRealRandomVariable { + using sym = sym_type; + + static constexpr bool symmetric{is_symmetric_v}; + + RealValue value{}; + double variance{}; // variance (sigma^2) of the error range, in p.u. +}; + +template struct IndependentRealRandomVariable { + using sym = sym_type; + + static constexpr bool symmetric{is_symmetric_v}; + + RealValue value{}; + RealValue variance{}; // variance (sigma^2) of the error range, in p.u. +}; + +// 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 struct UniformComplexRandomVariable { + using sym = sym_type; + + static constexpr bool symmetric{is_symmetric_v}; + + ComplexValue value{}; + double variance{}; // variance (sigma^2) of the error range, in p.u. +}; + +// Complex measured value of a sensor in p.u. modeled as separate real and imaginary components with uniform variances +// (circularly symmetric) +template struct DecomposedUniformComplexRandomVariable { + using sym = sym_type; + + static constexpr bool symmetric{is_symmetric_v}; + + UniformRealRandomVariable real_component; + UniformRealRandomVariable imag_component; + + ComplexValue value() const { return {real_component.value, imag_component.value}; } +}; + +// Complex measured value of a sensor in p.u. modeled as separate real and imaginary components with independent +// variances (circularly symmetric) +template struct DecomposedIndependentComplexRandomVariable { + using sym = sym_type; + + static constexpr bool symmetric{is_symmetric_v}; + + IndependentRealRandomVariable real_component; + IndependentRealRandomVariable imag_component; + + ComplexValue value() const { return {real_component.value, imag_component.value}; } +}; + +// Complex measured value of a sensor in p.u. in polar coordinates (magnitude and angle) +// (circularly symmetric) +template struct PolarComplexRandomVariable { + using sym = sym_type; + + static constexpr bool symmetric{is_symmetric_v}; + + UniformRealRandomVariable magnitude; + UniformRealRandomVariable angle; + + ComplexValue value() const { return magnitude.value * exp(1.0i * angle.value); } + + explicit operator UniformComplexRandomVariable() const { + return UniformComplexRandomVariable{ + .value = value(), .variance = magnitude.variance + magnitude.value * magnitude.value * angle.variance}; + } + explicit operator DecomposedIndependentComplexRandomVariable() 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 DecomposedIndependentComplexRandomVariable{ + .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}}; + } +}; +} // namespace power_grid_model 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 b32189e30d..933c698595 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 @@ -11,6 +11,8 @@ #include "../common/enum.hpp" #include "../common/exception.hpp" +#include + namespace power_grid_model { class GenericCurrentSensor : public Sensor { @@ -48,6 +50,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_; @@ -123,12 +129,22 @@ template class CurrentSensor : public Ge double i_sigma_{}; double i_angle_sigma_{}; - // TODO when filling the functions below take in mind that i_angle_sigma is optional + // TODO(mgovers) when filling the functions below take in mind that i_angle_sigma is optional CurrentSensorCalcParam sym_calc_param() const final { - CurrentSensorCalcParam calc_param{}; - // TODO - return calc_param; + if constexpr (is_asymmetric_v) { + return {}; // TODO + } else { + auto const i_polar = PolarComplexRandomVariable{ + .magnitude = {.value = i_measured_, .variance = i_sigma_ * i_sigma_}, + .angle = {.value = i_angle_measured_, .variance = i_angle_sigma_ * i_angle_sigma_}}; + auto const i_decomposed = static_cast>(i_polar); + + return CurrentSensorCalcParam{.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}; + } } CurrentSensorCalcParam asym_calc_param() const final { CurrentSensorCalcParam calc_param{}; diff --git a/power_grid_model_c/power_grid_model_c/src/handle.hpp b/power_grid_model_c/power_grid_model_c/src/handle.hpp index 29a2a6ced9..54c58a6c46 100644 --- a/power_grid_model_c/power_grid_model_c/src/handle.hpp +++ b/power_grid_model_c/power_grid_model_c/src/handle.hpp @@ -28,8 +28,8 @@ struct PGM_Handle { }; template -auto call_with_catch(PGM_Handle* handle, Functor func, PGM_Idx error_code, - std::string_view extra_msg = {}) -> std::invoke_result_t { +inline auto call_with_catch(PGM_Handle* handle, Functor func, PGM_Idx error_code, + std::string_view extra_msg = {}) -> std::invoke_result_t { if (handle) { PGM_clear_error(handle); } diff --git a/tests/cpp_unit_tests/test_current_sensor.cpp b/tests/cpp_unit_tests/test_current_sensor.cpp index 66384863cf..e4f304db4b 100644 --- a/tests/cpp_unit_tests/test_current_sensor.cpp +++ b/tests/cpp_unit_tests/test_current_sensor.cpp @@ -40,7 +40,7 @@ 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; @@ -60,10 +60,10 @@ 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(imag(sym_sensor_param.value) == doctest::Approx(0.0)); + CHECK(sym_sensor_param.i_real_variance == doctest::Approx(pow(1.0 / base_current, 2))); + CHECK(sym_sensor_param.i_imag_variance == doctest::Approx(pow(0.2 * 1.0e3 / base_current, 2))); + CHECK(real(sym_sensor_param.value) == doctest::Approx(1.0e3 / base_current)); + CHECK(imag(sym_sensor_param.value) == doctest::Approx(0.0 / base_current)); CHECK(is_nan(sym_sensor_output.id)); CHECK(is_nan(sym_sensor_output.energized)); @@ -95,8 +95,11 @@ TEST_CASE("Test current sensor") { } } SUBCASE("Symmetric calculation parameters") { + double const u_rated = 10.0e3; + double const base_current = base_power_3p / u_rated / sqrt3; + CurrentSensor sym_current_sensor{ - {1, 1, MeasuredTerminalType::branch3_1, AngleMeasurementType::local}, 1.0}; + {1, 1, MeasuredTerminalType::branch3_1, AngleMeasurementType::local}, u_rated}; SUBCASE("No phase shift") { sym_current_sensor.update( @@ -104,10 +107,10 @@ TEST_CASE("Test current sensor") { auto const sym_param = sym_current_sensor.calc_param(); CHECK(sym_param.angle_measurement_type == AngleMeasurementType::local); - CHECK(sym_param.i_real_variance == doctest::Approx(1.0)); - CHECK(sym_param.i_imag_variance == doctest::Approx(0.2)); - CHECK(real(sym_param.value) == doctest::Approx(1.0)); - CHECK(imag(sym_param.value) == doctest::Approx(0.0)); + 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") { @@ -116,10 +119,25 @@ TEST_CASE("Test current sensor") { auto const sym_param = sym_current_sensor.calc_param(); CHECK(sym_param.angle_measurement_type == AngleMeasurementType::local); - CHECK(sym_param.i_real_variance == doctest::Approx(0.2)); - CHECK(sym_param.i_imag_variance == doctest::Approx(1.0)); - CHECK(real(sym_param.value) == doctest::Approx(0.0)); - CHECK(imag(sym_param.value) == doctest::Approx(1.0)); + 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(); + + 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))); } } } From 4f63f86de56a2589e83c4bd90931f781eb23ed81 Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Tue, 4 Feb 2025 10:35:10 +0100 Subject: [PATCH 04/12] unit tests for statistics Signed-off-by: Martijn Govers --- tests/cpp_unit_tests/CMakeLists.txt | 1 + tests/cpp_unit_tests/test_statistics.cpp | 155 +++++++++++++++++++++++ 2 files changed, 156 insertions(+) create mode 100644 tests/cpp_unit_tests/test_statistics.cpp diff --git a/tests/cpp_unit_tests/CMakeLists.txt b/tests/cpp_unit_tests/CMakeLists.txt index 12c4d801a5..9bcf752167 100644 --- a/tests/cpp_unit_tests/CMakeLists.txt +++ b/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" diff --git a/tests/cpp_unit_tests/test_statistics.cpp b/tests/cpp_unit_tests/test_statistics.cpp new file mode 100644 index 0000000000..e45b273edb --- /dev/null +++ b/tests/cpp_unit_tests/test_statistics.cpp @@ -0,0 +1,155 @@ +// SPDX-FileCopyrightText: Contributors to the Power Grid Model project +// +// SPDX-License-Identifier: MPL-2.0 + +#include + +#include + +namespace power_grid_model { +namespace { +using std::numbers::sqrt2; +constexpr auto inv_sqrt2 = sqrt2 / 2; +} // namespace + +TEST_CASE("Test statistics") { + SUBCASE("PolarComplexRandomVariable") { + SUBCASE("Constructor") { + for (auto const [magnitude, magnitude_variance, angle, angle_variance] : + std::array{std::tuple{1.0, 1.0, 0.0, 0.2}, std::tuple{2.0, 3.0, 0.0, 0.2}, + std::tuple{1.0, 1.0, pi / 2, 0.2}, std::tuple{1.0, 1.0, pi / 4, 0.2}}) { + CAPTURE(magnitude); + CAPTURE(magnitude_variance); + CAPTURE(angle); + CAPTURE(angle_variance); + + PolarComplexRandomVariable const polar{ + .magnitude = {.value = magnitude, .variance = magnitude_variance}, + .angle = {.value = angle, .variance = angle_variance}}; + + CHECK(polar.magnitude.value == magnitude); + CHECK(polar.magnitude.variance == magnitude_variance); + CHECK(polar.angle.value == angle); + CHECK(polar.angle.variance == angle_variance); + } + } + SUBCASE("Aggregate value") { + for (auto const [magnitude, magnitude_variance, angle_variance] : + std::array{std::tuple{1.0, 1.0, 0.2}, std::tuple{2.0, 1.0, 0.2}, std::tuple{1.0, 3.0, 0.2}, + std::tuple{1.0, 2.0, 0.4}}) { + CAPTURE(magnitude); + CAPTURE(magnitude_variance); + CAPTURE(angle_variance); + + SUBCASE("No phase shift") { + PolarComplexRandomVariable const polar{ + .magnitude = {.value = magnitude, .variance = magnitude_variance}, + .angle = {.value = 0.0, .variance = angle_variance}}; + + CHECK(real(polar.value()) == doctest::Approx(polar.magnitude.value)); + CHECK(imag(polar.value()) == doctest::Approx(0.0)); + } + SUBCASE("Perpendicular phase shift") { + PolarComplexRandomVariable const polar{ + .magnitude = {.value = magnitude, .variance = magnitude_variance}, + .angle = {.value = pi / 2, .variance = angle_variance}}; + + CHECK(real(polar.value()) == doctest::Approx(0.0)); + CHECK(imag(polar.value()) == doctest::Approx(polar.magnitude.value)); + } + SUBCASE("45deg phase shift") { + PolarComplexRandomVariable const polar{ + .magnitude = {.value = magnitude, .variance = magnitude_variance}, + .angle = {.value = pi / 4, .variance = angle_variance}}; + + CHECK(real(polar.value()) == polar.magnitude.value * inv_sqrt2); + CHECK(imag(polar.value()) == real(polar.value())); + } + } + } + + SUBCASE("Conversion to UniformComplexRandomVariable") { + for (auto const [magnitude, magnitude_variance, angle, angle_variance] : + std::array{std::tuple{1.0, 1.0, 0.0, 0.2}, std::tuple{2.0, 3.0, 0.0, 0.2}, + std::tuple{1.0, 1.0, pi / 2, 0.2}, std::tuple{1.0, 1.0, pi / 4, 0.2}}) { + CAPTURE(magnitude); + CAPTURE(magnitude_variance); + CAPTURE(angle); + CAPTURE(angle_variance); + + PolarComplexRandomVariable const polar{ + .magnitude = {.value = magnitude, .variance = magnitude_variance}, + .angle = {.value = angle, .variance = angle_variance}}; + + auto const uniform = static_cast>(polar); + + CHECK(real(uniform.value) == doctest::Approx(real(polar.value()))); + CHECK(imag(uniform.value) == doctest::Approx(imag(polar.value()))); + CHECK(uniform.variance == + doctest::Approx(polar.magnitude.variance + magnitude * magnitude * polar.angle.variance)); + } + } + + SUBCASE("Conversion to DecomposedIndependentComplexRandomVariable") { + for (auto const [magnitude, magnitude_variance, angle_variance] : + std::array{std::tuple{1.0, 1.0, 0.2}, std::tuple{2.0, 1.0, 0.2}, std::tuple{1.0, 3.0, 0.2}, + std::tuple{1.0, 2.0, 0.4}}) { + CAPTURE(magnitude); + CAPTURE(magnitude_variance); + CAPTURE(angle_variance); + + SUBCASE("No phase shift") { + PolarComplexRandomVariable const polar{ + .magnitude = {.value = magnitude, .variance = magnitude}, + .angle = {.value = 0.0, .variance = angle_variance}}; + + auto const independent = + static_cast>(polar); + + CHECK(independent.real_component.value == doctest::Approx(polar.magnitude.value)); + CHECK(independent.imag_component.value == doctest::Approx(0.0)); + CHECK(independent.real_component.variance == doctest::Approx(polar.magnitude.variance)); + CHECK(independent.imag_component.variance == + doctest::Approx(magnitude * magnitude * polar.angle.variance)); + CHECK(real(independent.value()) == doctest::Approx(real(polar.value()))); + CHECK(imag(independent.value()) == doctest::Approx(imag(polar.value()))); + } + + SUBCASE("Perpendicular phase shift") { + PolarComplexRandomVariable const polar{ + .magnitude = {.value = magnitude, .variance = magnitude}, + .angle = {.value = pi / 2, .variance = angle_variance}}; + + auto const independent = + static_cast>(polar); + + CHECK(independent.real_component.value == doctest::Approx(0.0)); + CHECK(independent.imag_component.value == doctest::Approx(polar.magnitude.value)); + CHECK(independent.real_component.variance == + doctest::Approx(magnitude * magnitude * polar.angle.variance)); + CHECK(independent.imag_component.variance == doctest::Approx(polar.magnitude.variance)); + CHECK(real(independent.value()) == doctest::Approx(real(polar.value()))); + CHECK(imag(independent.value()) == doctest::Approx(imag(polar.value()))); + } + + SUBCASE("45deg phase shift") { + PolarComplexRandomVariable const polar{ + .magnitude = {.value = magnitude, .variance = magnitude}, + .angle = {.value = pi / 4, .variance = angle_variance}}; + + auto const independent = + static_cast>(polar); + auto const uniform = static_cast>(polar); + + CHECK(independent.real_component.value == doctest::Approx(real(uniform.value))); + CHECK(independent.imag_component.value == doctest::Approx(imag(uniform.value))); + CHECK(independent.real_component.variance == doctest::Approx(uniform.variance / 2)); + CHECK(independent.imag_component.variance == doctest::Approx(independent.real_component.variance)); + CHECK(real(independent.value()) == doctest::Approx(real(polar.value()))); + CHECK(imag(independent.value()) == doctest::Approx(imag(polar.value()))); + } + } + } + } +} +} // namespace power_grid_model From d5234f154965476a10a49ad0d31c4f47efd10d02 Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Tue, 4 Feb 2025 15:46:21 +0100 Subject: [PATCH 05/12] extend statistics functionality Signed-off-by: Martijn Govers --- .../calculation_parameters.hpp | 4 +- .../power_grid_model/common/statistics.hpp | 78 +++- .../component/current_sensor.hpp | 2 +- tests/cpp_unit_tests/test_current_sensor.cpp | 4 +- tests/cpp_unit_tests/test_statistics.cpp | 361 +++++++++++++++--- 5 files changed, 384 insertions(+), 65 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp index ffb2fcf991..b49718e52c 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp @@ -108,8 +108,8 @@ template struct CurrentSensorCalcParam { AngleMeasurementType angle_measurement_type{}; ComplexValue 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 i_real_variance{}; // variance (sigma^2) of the error range of real part of the current, in p.u. + RealValue i_imag_variance{}; // variance (sigma^2) of the error range of imaginary part of the current, in p.u. }; template diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp index 3599ca82cd..95819525e2 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp @@ -15,6 +15,17 @@ template struct UniformRealRandomVariable { RealValue value{}; double variance{}; // variance (sigma^2) of the error range, in p.u. + + explicit operator UniformRealRandomVariable() const + requires(is_symmetric_v) + { + return {.value = RealValue{std::piecewise_construct, value}, .variance = variance}; + } + explicit operator UniformRealRandomVariable() const + requires(is_asymmetric_v) + { + return {.value = mean_val(value), .variance = variance / 3.0}; + } }; template struct IndependentRealRandomVariable { @@ -24,10 +35,29 @@ template struct IndependentRealRandomVariable { RealValue value{}; RealValue variance{}; // variance (sigma^2) of the error range, in p.u. + + explicit operator UniformRealRandomVariable() const { + constexpr auto scale = is_asymmetric_v ? 3.0 : 1.0; + return {.value = mean_val(value), .variance = mean_val(variance) / scale}; + } + explicit operator UniformRealRandomVariable() const { + return {.value = value, .variance = mean_val(variance)}; + } + explicit operator IndependentRealRandomVariable() const + requires(is_symmetric_v) + { + return {.value = RealValue{std::piecewise_construct, value}, + .variance = RealValue{std::piecewise_construct, variance}}; + } + explicit operator IndependentRealRandomVariable() const + requires(is_asymmetric_v) + { + return {.value = mean_val(value), .variance = mean_val(variance) / 3.0}; + } }; // Complex measured value of a sensor in p.u. with a uniform variance across all phases and axes of the complex plane -// (circularly symmetric) +// (rotationally symmetric) template struct UniformComplexRandomVariable { using sym = sym_type; @@ -37,22 +67,31 @@ template struct UniformComplexRandomVariable { double variance{}; // variance (sigma^2) of the error range, in p.u. }; -// Complex measured value of a sensor in p.u. modeled as separate real and imaginary components with uniform variances -// (circularly symmetric) -template struct DecomposedUniformComplexRandomVariable { +inline UniformComplexRandomVariable pos_seq(UniformComplexRandomVariable const& var) { + return {.value = pos_seq(var.value), .variance = var.variance / 3.0}; +} +inline UniformComplexRandomVariable three_phase(UniformComplexRandomVariable const& var) { + return {.value = ComplexValue{var.value}, .variance = var.variance}; +} + +// Complex measured value of a sensor in p.u. with separate variances per phase (but rotationally symmetric in the +// complex plane) +template struct IndependentComplexRandomVariable { using sym = sym_type; static constexpr bool symmetric{is_symmetric_v}; - UniformRealRandomVariable real_component; - UniformRealRandomVariable imag_component; + ComplexValue value{}; + RealValue variance{}; // variance (sigma^2) of the error range, in p.u. - ComplexValue value() const { return {real_component.value, imag_component.value}; } + explicit operator UniformComplexRandomVariable() const { + return UniformComplexRandomVariable{.value = value, .variance = sum_val(variance)}; + } }; // Complex measured value of a sensor in p.u. modeled as separate real and imaginary components with independent -// variances (circularly symmetric) -template struct DecomposedIndependentComplexRandomVariable { +// variances (rotationally symmetric) +template struct DecomposedComplexRandomVariable { using sym = sym_type; static constexpr bool symmetric{is_symmetric_v}; @@ -61,10 +100,19 @@ template struct DecomposedIndependentComplexRandomVariab IndependentRealRandomVariable imag_component; ComplexValue value() const { return {real_component.value, imag_component.value}; } + + explicit operator UniformComplexRandomVariable() const { + return static_cast>( + static_cast>(*this)); + } + explicit operator IndependentComplexRandomVariable() const { + return IndependentComplexRandomVariable{.value = value(), + .variance = real_component.variance + imag_component.variance}; + } }; // Complex measured value of a sensor in p.u. in polar coordinates (magnitude and angle) -// (circularly symmetric) +// (rotationally symmetric) template struct PolarComplexRandomVariable { using sym = sym_type; @@ -76,15 +124,19 @@ template struct PolarComplexRandomVariable { ComplexValue value() const { return magnitude.value * exp(1.0i * angle.value); } explicit operator UniformComplexRandomVariable() const { - return UniformComplexRandomVariable{ + return static_cast>( + static_cast>(*this)); + } + explicit operator IndependentComplexRandomVariable() const { + return IndependentComplexRandomVariable{ .value = value(), .variance = magnitude.variance + magnitude.value * magnitude.value * angle.variance}; } - explicit operator DecomposedIndependentComplexRandomVariable() const { + explicit operator DecomposedComplexRandomVariable() 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 DecomposedIndependentComplexRandomVariable{ + return DecomposedComplexRandomVariable{ .real_component = {.value = real_component, .variance = magnitude.variance * cos_theta * cos_theta + imag_component * imag_component * angle.variance}, 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 933c698595..eabc1948e9 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 @@ -138,7 +138,7 @@ template class CurrentSensor : public Ge auto const i_polar = PolarComplexRandomVariable{ .magnitude = {.value = i_measured_, .variance = i_sigma_ * i_sigma_}, .angle = {.value = i_angle_measured_, .variance = i_angle_sigma_ * i_angle_sigma_}}; - auto const i_decomposed = static_cast>(i_polar); + auto const i_decomposed = static_cast>(i_polar); return CurrentSensorCalcParam{.angle_measurement_type = angle_measurement_type(), .value = i_decomposed.value(), diff --git a/tests/cpp_unit_tests/test_current_sensor.cpp b/tests/cpp_unit_tests/test_current_sensor.cpp index e4f304db4b..7aa168c940 100644 --- a/tests/cpp_unit_tests/test_current_sensor.cpp +++ b/tests/cpp_unit_tests/test_current_sensor.cpp @@ -71,8 +71,8 @@ TEST_CASE("Test current sensor") { CHECK(is_nan(sym_sensor_output.i_angle_residual)); // 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(asym_sensor_param.i_real_variance[0] == doctest::Approx(0.0)); + CHECK(asym_sensor_param.i_imag_variance[1] == 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)); diff --git a/tests/cpp_unit_tests/test_statistics.cpp b/tests/cpp_unit_tests/test_statistics.cpp index e45b273edb..35081384d0 100644 --- a/tests/cpp_unit_tests/test_statistics.cpp +++ b/tests/cpp_unit_tests/test_statistics.cpp @@ -8,11 +8,259 @@ namespace power_grid_model { namespace { +using AsymRealValue = RealValue; + using std::numbers::sqrt2; constexpr auto inv_sqrt2 = sqrt2 / 2; } // namespace TEST_CASE("Test statistics") { + SUBCASE("UniformRealRandomVariable") { + for (auto const [value, variance] : + std::array{std::tuple{1.0, 1.0}, std::tuple{2.0, 3.0}, std::tuple{0.0, 1.0}, std::tuple{2.0, 3.0}}) { + CAPTURE(value); + CAPTURE(variance); + + UniformRealRandomVariable const uniform{.value = value, .variance = variance}; + + SUBCASE("Constructor") { + CHECK(uniform.value == value); + CHECK(uniform.variance == variance); + } + SUBCASE("Conversion to UniformRealRandomVariable") { + auto const asymmetric = static_cast>(uniform); + + CHECK(asymmetric.value(0) == doctest::Approx(uniform.value)); + CHECK(asymmetric.value(1) == doctest::Approx(uniform.value)); + CHECK(asymmetric.value(2) == doctest::Approx(uniform.value)); + CHECK(asymmetric.variance == doctest::Approx(variance)); + } + } + } + SUBCASE("UniformRealRandomVariable") { + for (auto const [value_a, value_b, value_c, variance] : + std::array{std::tuple{1.0, 2.0, 3.0, 1.0}, std::tuple{2.0, 2.1, 2.2, 3.0}, std::tuple{0.0, 0.1, 0.2, 1.0}, + std::tuple{2.0, 0.0, 0.0, 3.0}}) { + CAPTURE(value_a); + CAPTURE(value_b); + CAPTURE(value_c); + CAPTURE(variance); + + UniformRealRandomVariable const uniform{.value = {value_a, value_b, value_c}, + .variance = variance}; + + SUBCASE("Constructor") { + CHECK(uniform.value(0) == value_a); + CHECK(uniform.value(1) == value_b); + CHECK(uniform.value(2) == value_c); + CHECK(uniform.variance == variance); + } + SUBCASE("Conversion to UniformRealRandomVariable") { + auto const symmetric = static_cast>(uniform); + + CHECK(symmetric.value == doctest::Approx(mean_val(uniform.value))); + CHECK(symmetric.variance == doctest::Approx(variance / 3)); + } + } + } + + SUBCASE("IndependentRealRandomVariable") { + for (auto const [value, variance] : + std::array{std::tuple{1.0, 1.0}, std::tuple{2.0, 3.0}, std::tuple{0.0, 1.0}, std::tuple{2.0, 3.0}}) { + CAPTURE(value); + CAPTURE(variance); + + IndependentRealRandomVariable const independent{.value = value, .variance = variance}; + + SUBCASE("Constructor") { + CHECK(independent.value == value); + CHECK(independent.variance == variance); + } + SUBCASE("Conversion to IndependentRealRandomVariable") { + auto const asymmetric = static_cast>(independent); + auto const symmetric = static_cast>(asymmetric); + + CHECK(asymmetric.value(0) == doctest::Approx(independent.value)); + CHECK(asymmetric.value(1) == doctest::Approx(independent.value)); + CHECK(asymmetric.value(2) == doctest::Approx(independent.value)); + CHECK(asymmetric.variance(0) == doctest::Approx(independent.variance)); + CHECK(asymmetric.variance(1) == doctest::Approx(independent.variance)); + CHECK(asymmetric.variance(2) == doctest::Approx(independent.variance)); + } + } + } + SUBCASE("IndependentRealRandomVariable") { + for (auto const [value_a, value_b, value_c, variance_a, variance_b, variance_c] : + std::array{std::tuple{1.0, 2.0, 3.0, 1.0, 2.0, 3.0}, std::tuple{2.0, 2.1, 2.2, 3.0, 1.0, 2.0}, + std::tuple{0.0, 0.1, 0.2, 1.0, 1.0, 1.0}, std::tuple{2.0, 0.0, 0.0, 3.0, 3.0, 3.0}}) { + CAPTURE(value_a); + CAPTURE(value_b); + CAPTURE(value_c); + CAPTURE(variance_a); + CAPTURE(variance_b); + CAPTURE(variance_c); + + IndependentRealRandomVariable const independent{ + .value = {value_a, value_b, value_c}, .variance = {variance_a, variance_b, variance_c}}; + + SUBCASE("Constructor") { + CHECK(independent.value(0) == value_a); + CHECK(independent.value(1) == value_b); + CHECK(independent.value(2) == value_c); + CHECK(independent.variance(0) == variance_a); + CHECK(independent.variance(1) == variance_b); + CHECK(independent.variance(2) == variance_c); + } + SUBCASE("Conversion to IndependentRealRandomVariable") { + auto const symmetric = static_cast>(independent); + auto const asymmetric = static_cast>(symmetric); + + CHECK(symmetric.value == doctest::Approx(mean_val(independent.value))); + CHECK(symmetric.variance == doctest::Approx(mean_val(independent.variance) / 3.0)); + } + SUBCASE("Conversion to UniformRealRandomVariable") { + auto const uniform = static_cast>(independent); + auto const via_asym_uniform = static_cast>( + static_cast>(independent)); + auto const via_sym_independent = static_cast>( + static_cast>(independent)); + + CHECK(uniform.value == doctest::Approx(via_asym_uniform.value)); + CHECK(uniform.variance == doctest::Approx(via_asym_uniform.variance)); + CHECK(uniform.value == doctest::Approx(via_sym_independent.value)); + CHECK(uniform.variance == doctest::Approx(via_sym_independent.variance)); + } + } + } + + SUBCASE("UniformComplexRandomVariable") { + for (auto const [real_value, imag_value, variance] : + std::array{std::tuple{1.0, 0.0, 1.0}, std::tuple{2.0, 0.0, 3.0}, std::tuple{0.0, 1.0, 1.0}, + std::tuple{0.0, 2.0, 1.0}, std::tuple{1.0, 1.0, 1.0}, std::tuple{2.0, 2.0, 3.0}}) { + CAPTURE(real_value); + CAPTURE(imag_value); + CAPTURE(variance); + + UniformComplexRandomVariable const uniform{.value = {real_value, imag_value}, + .variance = variance}; + + SUBCASE("Constructor") { + CHECK(real(uniform.value) == real_value); + CHECK(imag(uniform.value) == imag_value); + CHECK(uniform.variance == variance); + } + SUBCASE("To three-phase") { + auto const asymmetric = three_phase(uniform); + + CHECK(real(asymmetric.value(0)) == doctest::Approx(real(uniform.value))); + CHECK(imag(asymmetric.value(0)) == doctest::Approx(imag(uniform.value))); + CHECK(real(asymmetric.value(1)) == doctest::Approx(real(uniform.value * a2))); + CHECK(imag(asymmetric.value(1)) == doctest::Approx(imag(uniform.value * a2))); + CHECK(real(asymmetric.value(2)) == doctest::Approx(real(uniform.value * a))); + CHECK(imag(asymmetric.value(2)) == doctest::Approx(imag(uniform.value * a))); + CHECK(asymmetric.variance == variance); + } + } + } + + SUBCASE("UniformComplexRandomVariable") { + for (auto const [real_value, imag_value, variance] : + std::array{std::tuple{AsymRealValue{1.0, 2.0, 3.0}, AsymRealValue{0.0, 0.0, 0.0}, 1.0}, + std::tuple{AsymRealValue{2.0, 0.0, 0.0}, AsymRealValue{0.0, 3.0, 3.0}, 3.0}, + std::tuple{AsymRealValue{0.0, 0.0, 0.0}, AsymRealValue{1.0, 1.0, 1.0}, 1.0}, + std::tuple{AsymRealValue{0.0, -1.0, -2.0}, AsymRealValue{2.0, -1.0, -2.0}, 1.0}, + std::tuple{AsymRealValue{1.0, 1.0, 1.0}, AsymRealValue{1.0, 1.0, 1.0}, 1.0}, + std::tuple{AsymRealValue{2.0, 2.0, 2.0}, AsymRealValue{2.0, 2.0, 2.0}, 3.0}}) { + CAPTURE(real_value); + CAPTURE(imag_value); + CAPTURE(variance); + + UniformComplexRandomVariable const uniform{.value = {real_value, imag_value}, + .variance = variance}; + + SUBCASE("Constructor") { + CHECK(real(uniform.value(0)) == real_value(0)); + CHECK(imag(uniform.value(0)) == imag_value(0)); + CHECK(real(uniform.value(1)) == real_value(1)); + CHECK(imag(uniform.value(1)) == imag_value(1)); + CHECK(real(uniform.value(2)) == real_value(2)); + CHECK(imag(uniform.value(2)) == imag_value(2)); + CHECK(uniform.variance == variance); + } + SUBCASE("Positive sequence") { + auto const positive_sequence = pos_seq(uniform); + + CHECK(real(positive_sequence.value) == doctest::Approx(real(pos_seq(uniform.value)))); + CHECK(imag(positive_sequence.value) == doctest::Approx(imag(pos_seq(uniform.value)))); + CHECK(positive_sequence.variance == doctest::Approx(variance / 3.0)); + } + } + } + SUBCASE("IndependentComplexRandomVariable") { + for (auto const [real_value, imag_value, variance] : + std::array{std::tuple{1.0, 0.0, 1.0}, std::tuple{2.0, 0.0, 3.0}, std::tuple{0.0, 1.0, 1.0}, + std::tuple{0.0, 2.0, 1.0}, std::tuple{1.0, 1.0, 1.0}, std::tuple{2.0, 2.0, 3.0}}) { + CAPTURE(real_value); + CAPTURE(imag_value); + CAPTURE(variance); + + IndependentComplexRandomVariable const independent{.value = {real_value, imag_value}, + .variance = variance}; + + SUBCASE("Constructor") { + CHECK(real(independent.value) == real_value); + CHECK(imag(independent.value) == imag_value); + CHECK(independent.variance == variance); + } + SUBCASE("Conversion to UniformComplexRandomVariable") { + auto const uniform = static_cast>(independent); + + CHECK(real(uniform.value) == doctest::Approx(real(independent.value))); + CHECK(imag(uniform.value) == doctest::Approx(imag(independent.value))); + CHECK(uniform.variance == doctest::Approx(variance)); + } + } + } + SUBCASE("DecomposedComplexRandomVariable") { + for (auto const [real_value, real_variance, imag_value, imag_variance] : std::array{ + std::tuple{1.0, 1.0, 0.0, 0.2}, std::tuple{2.0, 3.0, 0.0, 0.2}, std::tuple{0.0, 1.0, 1.0, 0.2}, + std::tuple{0.0, 1.0, 2.0, 0.2}, std::tuple{1.0, 1.0, 1.0, 0.2}, std::tuple{2.0, 1.0, 2.0, 0.2}}) { + CAPTURE(real_value); + CAPTURE(real_variance); + CAPTURE(imag_value); + CAPTURE(imag_variance); + + DecomposedComplexRandomVariable const decomposed{ + .real_component = {.value = real_value, .variance = real_variance}, + .imag_component = {.value = imag_value, .variance = imag_variance}}; + + SUBCASE("Constructor") { + CHECK(decomposed.real_component.value == real_value); + CHECK(decomposed.real_component.variance == real_variance); + CHECK(decomposed.imag_component.value == imag_value); + CHECK(decomposed.imag_component.variance == imag_variance); + } + SUBCASE("Aggregate value") { + CHECK(real(decomposed.value()) == doctest::Approx(real_value)); + CHECK(imag(decomposed.value()) == doctest::Approx(imag_value)); + } + SUBCASE("Conversion to UniformComplexRandomVariable") { + auto const uniform = static_cast>(decomposed); + + CHECK(real(uniform.value) == doctest::Approx(real(decomposed.value()))); + CHECK(imag(uniform.value) == doctest::Approx(imag(decomposed.value()))); + CHECK(uniform.variance == doctest::Approx(real_variance + imag_variance)); + } + SUBCASE("Conversion to IndependentComplexRandomVariable") { + auto const independent = static_cast>(decomposed); + + CHECK(real(independent.value) == doctest::Approx(real(decomposed.value()))); + CHECK(imag(independent.value) == doctest::Approx(imag(decomposed.value()))); + CHECK(independent.variance == doctest::Approx(real_variance + imag_variance)); + } + } + } + SUBCASE("PolarComplexRandomVariable") { SUBCASE("Constructor") { for (auto const [magnitude, magnitude_variance, angle, angle_variance] : @@ -68,29 +316,7 @@ TEST_CASE("Test statistics") { } } - SUBCASE("Conversion to UniformComplexRandomVariable") { - for (auto const [magnitude, magnitude_variance, angle, angle_variance] : - std::array{std::tuple{1.0, 1.0, 0.0, 0.2}, std::tuple{2.0, 3.0, 0.0, 0.2}, - std::tuple{1.0, 1.0, pi / 2, 0.2}, std::tuple{1.0, 1.0, pi / 4, 0.2}}) { - CAPTURE(magnitude); - CAPTURE(magnitude_variance); - CAPTURE(angle); - CAPTURE(angle_variance); - - PolarComplexRandomVariable const polar{ - .magnitude = {.value = magnitude, .variance = magnitude_variance}, - .angle = {.value = angle, .variance = angle_variance}}; - - auto const uniform = static_cast>(polar); - - CHECK(real(uniform.value) == doctest::Approx(real(polar.value()))); - CHECK(imag(uniform.value) == doctest::Approx(imag(polar.value()))); - CHECK(uniform.variance == - doctest::Approx(polar.magnitude.variance + magnitude * magnitude * polar.angle.variance)); - } - } - - SUBCASE("Conversion to DecomposedIndependentComplexRandomVariable") { + SUBCASE("Conversion to DecomposedComplexRandomVariable") { for (auto const [magnitude, magnitude_variance, angle_variance] : std::array{std::tuple{1.0, 1.0, 0.2}, std::tuple{2.0, 1.0, 0.2}, std::tuple{1.0, 3.0, 0.2}, std::tuple{1.0, 2.0, 0.4}}) { @@ -103,16 +329,15 @@ TEST_CASE("Test statistics") { .magnitude = {.value = magnitude, .variance = magnitude}, .angle = {.value = 0.0, .variance = angle_variance}}; - auto const independent = - static_cast>(polar); + auto const decomposed = static_cast>(polar); - CHECK(independent.real_component.value == doctest::Approx(polar.magnitude.value)); - CHECK(independent.imag_component.value == doctest::Approx(0.0)); - CHECK(independent.real_component.variance == doctest::Approx(polar.magnitude.variance)); - CHECK(independent.imag_component.variance == + CHECK(decomposed.real_component.value == doctest::Approx(polar.magnitude.value)); + CHECK(decomposed.imag_component.value == doctest::Approx(0.0)); + CHECK(decomposed.real_component.variance == doctest::Approx(polar.magnitude.variance)); + CHECK(decomposed.imag_component.variance == doctest::Approx(magnitude * magnitude * polar.angle.variance)); - CHECK(real(independent.value()) == doctest::Approx(real(polar.value()))); - CHECK(imag(independent.value()) == doctest::Approx(imag(polar.value()))); + CHECK(real(decomposed.value()) == doctest::Approx(real(polar.value()))); + CHECK(imag(decomposed.value()) == doctest::Approx(imag(polar.value()))); } SUBCASE("Perpendicular phase shift") { @@ -120,16 +345,15 @@ TEST_CASE("Test statistics") { .magnitude = {.value = magnitude, .variance = magnitude}, .angle = {.value = pi / 2, .variance = angle_variance}}; - auto const independent = - static_cast>(polar); + auto const decomposed = static_cast>(polar); - CHECK(independent.real_component.value == doctest::Approx(0.0)); - CHECK(independent.imag_component.value == doctest::Approx(polar.magnitude.value)); - CHECK(independent.real_component.variance == + CHECK(decomposed.real_component.value == doctest::Approx(0.0)); + CHECK(decomposed.imag_component.value == doctest::Approx(polar.magnitude.value)); + CHECK(decomposed.real_component.variance == doctest::Approx(magnitude * magnitude * polar.angle.variance)); - CHECK(independent.imag_component.variance == doctest::Approx(polar.magnitude.variance)); - CHECK(real(independent.value()) == doctest::Approx(real(polar.value()))); - CHECK(imag(independent.value()) == doctest::Approx(imag(polar.value()))); + CHECK(decomposed.imag_component.variance == doctest::Approx(polar.magnitude.variance)); + CHECK(real(decomposed.value()) == doctest::Approx(real(polar.value()))); + CHECK(imag(decomposed.value()) == doctest::Approx(imag(polar.value()))); } SUBCASE("45deg phase shift") { @@ -137,19 +361,62 @@ TEST_CASE("Test statistics") { .magnitude = {.value = magnitude, .variance = magnitude}, .angle = {.value = pi / 4, .variance = angle_variance}}; - auto const independent = - static_cast>(polar); + auto const decomposed = static_cast>(polar); auto const uniform = static_cast>(polar); - CHECK(independent.real_component.value == doctest::Approx(real(uniform.value))); - CHECK(independent.imag_component.value == doctest::Approx(imag(uniform.value))); - CHECK(independent.real_component.variance == doctest::Approx(uniform.variance / 2)); - CHECK(independent.imag_component.variance == doctest::Approx(independent.real_component.variance)); - CHECK(real(independent.value()) == doctest::Approx(real(polar.value()))); - CHECK(imag(independent.value()) == doctest::Approx(imag(polar.value()))); + CHECK(decomposed.real_component.value == doctest::Approx(real(uniform.value))); + CHECK(decomposed.imag_component.value == doctest::Approx(imag(uniform.value))); + CHECK(decomposed.real_component.variance == doctest::Approx(uniform.variance / 2)); + CHECK(decomposed.imag_component.variance == doctest::Approx(decomposed.real_component.variance)); + CHECK(real(decomposed.value()) == doctest::Approx(real(polar.value()))); + CHECK(imag(decomposed.value()) == doctest::Approx(imag(polar.value()))); } } } + + SUBCASE("Conversion to IndependentComplexRandomVariable") { + for (auto const [magnitude, magnitude_variance, angle, angle_variance] : + std::array{std::tuple{1.0, 1.0, 0.0, 0.2}, std::tuple{2.0, 3.0, 0.0, 0.2}, + std::tuple{1.0, 1.0, pi / 2, 0.2}, std::tuple{1.0, 1.0, pi / 4, 0.2}}) { + CAPTURE(magnitude); + CAPTURE(magnitude_variance); + CAPTURE(angle); + CAPTURE(angle_variance); + + PolarComplexRandomVariable const polar{ + .magnitude = {.value = magnitude, .variance = magnitude_variance}, + .angle = {.value = angle, .variance = angle_variance}}; + + auto const independent = static_cast>(polar); + + CHECK(real(independent.value) == doctest::Approx(real(polar.value()))); + CHECK(imag(independent.value) == doctest::Approx(imag(polar.value()))); + CHECK(independent.variance == + doctest::Approx(polar.magnitude.variance + magnitude * magnitude * polar.angle.variance)); + } + } + + SUBCASE("Conversion to UniformComplexRandomVariable") { + for (auto const [magnitude, magnitude_variance, angle, angle_variance] : + std::array{std::tuple{1.0, 1.0, 0.0, 0.2}, std::tuple{2.0, 3.0, 0.0, 0.2}, + std::tuple{1.0, 1.0, pi / 2, 0.2}, std::tuple{1.0, 1.0, pi / 4, 0.2}}) { + CAPTURE(magnitude); + CAPTURE(magnitude_variance); + CAPTURE(angle); + CAPTURE(angle_variance); + + PolarComplexRandomVariable const polar{ + .magnitude = {.value = magnitude, .variance = magnitude_variance}, + .angle = {.value = angle, .variance = angle_variance}}; + + auto const uniform = static_cast>(polar); + + CHECK(real(uniform.value) == doctest::Approx(real(polar.value()))); + CHECK(imag(uniform.value) == doctest::Approx(imag(polar.value()))); + CHECK(uniform.variance == + doctest::Approx(polar.magnitude.variance + magnitude * magnitude * polar.angle.variance)); + } + } } } } // namespace power_grid_model From 909ebcf91699c006fe7a282060a6236f84e460b3 Mon Sep 17 00:00:00 2001 From: Nitish Bharambe Date: Tue, 18 Feb 2025 10:22:47 +0100 Subject: [PATCH 06/12] add implementation of current sensor Signed-off-by: Nitish Bharambe --- .../component/current_sensor.hpp | 38 ++++++++++--------- tests/cpp_unit_tests/test_current_sensor.cpp | 4 ++ 2 files changed, 24 insertions(+), 18 deletions(-) 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 eabc1948e9..0d7674af0e 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 @@ -10,6 +10,7 @@ #include "../common/common.hpp" #include "../common/enum.hpp" #include "../common/exception.hpp" +#include "../common/statistics.hpp" #include @@ -132,24 +133,22 @@ template class CurrentSensor : public Ge // TODO(mgovers) when filling the functions below take in mind that i_angle_sigma is optional CurrentSensorCalcParam sym_calc_param() const final { - if constexpr (is_asymmetric_v) { - return {}; // TODO - } else { - auto const i_polar = PolarComplexRandomVariable{ - .magnitude = {.value = i_measured_, .variance = i_sigma_ * i_sigma_}, - .angle = {.value = i_angle_measured_, .variance = i_angle_sigma_ * i_angle_sigma_}}; - auto const i_decomposed = static_cast>(i_polar); - - return CurrentSensorCalcParam{.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}; - } + auto const i_polar = PolarComplexRDV( + {i_measured_, i_sigma_ * i_sigma_}, {i_angle_measured_, i_angle_sigma_ * i_angle_sigma_}); + auto const i_decomposed = static_cast>(i_polar); + return CurrentSensorCalcParam{.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}; } CurrentSensorCalcParam asym_calc_param() const final { - CurrentSensorCalcParam calc_param{}; - // TODO - return calc_param; + auto const i_polar = PolarComplexRDV( + {i_measured_, i_sigma_ * i_sigma_}, {i_angle_measured_, i_angle_sigma_ * i_angle_sigma_}); + auto const i_decomposed = static_cast>(i_polar); + return CurrentSensorCalcParam{.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 get_sym_output(ComplexValue const& i) const final { return get_generic_output(i); @@ -160,8 +159,11 @@ template class CurrentSensor : public Ge template CurrentSensorOutput get_generic_output(ComplexValue const& i) const { CurrentSensorOutput output{}; - // TODO - (void)i; // Suppress unused variable warning + output.id = id(); + ComplexValue const i_residual{process_mean_val(i_measured_ - i) * base_current_}; + output.energized = 1; // power sensor is always energized + output.i_residual = cabs(i_residual); + output.i_angle_residual = arg(i_residual); return output; } }; diff --git a/tests/cpp_unit_tests/test_current_sensor.cpp b/tests/cpp_unit_tests/test_current_sensor.cpp index 7aa168c940..7b56de5cdb 100644 --- a/tests/cpp_unit_tests/test_current_sensor.cpp +++ b/tests/cpp_unit_tests/test_current_sensor.cpp @@ -6,6 +6,8 @@ #include +TEST_SUITE_BEGIN("test_current_sensor"); + namespace power_grid_model { namespace { auto const r_nan = RealValue{nan}; @@ -280,3 +282,5 @@ TEST_CASE("Test current sensor") { } } // namespace power_grid_model + +TEST_SUITE_END(); From 98839661da98328e75fe6704edf3c65219412175 Mon Sep 17 00:00:00 2001 From: Nitish Bharambe Date: Tue, 18 Feb 2025 16:52:01 +0100 Subject: [PATCH 07/12] add calc param Signed-off-by: Nitish Bharambe --- .../component/current_sensor.hpp | 28 ++++++++++++++----- 1 file changed, 21 insertions(+), 7 deletions(-) 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 0d7674af0e..97a51bdcf3 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 @@ -135,7 +135,14 @@ template class CurrentSensor : public Ge CurrentSensorCalcParam sym_calc_param() const final { auto const i_polar = PolarComplexRDV( {i_measured_, i_sigma_ * i_sigma_}, {i_angle_measured_, i_angle_sigma_ * i_angle_sigma_}); - auto const i_decomposed = static_cast>(i_polar); + auto const i_decomposed = [&i_polar]() { + if constexpr (is_symmetric_v) { + return static_cast>(i_polar); + } else { + return static_cast>( + static_cast>(i_polar)); + } + }(); return CurrentSensorCalcParam{.angle_measurement_type = angle_measurement_type(), .value = i_decomposed.value(), .i_real_variance = i_decomposed.real_component.variance, @@ -144,7 +151,14 @@ template class CurrentSensor : public Ge CurrentSensorCalcParam asym_calc_param() const final { auto const i_polar = PolarComplexRDV( {i_measured_, i_sigma_ * i_sigma_}, {i_angle_measured_, i_angle_sigma_ * i_angle_sigma_}); - auto const i_decomposed = static_cast>(i_polar); + auto const i_decomposed = [&i_polar]() { + if constexpr (is_symmetric_v) { + return static_cast>( + static_cast>(i_polar)); + } else { + return static_cast>(i_polar); + } + }(); return CurrentSensorCalcParam{.angle_measurement_type = angle_measurement_type(), .value = i_decomposed.value(), .i_real_variance = i_decomposed.real_component.variance, @@ -159,11 +173,11 @@ template class CurrentSensor : public Ge template CurrentSensorOutput get_generic_output(ComplexValue const& i) const { CurrentSensorOutput output{}; - output.id = id(); - ComplexValue const i_residual{process_mean_val(i_measured_ - i) * base_current_}; - output.energized = 1; // power sensor is always energized - output.i_residual = cabs(i_residual); - output.i_angle_residual = arg(i_residual); + // output.id = id(); + // ComplexValue const i_residual{process_mean_val(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; } }; From ea4fcd8cfacb4deb76bf8b2265a72c4787a63d9e Mon Sep 17 00:00:00 2001 From: Nitish Bharambe Date: Thu, 20 Feb 2025 10:00:59 +0100 Subject: [PATCH 08/12] wip current sensor Signed-off-by: Nitish Bharambe --- .../component/current_sensor.hpp | 12 +++--- tests/cpp_unit_tests/test_current_sensor.cpp | 39 ++++++++++--------- 2 files changed, 26 insertions(+), 25 deletions(-) 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 97a51bdcf3..756a9c73de 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 @@ -12,8 +12,6 @@ #include "../common/exception.hpp" #include "../common/statistics.hpp" -#include - namespace power_grid_model { class GenericCurrentSensor : public Sensor { @@ -173,11 +171,11 @@ template class CurrentSensor : public Ge template CurrentSensorOutput get_generic_output(ComplexValue const& i) const { CurrentSensorOutput output{}; - // output.id = id(); - // ComplexValue const i_residual{process_mean_val(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); + output.id = id(); + ComplexValue const i_residual{process_mean_val(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; } }; diff --git a/tests/cpp_unit_tests/test_current_sensor.cpp b/tests/cpp_unit_tests/test_current_sensor.cpp index 7b56de5cdb..3a7fa62cf0 100644 --- a/tests/cpp_unit_tests/test_current_sensor.cpp +++ b/tests/cpp_unit_tests/test_current_sensor.cpp @@ -32,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 sym_current_sensor_input{}; sym_current_sensor_input.id = 0; @@ -46,6 +46,8 @@ TEST_CASE("Test current sensor") { 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; ComplexValue const i_sym = (1.0 * 1e3 + 1i * 0.0) / base_current; ComplexValue const i_asym = i_sym * RealValue{1.0}; @@ -62,26 +64,27 @@ 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(pow(1.0 / base_current, 2))); - CHECK(sym_sensor_param.i_imag_variance == doctest::Approx(pow(0.2 * 1.0e3 / base_current, 2))); - CHECK(real(sym_sensor_param.value) == doctest::Approx(1.0e3 / base_current)); - CHECK(imag(sym_sensor_param.value) == doctest::Approx(0.0 / base_current)); + CHECK(sym_sensor_param.i_real_variance == doctest::Approx(pow(i_sigma_pu, 2))); + CHECK(sym_sensor_param.i_imag_variance == doctest::Approx(pow(0.2 * i_pu, 2))); + 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[0] == doctest::Approx(0.0)); - CHECK(asym_sensor_param.i_imag_variance[1] == 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_sigma_pu * 3.0)); + // CHECK(asym_sensor_param.i_imag_variance[1] == doctest::Approx((i_sigma_pu * sin(deg_240) * sin(deg_240) + + // i_pu * i_pu * 0.2 * cos(deg_240) * cos(deg_240)) /3.0)); + CHECK(real(asym_sensor_param.value[0]) == doctest::Approx(i_pu / 3.0)); + CHECK(imag(asym_sensor_param.value[1]) == doctest::Approx(i_pu * sin(deg_240) / 3.0)); + + 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); From decc9e4128a7a92ec389b7d26ec7b894891aad70 Mon Sep 17 00:00:00 2001 From: Nitish Bharambe Date: Mon, 24 Feb 2025 10:43:31 +0100 Subject: [PATCH 09/12] finish all tests Signed-off-by: Nitish Bharambe --- .../component/current_sensor.hpp | 38 +++++-------------- tests/cpp_unit_tests/test_current_sensor.cpp | 17 +++++---- 2 files changed, 19 insertions(+), 36 deletions(-) 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 756a9c73de..d39fdaa85c 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 @@ -130,37 +130,17 @@ template class CurrentSensor : public Ge // TODO(mgovers) when filling the functions below take in mind that i_angle_sigma is optional - CurrentSensorCalcParam sym_calc_param() const final { - auto const i_polar = PolarComplexRDV( - {i_measured_, i_sigma_ * i_sigma_}, {i_angle_measured_, i_angle_sigma_ * i_angle_sigma_}); - auto const i_decomposed = [&i_polar]() { - if constexpr (is_symmetric_v) { - return static_cast>(i_polar); - } else { - return static_cast>( - static_cast>(i_polar)); - } - }(); - return CurrentSensorCalcParam{.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}; - } - CurrentSensorCalcParam asym_calc_param() const final { + CurrentSensorCalcParam sym_calc_param() const final { return calc_decomposed_param(); } + CurrentSensorCalcParam asym_calc_param() const final { return calc_decomposed_param(); } + + template CurrentSensorCalcParam calc_decomposed_param() const { auto const i_polar = PolarComplexRDV( {i_measured_, i_sigma_ * i_sigma_}, {i_angle_measured_, i_angle_sigma_ * i_angle_sigma_}); - auto const i_decomposed = [&i_polar]() { - if constexpr (is_symmetric_v) { - return static_cast>( - static_cast>(i_polar)); - } else { - return static_cast>(i_polar); - } - }(); - return CurrentSensorCalcParam{.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}; + auto const i_decomposed = DecomposedComplexRDV(i_polar); + return CurrentSensorCalcParam{.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 get_sym_output(ComplexValue const& i) const final { return get_generic_output(i); diff --git a/tests/cpp_unit_tests/test_current_sensor.cpp b/tests/cpp_unit_tests/test_current_sensor.cpp index 3a7fa62cf0..f4bc58fd80 100644 --- a/tests/cpp_unit_tests/test_current_sensor.cpp +++ b/tests/cpp_unit_tests/test_current_sensor.cpp @@ -48,6 +48,8 @@ TEST_CASE("Test current sensor") { 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 const i_sym = (1.0 * 1e3 + 1i * 0.0) / base_current; ComplexValue const i_asym = i_sym * RealValue{1.0}; @@ -64,8 +66,8 @@ 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(pow(i_sigma_pu, 2))); - CHECK(sym_sensor_param.i_imag_variance == doctest::Approx(pow(0.2 * i_pu, 2))); + 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)); @@ -75,11 +77,12 @@ TEST_CASE("Test current sensor") { 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[0] == doctest::Approx(i_sigma_pu * 3.0)); - // CHECK(asym_sensor_param.i_imag_variance[1] == doctest::Approx((i_sigma_pu * sin(deg_240) * sin(deg_240) + - // i_pu * i_pu * 0.2 * cos(deg_240) * cos(deg_240)) /3.0)); - CHECK(real(asym_sensor_param.value[0]) == doctest::Approx(i_pu / 3.0)); - CHECK(imag(asym_sensor_param.value[1]) == doctest::Approx(i_pu * sin(deg_240) / 3.0)); + 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); From 4d37735354853985335f7bd324bed0ce2d34ec25 Mon Sep 17 00:00:00 2001 From: Nitish Bharambe <78108900+nitbharambe@users.noreply.github.com> Date: Tue, 25 Feb 2025 11:17:58 +0100 Subject: [PATCH 10/12] Update power_grid_model_c/power_grid_model/include/power_grid_model/component/current_sensor.hpp Signed-off-by: Nitish Bharambe <78108900+nitbharambe@users.noreply.github.com> --- .../include/power_grid_model/component/current_sensor.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 d39fdaa85c..57caec78dc 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 @@ -134,9 +134,9 @@ template class CurrentSensor : public Ge CurrentSensorCalcParam asym_calc_param() const final { return calc_decomposed_param(); } template CurrentSensorCalcParam calc_decomposed_param() const { - auto const i_polar = PolarComplexRDV( + auto const i_polar = PolarComplexRandVar( {i_measured_, i_sigma_ * i_sigma_}, {i_angle_measured_, i_angle_sigma_ * i_angle_sigma_}); - auto const i_decomposed = DecomposedComplexRDV(i_polar); + auto const i_decomposed = DecomposedComplexRandVar(i_polar); return CurrentSensorCalcParam{.angle_measurement_type = angle_measurement_type(), .value = i_decomposed.value(), .i_real_variance = i_decomposed.real_component.variance, From 273bbcd699974c7205542e43efc6444f93f2b1b5 Mon Sep 17 00:00:00 2001 From: Nitish Bharambe Date: Thu, 27 Feb 2025 09:54:11 +0100 Subject: [PATCH 11/12] change order of initialization Signed-off-by: Nitish Bharambe --- .../include/power_grid_model/component/current_sensor.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 57caec78dc..c9e46b7a10 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 @@ -123,10 +123,10 @@ template class CurrentSensor : public Ge private: double base_current_{}; double base_current_inv_{}; - RealValue i_measured_{}; RealValue i_angle_measured_{}; - double i_sigma_{}; double i_angle_sigma_{}; + double i_sigma_{}; + RealValue i_measured_{}; // TODO(mgovers) when filling the functions below take in mind that i_angle_sigma is optional From 83838ec29be8649af3f5c338db76755ca80d0c93 Mon Sep 17 00:00:00 2001 From: Nitish Bharambe Date: Thu, 27 Feb 2025 15:12:30 +0100 Subject: [PATCH 12/12] remove todos Signed-off-by: Nitish Bharambe --- .../include/power_grid_model/component/current_sensor.hpp | 2 -- tests/cpp_unit_tests/test_current_sensor.cpp | 2 +- 2 files changed, 1 insertion(+), 3 deletions(-) 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 c9e46b7a10..a6c125bc48 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 @@ -128,8 +128,6 @@ template class CurrentSensor : public Ge double i_sigma_{}; RealValue i_measured_{}; - // TODO(mgovers) when filling the functions below take in mind that i_angle_sigma is optional - CurrentSensorCalcParam sym_calc_param() const final { return calc_decomposed_param(); } CurrentSensorCalcParam asym_calc_param() const final { return calc_decomposed_param(); } diff --git a/tests/cpp_unit_tests/test_current_sensor.cpp b/tests/cpp_unit_tests/test_current_sensor.cpp index f4bc58fd80..4e9e5df487 100644 --- a/tests/cpp_unit_tests/test_current_sensor.cpp +++ b/tests/cpp_unit_tests/test_current_sensor.cpp @@ -32,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); // TODO check why this inerferes with normal output + CAPTURE(terminal_type); CurrentSensorInput sym_current_sensor_input{}; sym_current_sensor_input.id = 0;