From 32d404e98f8eccf6f5b72ec5c4a0ecebf539714d Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Tue, 1 Apr 2025 15:44:49 +0200 Subject: [PATCH 01/25] add basic current measurement support to measured values Signed-off-by: Martijn Govers --- .../calculation_parameters.hpp | 2 + .../math_solver/measured_values.hpp | 83 ++++++++++++++----- tests/cpp_unit_tests/test_observability.cpp | 45 +++++++--- 3 files changed, 94 insertions(+), 36 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 772b9edf93..5929941044 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 @@ -235,6 +235,8 @@ template struct StateEstimationInput { std::vector> measured_branch_from_power; std::vector> measured_branch_to_power; std::vector> measured_bus_injection; + std::vector> measured_branch_from_current; + std::vector> measured_branch_to_current; }; struct ShortCircuitInput { diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp index 2c69821f06..49dfb17629 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp @@ -72,8 +72,12 @@ template class MeasuredValues { constexpr bool has_voltage(Idx bus) const { return idx_voltage_[bus] >= 0; } constexpr bool has_angle_measurement(Idx bus) const { return !is_nan(imag(voltage(bus))); } constexpr bool has_bus_injection(Idx bus) const { return bus_injection_[bus].idx_bus_injection >= 0; } - constexpr bool has_branch_from(Idx branch) const { return idx_branch_from_power_[branch] >= 0; } - constexpr bool has_branch_to(Idx branch) const { return idx_branch_to_power_[branch] >= 0; } + constexpr bool has_branch_from(Idx branch) const { + return idx_branch_from_power_[branch] >= 0 || idx_branch_from_current_[branch]; + } + constexpr bool has_branch_to(Idx branch) const { + return idx_branch_to_power_[branch] >= 0 || idx_branch_to_current_[branch] >= 0; + } constexpr bool has_shunt(Idx shunt) const { return idx_shunt_power_[shunt] >= 0; } constexpr bool has_load_gen(Idx load_gen) const { return idx_load_gen_power_[load_gen] >= 0; } constexpr bool has_source(Idx source) const { return idx_source_power_[source] >= 0; } @@ -91,6 +95,12 @@ template class MeasuredValues { return power_main_value_[idx_branch_from_power_[branch]]; } constexpr auto const& branch_to_power(Idx branch) const { return power_main_value_[idx_branch_to_power_[branch]]; } + constexpr auto const& branch_from_current(Idx branch) const { + return current_main_value_[idx_branch_from_current_[branch]]; + } + constexpr auto const& branch_to_current(Idx branch) const { + return current_main_value_[idx_branch_to_current_[branch]]; + } constexpr auto const& shunt_power(Idx shunt) const { return power_main_value_[idx_shunt_power_[shunt]]; } constexpr auto const& load_gen_power(Idx load_gen) const { return extra_value_[idx_load_gen_power_[load_gen]]; } constexpr auto const& source_power(Idx source) const { return extra_value_[idx_source_power_[source]]; } @@ -177,6 +187,7 @@ template class MeasuredValues { // branch/shunt flow, bus voltage, injection flow std::vector> voltage_main_value_; std::vector> power_main_value_; + std::vector> current_main_value_; // flat array of all the load_gen/source measurement // not relevant for the main calculation, as extra data for load_gen/source calculation @@ -197,6 +208,9 @@ template class MeasuredValues { // relevant for extra value IdxVector idx_load_gen_power_; IdxVector idx_source_power_; + // current measurement + IdxVector idx_branch_from_current_; + IdxVector idx_branch_to_current_; Idx n_voltage_measurements_{}; Idx n_voltage_angle_measurements_{}; @@ -211,7 +225,7 @@ template class MeasuredValues { void process_bus_related_measurements(StateEstimationInput const& input) { /* - The main purpose of this function is to aggregate all voltage and power sensor values to + The main purpose of this function is to aggregate all voltage and power/current sensor values to one voltage sensor value per bus. one injection power sensor value per bus. one power sensor value per shunt (in injection reference direction, note shunt itself is not considered as @@ -383,36 +397,44 @@ template class MeasuredValues { void process_branch_measurements(StateEstimationInput const& input) { /* - The main purpose of this function is to aggregate all power sensor values to one power sensor value per branch - side. + The main purpose of this function is to aggregate all power/current sensor values to one power/current sensor + value per branch side. This function loops through all branches. The branch_bus_idx contains the from and to bus indexes of the branch, or disconnected if the branch is not connected at that side. For each branch the checker checks if the from and to side are connected by checking if branch_bus_idx = disconnected. - If the branch_bus_idx = disconnected, idx_branch_to_power_/idx_branch_from_power_ is set to disconnected. + If the branch_bus_idx = disconnected, idx_branch_(to|from)_(power|current)_ is set to disconnected. If the side is connected, but there are no measurements in this branch side - idx_branch_to_power_/idx_branch_from_power_ is set to disconnected. - Else, idx_branch_to_power_/idx_branch_from_power_ is set to the index of the aggregated data in - power_main_value_. + idx_branch_(to|from)_(power|current)_ is set to disconnected. + Else, idx_branch_(to|from)_(power|current)_ is set to the index of the aggregated data in + power/current_main_value_. All measurement values for a single side of a branch are combined in a weighted average, which is appended to - power_main_value_. The values in power_main_value_ can be found using - idx_branch_to_power_/idx_branch_from_power_. + power/current_main_value_. The values in power/current_main_value_ can be found using + idx_branch_(to|from)_(power|current)_. */ MathModelTopology const& topo = math_topology(); static constexpr auto branch_from_checker = [](BranchIdx x) { return x[0] != -1; }; static constexpr auto branch_to_checker = [](BranchIdx x) { return x[1] != -1; }; for (Idx const branch : IdxRange{topo.n_branch()}) { - // from side + // from side power idx_branch_from_power_[branch] = process_one_object(branch, topo.power_sensors_per_branch_from, topo.branch_bus_idx, input.measured_branch_from_power, power_main_value_, branch_from_checker); - // to side + // to side current idx_branch_to_power_[branch] = process_one_object(branch, topo.power_sensors_per_branch_to, topo.branch_bus_idx, input.measured_branch_to_power, power_main_value_, branch_to_checker); + // from side power + idx_branch_from_current_[branch] = + process_one_object(branch, topo.current_sensors_per_branch_from, topo.branch_bus_idx, + input.measured_branch_from_current, current_main_value_, branch_from_checker); + // to side current + idx_branch_to_current_[branch] = + process_one_object(branch, topo.current_sensors_per_branch_to, topo.branch_bus_idx, + input.measured_branch_to_current, current_main_value_, branch_to_checker); } } @@ -456,8 +478,9 @@ template class MeasuredValues { // Since p and q are entirely decoupled, the real and imaginary components accumulate separately template requires(!only_magnitude) - static PowerSensorCalcParam combine_measurements(std::vector> const& data, - IdxRange const& sensors) { + static DecomposedComplexRandVar combine_measurements(std::ranges::range auto const& measurements) + requires std::same_as, DecomposedComplexRandVar> + { struct AccumulatedValues { RealValue inverse_p_variance{}; RealValue inverse_q_variance{}; @@ -466,8 +489,7 @@ template class MeasuredValues { }; AccumulatedValues accumulated; - std::ranges::for_each(sensors, [&data, &accumulated](auto pos) { - DecomposedComplexRandVar const& measurement = data[pos]; + std::ranges::for_each(measurements, [&accumulated](auto const& measurement) { accumulated.inverse_p_variance += RealValue{1.0} / measurement.real_component.variance; accumulated.inverse_q_variance += RealValue{1.0} / measurement.imag_component.variance; @@ -476,18 +498,34 @@ template class MeasuredValues { }); if (is_normal(accumulated.inverse_p_variance) && is_normal(accumulated.inverse_q_variance)) { - return PowerSensorCalcParam{ + return DecomposedComplexRandVar{ .real_component = {.value = accumulated.p_value / accumulated.inverse_p_variance, .variance = RealValue{1.0} / accumulated.inverse_p_variance}, .imag_component = {.value = accumulated.q_value / accumulated.inverse_q_variance, .variance = RealValue{1.0} / accumulated.inverse_q_variance}}; } - return PowerSensorCalcParam{ + return DecomposedComplexRandVar{ .real_component = {.value = accumulated.p_value, .variance = RealValue{std::numeric_limits::infinity()}}, .imag_component = {.value = accumulated.q_value, .variance = RealValue{std::numeric_limits::infinity()}}}; } + template + requires(!only_magnitude) + static PowerSensorCalcParam combine_measurements(std::vector> const& data, + IdxRange const& sensors) { + return combine_measurements( + std::views::transform(sensors, [&](Idx pos) -> auto& { return data[pos]; })); + } + template + requires(!only_magnitude) + static CurrentSensorCalcParam combine_measurements(std::vector> const& data, + IdxRange const& sensors) { + return {.angle_measurement_type = + data[sensors.front()].angle_measurement_type, // TODO(mgovers): this should be fixed + .measurement = combine_measurements( + std::views::transform(sensors, [&](Idx pos) -> auto& { return data[pos].measurement; }))}; + } template static auto combine_measurements(std::vector const& data) { @@ -512,11 +550,10 @@ template class MeasuredValues { static constexpr DefaultStatusChecker default_status_checker{}; - template + template static Idx process_one_object(Idx const object, grouped_idx_vector_type auto const& sensors_per_object, - std::vector const& object_status, - std::vector> const& input_data, - std::vector>& result_data, + std::vector const& object_status, std::vector const& input_data, + std::vector& result_data, StatusChecker status_checker = default_status_checker) { if (!status_checker(object_status[object])) { return disconnected; diff --git a/tests/cpp_unit_tests/test_observability.cpp b/tests/cpp_unit_tests/test_observability.cpp index fec8450788..628f3907ab 100644 --- a/tests/cpp_unit_tests/test_observability.cpp +++ b/tests/cpp_unit_tests/test_observability.cpp @@ -11,16 +11,32 @@ namespace power_grid_model { namespace { -void check_not_observable(MathModelTopology const& topo, MathModelParam const& param, - StateEstimationInput const& se_input) { +void check_whether_observable(bool is_observable, MathModelTopology const& topo, + MathModelParam const& param, + StateEstimationInput const& se_input) { auto topo_ptr = std::make_shared(topo); auto param_ptr = std::make_shared const>(param); YBus const y_bus{topo_ptr, param_ptr}; math_solver::MeasuredValues const measured_values{y_bus.shared_topology(), se_input}; - CHECK_THROWS_AS( - math_solver::necessary_observability_check(measured_values, y_bus.math_topology(), y_bus.y_bus_structure()), - NotObservableError); + if (is_observable) { + CHECK_NOTHROW(math_solver::necessary_observability_check(measured_values, y_bus.math_topology(), + y_bus.y_bus_structure())); + } else { + CHECK_THROWS_AS( + math_solver::necessary_observability_check(measured_values, y_bus.math_topology(), y_bus.y_bus_structure()), + NotObservableError); + } +} + +void check_observable(MathModelTopology const& topo, MathModelParam const& param, + StateEstimationInput const& se_input) { + check_whether_observable(true, topo, param, se_input); +} + +void check_not_observable(MathModelTopology const& topo, MathModelParam const& param, + StateEstimationInput const& se_input) { + check_whether_observable(false, topo, param, se_input); } } // namespace @@ -45,12 +61,13 @@ TEST_CASE("Necessary observability check") { topo.power_sensors_per_shunt = {from_sparse, {0}}; topo.power_sensors_per_branch_from = {from_sparse, {0, 1, 1, 1}}; topo.power_sensors_per_branch_to = {from_sparse, {0, 0, 0, 0}}; + topo.current_sensors_per_branch_from = {from_sparse, {0, 0, 0, 0}}; + topo.current_sensors_per_branch_to = {from_sparse, {0, 0, 0, 0}}; topo.voltage_sensors_per_bus = {from_sparse, {0, 1, 1, 1}}; MathModelParam param; param.source_param = {SourceCalcParam{.y1 = 10.0 - 50.0i, .y0 = 10.0 - 50.0i}}; param.branch_param = {{1.0, -1.0, -1.0, 1.0}, {1.0, -1.0, -1.0, 1.0}, {1.0, -1.0, -1.0, 1.0}}; - auto param_ptr = std::make_shared const>(param); StateEstimationInput se_input; se_input.source_status = {1}; @@ -61,13 +78,7 @@ TEST_CASE("Necessary observability check") { se_input.measured_branch_from_power = { {.real_component = {.value = 1.0, .variance = 1.0}, .imag_component = {.value = 0.0, .variance = 1.0}}}; - SUBCASE("Observable grid") { - auto topo_ptr = std::make_shared(topo); - YBus const y_bus{topo_ptr, param_ptr}; - math_solver::MeasuredValues const measured_values{y_bus.shared_topology(), se_input}; - CHECK_NOTHROW(math_solver::necessary_observability_check(measured_values, y_bus.math_topology(), - y_bus.y_bus_structure())); - } + SUBCASE("Observable grid") { check_observable(topo, param, se_input); } SUBCASE("No voltage sensor") { topo.voltage_sensors_per_bus = {from_sparse, {0, 0, 0, 0}}; @@ -107,6 +118,14 @@ TEST_CASE("Necessary observability check") { // this will throw NotObservableError check_not_observable(topo, param, se_input); } + SUBCASE("Power sensor and current sensor are equivalent") { + // remove the power sensor -> not observable + topo.power_sensors_per_branch_from = {from_sparse, {0, 0, 0, 0}}; + check_not_observable(topo, param, se_input); + // add the current sensor -> observable + topo.current_sensors_per_branch_from = {from_sparse, {0, 1, 1, 1}}; + check_observable(topo, param, se_input); + } } } // namespace power_grid_model From 968786d72af77a729336783562328be06fa32d8b Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Tue, 1 Apr 2025 16:10:40 +0200 Subject: [PATCH 02/25] separate accumulation for real and imag part Signed-off-by: Martijn Govers --- .../math_solver/measured_values.hpp | 63 ++++++++++--------- 1 file changed, 33 insertions(+), 30 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp index 49dfb17629..96fdc17e70 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp @@ -476,46 +476,49 @@ template class MeasuredValues { // set imag part to nan, to signal this is a magnitude only measurement // // Since p and q are entirely decoupled, the real and imaginary components accumulate separately + template + requires(!only_magnitude) + static IndependentRealRandVar combine_measurements(std::ranges::range auto const& measurements) + requires std::same_as, IndependentRealRandVar> + { + RealValue accumulated_inverse_variance{}; + RealValue accumulated_value{}; + + std::ranges::for_each(measurements, + [&accumulated_inverse_variance, &accumulated_value](auto const& measurement) { + accumulated_inverse_variance += RealValue{1.0} / measurement.variance; + accumulated_value += measurement.value / measurement.variance; + }); + + if (is_normal(accumulated_inverse_variance)) { + return IndependentRealRandVar{.value = accumulated_value / accumulated_inverse_variance, + .variance = RealValue{1.0} / accumulated_inverse_variance}; + } + return {.value = accumulated_value, .variance = RealValue{std::numeric_limits::infinity()}}; + } template requires(!only_magnitude) static DecomposedComplexRandVar combine_measurements(std::ranges::range auto const& measurements) requires std::same_as, DecomposedComplexRandVar> { - struct AccumulatedValues { - RealValue inverse_p_variance{}; - RealValue inverse_q_variance{}; - RealValue p_value{}; - RealValue q_value{}; - }; - - AccumulatedValues accumulated; - std::ranges::for_each(measurements, [&accumulated](auto const& measurement) { - accumulated.inverse_p_variance += RealValue{1.0} / measurement.real_component.variance; - accumulated.inverse_q_variance += RealValue{1.0} / measurement.imag_component.variance; - - accumulated.p_value += measurement.real_component.value / measurement.real_component.variance; - accumulated.q_value += measurement.imag_component.value / measurement.imag_component.variance; - }); - - if (is_normal(accumulated.inverse_p_variance) && is_normal(accumulated.inverse_q_variance)) { - return DecomposedComplexRandVar{ - .real_component = {.value = accumulated.p_value / accumulated.inverse_p_variance, - .variance = RealValue{1.0} / accumulated.inverse_p_variance}, - .imag_component = {.value = accumulated.q_value / accumulated.inverse_q_variance, - .variance = RealValue{1.0} / accumulated.inverse_q_variance}}; + auto result = DecomposedComplexRandVar{ + .real_component = combine_measurements( + std::views::transform(measurements, [](auto const& x) -> auto& { return x.real_component; })), + .imag_component = combine_measurements( + std::views::transform(measurements, [](auto const& x) -> auto& { return x.imag_component; }))}; + + if (!(is_normal(result.real_component.variance) && is_normal(result.imag_component.variance))) { + result.real_component.variance = RealValue{std::numeric_limits::infinity()}; + result.imag_component.variance = RealValue{std::numeric_limits::infinity()}; } - return DecomposedComplexRandVar{ - .real_component = {.value = accumulated.p_value, - .variance = RealValue{std::numeric_limits::infinity()}}, - .imag_component = {.value = accumulated.q_value, - .variance = RealValue{std::numeric_limits::infinity()}}}; + return result; } template requires(!only_magnitude) static PowerSensorCalcParam combine_measurements(std::vector> const& data, IdxRange const& sensors) { - return combine_measurements( - std::views::transform(sensors, [&](Idx pos) -> auto& { return data[pos]; })); + return combine_measurements(sensors | + std::views::transform([&](Idx pos) -> auto& { return data[pos]; })); } template requires(!only_magnitude) @@ -524,7 +527,7 @@ template class MeasuredValues { return {.angle_measurement_type = data[sensors.front()].angle_measurement_type, // TODO(mgovers): this should be fixed .measurement = combine_measurements( - std::views::transform(sensors, [&](Idx pos) -> auto& { return data[pos].measurement; }))}; + sensors | std::views::transform([&](Idx pos) -> auto& { return data[pos].measurement; }))}; } template From 19b6febad902254d2a305ad66dd2462b8fdc9716 Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Tue, 1 Apr 2025 16:13:42 +0200 Subject: [PATCH 03/25] minor Signed-off-by: Martijn Govers --- .../include/power_grid_model/math_solver/measured_values.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp index 96fdc17e70..ad6811895e 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp @@ -503,9 +503,9 @@ template class MeasuredValues { { auto result = DecomposedComplexRandVar{ .real_component = combine_measurements( - std::views::transform(measurements, [](auto const& x) -> auto& { return x.real_component; })), + measurements | std::views::transform([](auto const& x) -> auto& { return x.real_component; })), .imag_component = combine_measurements( - std::views::transform(measurements, [](auto const& x) -> auto& { return x.imag_component; }))}; + measurements | std::views::transform([](auto const& x) -> auto& { return x.imag_component; }))}; if (!(is_normal(result.real_component.variance) && is_normal(result.imag_component.variance))) { result.real_component.variance = RealValue{std::numeric_limits::infinity()}; From 086403d6421efe400a128ba7ec9e74b6e59b3949 Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Thu, 3 Apr 2025 15:57:38 +0200 Subject: [PATCH 04/25] minor Signed-off-by: Martijn Govers --- .../include/power_grid_model/math_solver/measured_values.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp index ad6811895e..6dacafd06d 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp @@ -478,7 +478,7 @@ template class MeasuredValues { // Since p and q are entirely decoupled, the real and imaginary components accumulate separately template requires(!only_magnitude) - static IndependentRealRandVar combine_measurements(std::ranges::range auto const& measurements) + static IndependentRealRandVar combine_measurements(std::ranges::view auto measurements) requires std::same_as, IndependentRealRandVar> { RealValue accumulated_inverse_variance{}; @@ -498,7 +498,7 @@ template class MeasuredValues { } template requires(!only_magnitude) - static DecomposedComplexRandVar combine_measurements(std::ranges::range auto const& measurements) + static DecomposedComplexRandVar combine_measurements(std::ranges::view auto measurements) requires std::same_as, DecomposedComplexRandVar> { auto result = DecomposedComplexRandVar{ From 4af2e065c8aa81dbb2e37dd6f72a100cd7342136 Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Thu, 3 Apr 2025 16:40:47 +0200 Subject: [PATCH 05/25] combine measurements in statistics.hpp Signed-off-by: Martijn Govers --- .../power_grid_model/common/statistics.hpp | 45 +++++++++++++++++ .../math_solver/measured_values.hpp | 50 ++----------------- 2 files changed, 48 insertions(+), 47 deletions(-) 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 a4fdab9695..4345475af4 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 @@ -7,6 +7,8 @@ #include "common.hpp" #include "three_phase_tensor.hpp" +#include + /** * @file statistics.hpp * @brief This file contains various structures and functions for handling statistical representations of @@ -204,4 +206,47 @@ template struct PolarComplexRandVar { 9.0}}; } }; + +namespace statistics { +// combine multiple measurements of one quantity using Kalman filter +constexpr auto combine_measurements(std::ranges::view auto measurements) + requires std::same_as, + IndependentRealRandVar::sym>> +{ + using sym = std::ranges::range_value_t::sym; + + RealValue accumulated_inverse_variance{}; + RealValue accumulated_value{}; + + std::ranges::for_each(measurements, [&accumulated_inverse_variance, &accumulated_value](auto const& measurement) { + accumulated_inverse_variance += RealValue{1.0} / measurement.variance; + accumulated_value += measurement.value / measurement.variance; + }); + + if (is_normal(accumulated_inverse_variance)) { + return IndependentRealRandVar{.value = accumulated_value / accumulated_inverse_variance, + .variance = RealValue{1.0} / accumulated_inverse_variance}; + } + return IndependentRealRandVar{.value = accumulated_value, + .variance = RealValue{std::numeric_limits::infinity()}}; +} +constexpr auto combine_measurements(std::ranges::view auto measurements) + requires std::same_as, + DecomposedComplexRandVar::sym>> +{ + using sym = std::ranges::range_value_t::sym; + + DecomposedComplexRandVar result{ + .real_component = combine_measurements( + measurements | std::views::transform([](auto const& x) -> auto& { return x.real_component; })), + .imag_component = combine_measurements( + measurements | std::views::transform([](auto const& x) -> auto& { return x.imag_component; }))}; + + if (!(is_normal(result.real_component.variance) && is_normal(result.imag_component.variance))) { + result.real_component.variance = RealValue{std::numeric_limits::infinity()}; + result.imag_component.variance = RealValue{std::numeric_limits::infinity()}; + } + return result; +} +} // namespace statistics } // namespace power_grid_model diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp index 6dacafd06d..b271b29847 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp @@ -469,56 +469,12 @@ template class MeasuredValues { return VoltageSensorCalcParam{accumulated_value / accumulated_inverse_variance, 1.0 / accumulated_inverse_variance}; } - - // combine multiple measurements of one quantity - // using Kalman filter - // if only_magnitude = true, combine the abs value of the individual data - // set imag part to nan, to signal this is a magnitude only measurement - // - // Since p and q are entirely decoupled, the real and imaginary components accumulate separately - template - requires(!only_magnitude) - static IndependentRealRandVar combine_measurements(std::ranges::view auto measurements) - requires std::same_as, IndependentRealRandVar> - { - RealValue accumulated_inverse_variance{}; - RealValue accumulated_value{}; - - std::ranges::for_each(measurements, - [&accumulated_inverse_variance, &accumulated_value](auto const& measurement) { - accumulated_inverse_variance += RealValue{1.0} / measurement.variance; - accumulated_value += measurement.value / measurement.variance; - }); - - if (is_normal(accumulated_inverse_variance)) { - return IndependentRealRandVar{.value = accumulated_value / accumulated_inverse_variance, - .variance = RealValue{1.0} / accumulated_inverse_variance}; - } - return {.value = accumulated_value, .variance = RealValue{std::numeric_limits::infinity()}}; - } - template - requires(!only_magnitude) - static DecomposedComplexRandVar combine_measurements(std::ranges::view auto measurements) - requires std::same_as, DecomposedComplexRandVar> - { - auto result = DecomposedComplexRandVar{ - .real_component = combine_measurements( - measurements | std::views::transform([](auto const& x) -> auto& { return x.real_component; })), - .imag_component = combine_measurements( - measurements | std::views::transform([](auto const& x) -> auto& { return x.imag_component; }))}; - - if (!(is_normal(result.real_component.variance) && is_normal(result.imag_component.variance))) { - result.real_component.variance = RealValue{std::numeric_limits::infinity()}; - result.imag_component.variance = RealValue{std::numeric_limits::infinity()}; - } - return result; - } template requires(!only_magnitude) static PowerSensorCalcParam combine_measurements(std::vector> const& data, IdxRange const& sensors) { - return combine_measurements(sensors | - std::views::transform([&](Idx pos) -> auto& { return data[pos]; })); + return statistics::combine_measurements(sensors | + std::views::transform([&](Idx pos) -> auto& { return data[pos]; })); } template requires(!only_magnitude) @@ -526,7 +482,7 @@ template class MeasuredValues { IdxRange const& sensors) { return {.angle_measurement_type = data[sensors.front()].angle_measurement_type, // TODO(mgovers): this should be fixed - .measurement = combine_measurements( + .measurement = statistics::combine_measurements( sensors | std::views::transform([&](Idx pos) -> auto& { return data[pos].measurement; }))}; } From 286dc2710dc8982952d46007758aca007c2a6205 Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Thu, 3 Apr 2025 17:27:04 +0200 Subject: [PATCH 06/25] measured values basic support for local current measurement Signed-off-by: Martijn Govers --- .../power_grid_model/common/statistics.hpp | 24 ++++++---- .../math_solver/measured_values.hpp | 46 +++++++++---------- tests/cpp_unit_tests/test_observability.cpp | 8 +++- 3 files changed, 43 insertions(+), 35 deletions(-) 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 4345475af4..784f81909c 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 @@ -211,24 +211,30 @@ namespace statistics { // combine multiple measurements of one quantity using Kalman filter constexpr auto combine_measurements(std::ranges::view auto measurements) requires std::same_as, - IndependentRealRandVar::sym>> + UniformRealRandVar::sym>> || + std::same_as, + IndependentRealRandVar::sym>> || + std::same_as, + UniformComplexRandVar::sym>> { - using sym = std::ranges::range_value_t::sym; + using RandVarType = std::ranges::range_value_t; + using ValueType = decltype(RandVarType::value); + using VarianceType = decltype(RandVarType::variance); + using sym = RandVarType::sym; - RealValue accumulated_inverse_variance{}; - RealValue accumulated_value{}; + VarianceType accumulated_inverse_variance{}; + ValueType accumulated_value{}; std::ranges::for_each(measurements, [&accumulated_inverse_variance, &accumulated_value](auto const& measurement) { - accumulated_inverse_variance += RealValue{1.0} / measurement.variance; + accumulated_inverse_variance += VarianceType{1.0} / measurement.variance; accumulated_value += measurement.value / measurement.variance; }); if (is_normal(accumulated_inverse_variance)) { - return IndependentRealRandVar{.value = accumulated_value / accumulated_inverse_variance, - .variance = RealValue{1.0} / accumulated_inverse_variance}; + return RandVarType{.value = accumulated_value / accumulated_inverse_variance, + .variance = VarianceType{1.0} / accumulated_inverse_variance}; } - return IndependentRealRandVar{.value = accumulated_value, - .variance = RealValue{std::numeric_limits::infinity()}}; + return RandVarType{.value = accumulated_value, .variance = VarianceType{std::numeric_limits::infinity()}}; } constexpr auto combine_measurements(std::ranges::view auto measurements) requires std::same_as, diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp index b271b29847..715d735d1b 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp @@ -54,6 +54,8 @@ template class MeasuredValues { idx_shunt_power_(math_topology().n_shunt()), idx_load_gen_power_(math_topology().n_load_gen()), idx_source_power_(math_topology().n_source()), + idx_branch_from_current_(math_topology().n_branch()), + idx_branch_to_current_(math_topology().n_branch()), // default angle shift // sym: 0 // asym: 0, -120deg, -240deg @@ -73,7 +75,7 @@ template class MeasuredValues { constexpr bool has_angle_measurement(Idx bus) const { return !is_nan(imag(voltage(bus))); } constexpr bool has_bus_injection(Idx bus) const { return bus_injection_[bus].idx_bus_injection >= 0; } constexpr bool has_branch_from(Idx branch) const { - return idx_branch_from_power_[branch] >= 0 || idx_branch_from_current_[branch]; + return idx_branch_from_power_[branch] >= 0 || idx_branch_from_current_[branch] >= 0; } constexpr bool has_branch_to(Idx branch) const { return idx_branch_to_power_[branch] >= 0 || idx_branch_to_current_[branch] >= 0; @@ -445,29 +447,24 @@ template class MeasuredValues { template static VoltageSensorCalcParam combine_measurements(std::vector> const& data, IdxRange const& sensors) { - double accumulated_inverse_variance{}; - ComplexValue accumulated_value{}; - for (auto pos : sensors) { - auto const& measurement = data[pos]; - auto const inv_variance = 1.0 / measurement.variance; - - accumulated_inverse_variance += inv_variance; - if constexpr (only_magnitude) { - ComplexValue abs_value = piecewise_complex_value(DoubleComplex{0.0, nan}); - abs_value += detail::cabs_or_real(measurement.value); - accumulated_value += abs_value * inv_variance; - } else { - // accumulate value - accumulated_value += measurement.value * inv_variance; - } - } - - if (!is_normal(accumulated_inverse_variance)) { - return VoltageSensorCalcParam{accumulated_value, std::numeric_limits::infinity()}; + auto complex_measurements = sensors | std::views::transform([&](Idx pos) -> auto& { return data[pos]; }); + if constexpr (only_magnitude) { + auto const combined_magnitude_measurement = statistics::combine_measurements( + complex_measurements | std::views::transform([](auto const& measurement) { + return UniformRealRandVar{.value = detail::cabs_or_real(measurement.value), + .variance = measurement.variance}; + })); + return UniformComplexRandVar{.value = + [&combined_magnitude_measurement]() { + ComplexValue abs_value = + piecewise_complex_value(DoubleComplex{0.0, nan}); + abs_value += combined_magnitude_measurement.value; + return abs_value; + }(), + .variance = combined_magnitude_measurement.variance}; + } else { + return statistics::combine_measurements(complex_measurements); } - - return VoltageSensorCalcParam{accumulated_value / accumulated_inverse_variance, - 1.0 / accumulated_inverse_variance}; } template requires(!only_magnitude) @@ -480,8 +477,7 @@ template class MeasuredValues { requires(!only_magnitude) static CurrentSensorCalcParam combine_measurements(std::vector> const& data, IdxRange const& sensors) { - return {.angle_measurement_type = - data[sensors.front()].angle_measurement_type, // TODO(mgovers): this should be fixed + return {.angle_measurement_type = AngleMeasurementType::local, // TODO(mgovers): this should be fixed .measurement = statistics::combine_measurements( sensors | std::views::transform([&](Idx pos) -> auto& { return data[pos].measurement; }))}; } diff --git a/tests/cpp_unit_tests/test_observability.cpp b/tests/cpp_unit_tests/test_observability.cpp index 628f3907ab..a4f0571b29 100644 --- a/tests/cpp_unit_tests/test_observability.cpp +++ b/tests/cpp_unit_tests/test_observability.cpp @@ -121,11 +121,17 @@ TEST_CASE("Necessary observability check") { SUBCASE("Power sensor and current sensor are equivalent") { // remove the power sensor -> not observable topo.power_sensors_per_branch_from = {from_sparse, {0, 0, 0, 0}}; + se_input.measured_branch_from_power = {}; check_not_observable(topo, param, se_input); + // add the current sensor -> observable topo.current_sensors_per_branch_from = {from_sparse, {0, 1, 1, 1}}; + se_input.measured_branch_from_current = {{.angle_measurement_type = AngleMeasurementType::local, + .measurement = {.real_component = {.value = 1.0, .variance = 1.0}, + .imag_component = {.value = 0.0, .variance = 1.0}}}}; + check_observable(topo, param, se_input); - } + } // namespace power_grid_model } } // namespace power_grid_model From b14b5c144fafcc0adede1fdb326adb48b78a54cc Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Mon, 7 Apr 2025 09:27:22 +0200 Subject: [PATCH 07/25] local_angle, global_angle + add error handling current sensor Signed-off-by: Martijn Govers --- .../power_grid_model/common/exception.hpp | 41 +++++++++------- .../math_solver/measured_values.hpp | 26 ++++++++-- src/power_grid_model/_core/error_handling.py | 3 ++ src/power_grid_model/errors.py | 4 ++ tests/cpp_unit_tests/test_observability.cpp | 2 +- tests/unit/test_error_handling.py | 49 ++++++++++++++++++- 6 files changed, 102 insertions(+), 23 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/common/exception.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/common/exception.hpp index 084108a4be..a5ed7eec80 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/common/exception.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/common/exception.hpp @@ -39,13 +39,13 @@ class InvalidArguments : public PowerGridError { }; template ... Options> - InvalidArguments(std::string const& method, std::string const& arguments) { + InvalidArguments(std::string_view method, std::string_view arguments) { append_msg(std::format("{} is not implemented for {}!\n", method, arguments)); } template requires(std::same_as, TypeValuePair> && ...) - InvalidArguments(std::string const& method, Options const&... options) + InvalidArguments(std::string_view method, Options const&... options) : InvalidArguments{method, "the following combination of options"} { (append_msg(std::format(" {}: {}\n", options.name, options.value)), ...); } @@ -54,7 +54,7 @@ class InvalidArguments : public PowerGridError { class MissingCaseForEnumError : public InvalidArguments { public: template - MissingCaseForEnumError(std::string const& method, T const& value) + MissingCaseForEnumError(std::string_view method, T const& value) : InvalidArguments{method, std::format("{} #{}", typeid(T).name(), detail::to_string(static_cast(value)))} {} }; @@ -97,7 +97,7 @@ class InvalidTransformerClock : public PowerGridError { class SparseMatrixError : public PowerGridError { public: - SparseMatrixError(Idx err, std::string const& msg = "") { + SparseMatrixError(Idx err, std::string_view msg = "") { append_msg( std::format("Sparse matrix error with error code #{} (possibly singular)\n", detail::to_string(err))); if (!msg.empty()) { @@ -117,7 +117,7 @@ class SparseMatrixError : public PowerGridError { class NotObservableError : public PowerGridError { public: - NotObservableError(std::string const& msg = "") { + NotObservableError(std::string_view msg = "") { append_msg("Not enough measurements available for state estimation.\n"); if (!msg.empty()) { append_msg(std::format("{}\n", msg)); @@ -137,7 +137,7 @@ class IterationDiverge : public PowerGridError { class MaxIterationReached : public IterationDiverge { public: - MaxIterationReached(std::string const& msg = "") { + MaxIterationReached(std::string_view msg = "") { append_msg(std::format("Maximum number of iterations reached{}\n", msg)); } }; @@ -161,14 +161,14 @@ class Idx2DNotFound : public PowerGridError { class InvalidMeasuredObject : public PowerGridError { public: - InvalidMeasuredObject(std::string const& object, std::string const& sensor) { + InvalidMeasuredObject(std::string_view object, std::string_view sensor) { append_msg(std::format("{} measurement is not supported for object of type {}", sensor, object)); } }; class InvalidMeasuredTerminalType : public PowerGridError { public: - InvalidMeasuredTerminalType(MeasuredTerminalType const terminal_type, std::string const& sensor) { + InvalidMeasuredTerminalType(MeasuredTerminalType const terminal_type, std::string_view sensor) { append_msg(std::format("{} measurement is not supported for object of type {}", sensor, detail::to_string(static_cast(terminal_type)))); } @@ -176,10 +176,10 @@ class InvalidMeasuredTerminalType : public PowerGridError { class InvalidRegulatedObject : public PowerGridError { public: - InvalidRegulatedObject(std::string const& object, std::string const& regulator) { + InvalidRegulatedObject(std::string_view object, std::string_view regulator) { append_msg(std::format("{} regulator is not supported for object of type {}", regulator, object)); } - InvalidRegulatedObject(ID id, std::string const& regulator) { + InvalidRegulatedObject(ID id, std::string_view regulator) { append_msg( std::format("{} regulator is not supported for object with ID {}", regulator, detail::to_string(id))); } @@ -203,7 +203,7 @@ class AutomaticTapCalculationError : public PowerGridError { class AutomaticTapInputError : public PowerGridError { public: - AutomaticTapInputError(std::string const& msg) { + AutomaticTapInputError(std::string_view msg) { append_msg(std::format("Automatic tap changer has invalid configuration. {}", msg)); } }; @@ -215,14 +215,21 @@ class IDWrongType : public PowerGridError { } }; +class ConflictingAngleMeasurementType : public PowerGridError { + public: + ConflictingAngleMeasurementType(std::string_view msg) { + append_msg(std::format("Conflicting angle measurement type. {}", msg)); + } +}; + class CalculationError : public PowerGridError { public: - explicit CalculationError(std::string const& msg) { append_msg(msg); } + explicit CalculationError(std::string_view msg) { append_msg(msg); } }; class BatchCalculationError : public CalculationError { public: - BatchCalculationError(std::string const& msg, IdxVector failed_scenarios, std::vector err_msgs) + BatchCalculationError(std::string_view msg, IdxVector failed_scenarios, std::vector err_msgs) : CalculationError(msg), failed_scenarios_{std::move(failed_scenarios)}, err_msgs_(std::move(err_msgs)) {} IdxVector const& failed_scenarios() const { return failed_scenarios_; } @@ -270,12 +277,12 @@ class InvalidShortCircuitPhaseOrType : public PowerGridError { class SerializationError : public PowerGridError { public: - explicit SerializationError(std::string const& msg) { append_msg(msg); } + explicit SerializationError(std::string_view msg) { append_msg(msg); } }; class DatasetError : public PowerGridError { public: - explicit DatasetError(std::string const& msg) { append_msg(std::format("Dataset error: {}", msg)); } + explicit DatasetError(std::string_view msg) { append_msg(std::format("Dataset error: {}", msg)); } }; class ExperimentalFeature : public InvalidArguments { @@ -289,7 +296,7 @@ class NotImplementedError : public PowerGridError { class UnreachableHit : public PowerGridError { public: - UnreachableHit(std::string const& method, std::string const& reason_for_assumption) { + UnreachableHit(std::string_view method, std::string_view reason_for_assumption) { append_msg(std::format("Unreachable code hit when executing {}.\n The following assumption for unreachability " "was not met: {}.\n This may be a bug in the library\n", method, reason_for_assumption)); @@ -299,7 +306,7 @@ class UnreachableHit : public PowerGridError { class TapSearchStrategyIncompatibleError : public InvalidArguments { public: template - TapSearchStrategyIncompatibleError(std::string const& method, T1 const& value1, T2 const& value2) + TapSearchStrategyIncompatibleError(std::string_view method, T1 const& value1, T2 const& value2) : InvalidArguments{method, std::format("{} #{} and {} #{}", typeid(T1).name(), detail::to_string(static_cast(value1)), typeid(T2).name(), detail::to_string(static_cast(value2)))} {} diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp index 715d735d1b..0719dfab4f 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp @@ -9,6 +9,7 @@ Collect all measured Values */ #include "../calculation_parameters.hpp" +#include "../common/exception.hpp" #include "../common/three_phase_tensor.hpp" #include @@ -425,11 +426,11 @@ template class MeasuredValues { idx_branch_from_power_[branch] = process_one_object(branch, topo.power_sensors_per_branch_from, topo.branch_bus_idx, input.measured_branch_from_power, power_main_value_, branch_from_checker); - // to side current + // to side power idx_branch_to_power_[branch] = process_one_object(branch, topo.power_sensors_per_branch_to, topo.branch_bus_idx, input.measured_branch_to_power, power_main_value_, branch_to_checker); - // from side power + // from side current idx_branch_from_current_[branch] = process_one_object(branch, topo.current_sensors_per_branch_from, topo.branch_bus_idx, input.measured_branch_from_current, current_main_value_, branch_from_checker); @@ -437,6 +438,13 @@ template class MeasuredValues { idx_branch_to_current_[branch] = process_one_object(branch, topo.current_sensors_per_branch_to, topo.branch_bus_idx, input.measured_branch_to_current, current_main_value_, branch_to_checker); + + if ((!has_angle()) && std::ranges::any_of(current_main_value_, [](auto const& measurement) { + return measurement.angle_measurement_type == AngleMeasurementType::global_angle; + })) { + throw ConflictingAngleMeasurementType{ + "Global angle current measurements require a voltage angle measurement in the connected subgrid."}; + } } } @@ -477,9 +485,19 @@ template class MeasuredValues { requires(!only_magnitude) static CurrentSensorCalcParam combine_measurements(std::vector> const& data, IdxRange const& sensors) { - return {.angle_measurement_type = AngleMeasurementType::local, // TODO(mgovers): this should be fixed + auto const params = sensors | std::views::transform([&](Idx pos) -> auto& { return data[pos]; }); + auto const angle_measurement_type = sensors.empty() ? AngleMeasurementType::local_angle // fallback + : params.front().angle_measurement_type; + if (std::ranges::any_of(params, [angle_measurement_type](auto const& params) { + return params.angle_measurement_type != angle_measurement_type; + })) { + throw ConflictingAngleMeasurementType{ + "Cannot mix local and global angle current measurements on the same terminal."}; + } + + return {.angle_measurement_type = angle_measurement_type, .measurement = statistics::combine_measurements( - sensors | std::views::transform([&](Idx pos) -> auto& { return data[pos].measurement; }))}; + params | std::views::transform([&](auto const& params) -> auto& { return params.measurement; }))}; } template diff --git a/src/power_grid_model/_core/error_handling.py b/src/power_grid_model/_core/error_handling.py index 17d544fd01..17248926b3 100644 --- a/src/power_grid_model/_core/error_handling.py +++ b/src/power_grid_model/_core/error_handling.py @@ -16,6 +16,7 @@ AutomaticTapCalculationError, AutomaticTapInputError, ConflictID, + ConflictingAngleMeasurementType, ConflictVoltage, IDNotFound, IDWrongType, @@ -75,6 +76,7 @@ _AUTOMATIC_TAP_INPUT_ERROR_RE = re.compile(r"Automatic tap changer has invalid configuration") _ID_WRONG_TYPE_RE = re.compile(r"Wrong type for object with id (-?\d+)\n") +_CONFLICTING_ANGLE_MEASUREMENT_TYPE_RE = re.compile(r"Conflicting angle measurement type") _INVALID_CALCULATION_METHOD_RE = re.compile(r"The calculation method is invalid for this calculation!") _INVALID_SHORT_CIRCUIT_PHASE_OR_TYPE_RE = re.compile(r"short circuit type") # multiple different flavors _POWER_GRID_DATASET_ERROR_RE = re.compile(r"Dataset error: ") # multiple different flavors @@ -100,6 +102,7 @@ _AUTOMATIC_TAP_CALCULATION_ERROR_RE: AutomaticTapCalculationError, _AUTOMATIC_TAP_INPUT_ERROR_RE: AutomaticTapInputError, _ID_WRONG_TYPE_RE: IDWrongType, + _CONFLICTING_ANGLE_MEASUREMENT_TYPE_RE: ConflictingAngleMeasurementType, _INVALID_CALCULATION_METHOD_RE: InvalidCalculationMethod, _INVALID_SHORT_CIRCUIT_PHASE_OR_TYPE_RE: InvalidShortCircuitPhaseOrType, _POWER_GRID_DATASET_ERROR_RE: PowerGridDatasetError, diff --git a/src/power_grid_model/errors.py b/src/power_grid_model/errors.py index 4576253906..cbca6bf97c 100644 --- a/src/power_grid_model/errors.py +++ b/src/power_grid_model/errors.py @@ -101,6 +101,10 @@ class AutomaticTapInputError(PowerGridError): """Automatic tap changer has invalid configuration.""" +class ConflictingAngleMeasurementType(PowerGridError): + """Conflicting angle measurement types found.""" + + class InvalidShortCircuitPhaseOrType(PowerGridError): """Invalid (combination of) short circuit types and phase(s) provided.""" diff --git a/tests/cpp_unit_tests/test_observability.cpp b/tests/cpp_unit_tests/test_observability.cpp index a4f0571b29..efebf3b3ed 100644 --- a/tests/cpp_unit_tests/test_observability.cpp +++ b/tests/cpp_unit_tests/test_observability.cpp @@ -126,7 +126,7 @@ TEST_CASE("Necessary observability check") { // add the current sensor -> observable topo.current_sensors_per_branch_from = {from_sparse, {0, 1, 1, 1}}; - se_input.measured_branch_from_current = {{.angle_measurement_type = AngleMeasurementType::local, + se_input.measured_branch_from_current = {{.angle_measurement_type = AngleMeasurementType::local_angle, .measurement = {.real_component = {.value = 1.0, .variance = 1.0}, .imag_component = {.value = 0.0, .variance = 1.0}}}}; diff --git a/tests/unit/test_error_handling.py b/tests/unit/test_error_handling.py index b322757e5a..92ed97a322 100644 --- a/tests/unit/test_error_handling.py +++ b/tests/unit/test_error_handling.py @@ -8,10 +8,17 @@ from power_grid_model import PowerGridModel from power_grid_model._core.power_grid_meta import initialize_array -from power_grid_model.enum import CalculationMethod, LoadGenType, MeasuredTerminalType, TapChangingStrategy +from power_grid_model.enum import ( + CalculationMethod, + LoadGenType, + MeasuredTerminalType, + TapChangingStrategy, + _ExperimentalFeatures, +) from power_grid_model.errors import ( AutomaticTapInputError, ConflictID, + ConflictingAngleMeasurementType, ConflictVoltage, IDNotFound, IDWrongType, @@ -349,6 +356,46 @@ def test_automatic_tap_changing(): model.calculate_power_flow(tap_changing_strategy=TapChangingStrategy.min_voltage_tap) +def test_conflicting_angle_measurement_type(): + node_input = initialize_array("input", "node", 2) + node_input["id"] = [0, 1] + node_input["u_rated"] = [10e3, 10e3] + + line_input = initialize_array("input", "transformer", 1) + line_input["id"] = [2] + line_input["from_node"] = [0] + line_input["to_node"] = [1] + + source_input = initialize_array("input", "source", 1) + source_input["id"] = [3] + source_input["node"] = [0] + source_input["status"] = [1] + source_input["u_ref"] = [1.0] + + sym_current_sensor_input = initialize_array("input", "sym_current_sensor", 1) + sym_current_sensor_input["id"] = [4, 5] + sym_current_sensor_input["measured_object"] = [2, 2] + sym_current_sensor_input["measured_terminal_type"] = [MeasuredTerminalType.branch_from, MeasuredTerminalType.branch_from] + sym_current_sensor_input["angle_measurement_type"] = [AngleMeasurementType.local, AngleMeasurementType.global] + sym_current_sensor_input["status"] = [1, 1] + sym_current_sensor_input["i_sigma"] = [1.0, 1.0] + sym_current_sensor_input["i_angle_sigma"] = [1.0, 1.0] + sym_current_sensor_input["i_measured"] = [0.0, 0.0] + sym_current_sensor_input["i_angle_measured"] = [0.0, 0.0] + + model = PowerGridModel( + input_data={ + "node": node_input, + "line": line_input, + "source": source_input, + "sym_current_sensor": sym_current_sensor_input, + } + ) + + with pytest.raises(ConflictingAngleMeasurementType): + model._calculate_state_estimation(decode_error=True, experimental_features=_ExperimentalFeatures.enabled) + + @pytest.mark.skip(reason="TODO") def test_handle_power_grid_dataset_error(): pass From 2ed562c6e0386f453e11ffbe1878e7b5aac38fe1 Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Mon, 7 Apr 2025 12:44:30 +0200 Subject: [PATCH 08/25] test statistics combine measurements Signed-off-by: Martijn Govers --- .../power_grid_model/common/statistics.hpp | 16 +- .../math_solver/measured_values.hpp | 9 +- tests/cpp_unit_tests/test_statistics.cpp | 322 ++++++++++++++++++ 3 files changed, 335 insertions(+), 12 deletions(-) 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 784f81909c..239c981a68 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 @@ -209,13 +209,15 @@ template struct PolarComplexRandVar { namespace statistics { // combine multiple measurements of one quantity using Kalman filter -constexpr auto combine_measurements(std::ranges::view auto measurements) +constexpr auto accumulate(std::ranges::view auto measurements) requires std::same_as, UniformRealRandVar::sym>> || std::same_as, IndependentRealRandVar::sym>> || std::same_as, - UniformComplexRandVar::sym>> + UniformComplexRandVar::sym>> || + std::same_as, + IndependentComplexRandVar::sym>> { using RandVarType = std::ranges::range_value_t; using ValueType = decltype(RandVarType::value); @@ -236,17 +238,17 @@ constexpr auto combine_measurements(std::ranges::view auto measurements) } return RandVarType{.value = accumulated_value, .variance = VarianceType{std::numeric_limits::infinity()}}; } -constexpr auto combine_measurements(std::ranges::view auto measurements) +constexpr auto accumulate(std::ranges::view auto measurements) requires std::same_as, DecomposedComplexRandVar::sym>> { using sym = std::ranges::range_value_t::sym; DecomposedComplexRandVar result{ - .real_component = combine_measurements( - measurements | std::views::transform([](auto const& x) -> auto& { return x.real_component; })), - .imag_component = combine_measurements( - measurements | std::views::transform([](auto const& x) -> auto& { return x.imag_component; }))}; + .real_component = + accumulate(measurements | std::views::transform([](auto const& x) -> auto& { return x.real_component; })), + .imag_component = + accumulate(measurements | std::views::transform([](auto const& x) -> auto& { return x.imag_component; }))}; if (!(is_normal(result.real_component.variance) && is_normal(result.imag_component.variance))) { result.real_component.variance = RealValue{std::numeric_limits::infinity()}; diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp index 0719dfab4f..887baeafdb 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp @@ -457,7 +457,7 @@ template class MeasuredValues { IdxRange const& sensors) { auto complex_measurements = sensors | std::views::transform([&](Idx pos) -> auto& { return data[pos]; }); if constexpr (only_magnitude) { - auto const combined_magnitude_measurement = statistics::combine_measurements( + auto const combined_magnitude_measurement = statistics::accumulate( complex_measurements | std::views::transform([](auto const& measurement) { return UniformRealRandVar{.value = detail::cabs_or_real(measurement.value), .variance = measurement.variance}; @@ -471,15 +471,14 @@ template class MeasuredValues { }(), .variance = combined_magnitude_measurement.variance}; } else { - return statistics::combine_measurements(complex_measurements); + return statistics::accumulate(complex_measurements); } } template requires(!only_magnitude) static PowerSensorCalcParam combine_measurements(std::vector> const& data, IdxRange const& sensors) { - return statistics::combine_measurements(sensors | - std::views::transform([&](Idx pos) -> auto& { return data[pos]; })); + return statistics::accumulate(sensors | std::views::transform([&](Idx pos) -> auto& { return data[pos]; })); } template requires(!only_magnitude) @@ -496,7 +495,7 @@ template class MeasuredValues { } return {.angle_measurement_type = angle_measurement_type, - .measurement = statistics::combine_measurements( + .measurement = statistics::accumulate( params | std::views::transform([&](auto const& params) -> auto& { return params.measurement; }))}; } diff --git a/tests/cpp_unit_tests/test_statistics.cpp b/tests/cpp_unit_tests/test_statistics.cpp index 6a8462f73a..750a5f7e4e 100644 --- a/tests/cpp_unit_tests/test_statistics.cpp +++ b/tests/cpp_unit_tests/test_statistics.cpp @@ -815,6 +815,328 @@ TEST_CASE("Test statistics") { } } +TEST_CASE("Test statistics - accumulate") { + using statistics::accumulate; + using std::views::take; + + SUBCASE("UniformRealRandVar | IndependentRealRandVar") { + // using a template lambda to avoid code duplication and to avoid having to create a separate test case + auto const check = []() { + std::vector const measurements{ + {.value = 1.0, .variance = 0.2}, {.value = 2.0, .variance = 0.3}, {.value = 3.0, .variance = 0.6}}; + + CHECK(accumulate(measurements | take(0)).value == 0.0); + CHECK(is_inf(accumulate(measurements | take(0)).variance)); + + CHECK(accumulate(measurements | take(1)).value == measurements.front().value); + CHECK(accumulate(measurements | take(1)).variance == measurements.front().variance); + + CHECK(accumulate(measurements | take(2)).value == doctest::Approx(7.0 / 5.0)); + CHECK(accumulate(measurements | take(2)).variance == doctest::Approx(3.0 / 25.0)); + + CHECK(accumulate(measurements | take(3)).value == doctest::Approx(5.0 / 3.0)); + CHECK(accumulate(measurements | take(3)).variance == doctest::Approx(1.0 / 10.0)); + }; + SUBCASE("UniformRealRandVar") { check.template operator()>(); } + SUBCASE("IndependentRealRandVar") { + check.template operator()>(); + } + } + + SUBCASE("UniformRealRandVar") { + std::vector> const measurements{{.value = {1.0, 2.0, -1.0}, .variance = 0.2}, + {.value = {2.0, 4.0, 3.0}, .variance = 0.3}, + {.value = {4.0, 5.0, 6.0}, .variance = 0.6}}; + + CHECK(accumulate(measurements | take(0)).value(0) == 0.0); + CHECK(accumulate(measurements | take(0)).value(1) == 0.0); + CHECK(accumulate(measurements | take(0)).value(2) == 0.0); + CHECK(is_inf(accumulate(measurements | take(0)).variance)); + + CHECK(accumulate(measurements | take(1)).value(0) == measurements.front().value(0)); + CHECK(accumulate(measurements | take(1)).value(1) == measurements.front().value(1)); + CHECK(accumulate(measurements | take(1)).value(2) == measurements.front().value(2)); + CHECK(accumulate(measurements | take(1)).variance == measurements.front().variance); + + CHECK(accumulate(measurements | take(2)).value(0) == doctest::Approx(7.0 / 5.0)); + CHECK(accumulate(measurements | take(2)).value(1) == doctest::Approx(14.0 / 5.0)); + CHECK(accumulate(measurements | take(2)).value(2) == doctest::Approx(3.0 / 5.0)); + CHECK(accumulate(measurements | take(2)).variance == doctest::Approx(3.0 / 25.0)); + + CHECK(accumulate(measurements | take(3)).value(0) == doctest::Approx(11.0 / 6.0)); + CHECK(accumulate(measurements | take(3)).value(1) == doctest::Approx(19.0 / 6.0)); + CHECK(accumulate(measurements | take(3)).value(2) == doctest::Approx(3.0 / 2.0)); + CHECK(accumulate(measurements | take(3)).variance == doctest::Approx(1.0 / 10.0)); + } + + SUBCASE("IndependentRealRandVar") { + std::vector> const measurements{ + {.value = {1.0, 2.0, -1.0}, .variance = {0.2, 0.3, 0.4}}, + {.value = {2.0, 4.0, 3.0}, .variance = {0.3, 0.4, 0.5}}, + {.value = {4.0, 5.0, 6.0}, .variance = {0.6, 0.7, 0.8}}}; + + CHECK(accumulate(measurements | take(0)).value(0) == 0.0); + CHECK(accumulate(measurements | take(0)).value(1) == 0.0); + CHECK(accumulate(measurements | take(0)).value(2) == 0.0); + CHECK(is_inf(accumulate(measurements | take(0)).variance(0))); + CHECK(is_inf(accumulate(measurements | take(0)).variance(1))); + CHECK(is_inf(accumulate(measurements | take(0)).variance(2))); + + CHECK(accumulate(measurements | take(1)).value(0) == measurements.front().value(0)); + CHECK(accumulate(measurements | take(1)).value(1) == measurements.front().value(1)); + CHECK(accumulate(measurements | take(1)).value(2) == measurements.front().value(2)); + CHECK(accumulate(measurements | take(1)).variance(0) == measurements.front().variance(0)); + CHECK(accumulate(measurements | take(1)).variance(1) == measurements.front().variance(1)); + CHECK(accumulate(measurements | take(1)).variance(2) == measurements.front().variance(2)); + + CHECK(accumulate(measurements | take(2)).value(0) == doctest::Approx(7.0 / 5.0)); + CHECK(accumulate(measurements | take(2)).value(1) == doctest::Approx(20.0 / 7.0)); + CHECK(accumulate(measurements | take(2)).value(2) == doctest::Approx(7.0 / 9.0)); + CHECK(accumulate(measurements | take(2)).variance(0) == doctest::Approx(3.0 / 25.0)); + CHECK(accumulate(measurements | take(2)).variance(1) == doctest::Approx(6.0 / 35.0)); + CHECK(accumulate(measurements | take(2)).variance(2) == doctest::Approx(2.0 / 9.0)); + + CHECK(accumulate(measurements | take(3)).value(0) == doctest::Approx(11.0 / 6.0)); + CHECK(accumulate(measurements | take(3)).value(1) == doctest::Approx(200.0 / 61.0)); + CHECK(accumulate(measurements | take(3)).value(2) == doctest::Approx(44.0 / 23.0)); + CHECK(accumulate(measurements | take(3)).variance(0) == doctest::Approx(1.0 / 10.0)); + CHECK(accumulate(measurements | take(3)).variance(1) == doctest::Approx(42.0 / 305.0)); + CHECK(accumulate(measurements | take(3)).variance(2) == doctest::Approx(4.0 / 23.0)); + } + + SUBCASE("UniformComplexRandVar | IndependentComplexRandVar") { + // using a template lambda to avoid code duplication and to avoid having to create a separate test case + auto const check = []() { + std::vector const measurements{T{.value = 1.0 + 5.0i, .variance = 0.2}, + T{.value = 2.0 + 6.0i, .variance = 0.3}, + T{.value = 4.0 + 3.0i, .variance = 0.6}}; + + CHECK(accumulate(measurements | take(0)).value.real() == 0.0); + CHECK(accumulate(measurements | take(0)).value.imag() == 0.0); + CHECK(is_inf(accumulate(measurements | take(0)).variance)); + + CHECK(accumulate(measurements | take(1)).value.real() == measurements.front().value.real()); + CHECK(accumulate(measurements | take(1)).value.imag() == measurements.front().value.imag()); + CHECK(accumulate(measurements | take(1)).variance == measurements.front().variance); + + CHECK(accumulate(measurements | take(2)).value.real() == doctest::Approx(7.0 / 5.0)); + CHECK(accumulate(measurements | take(2)).value.imag() == doctest::Approx(27.0 / 5.0)); + CHECK(accumulate(measurements | take(2)).variance == doctest::Approx(3.0 / 25.0)); + + CHECK(accumulate(measurements | take(3)).value.real() == doctest::Approx(11.0 / 6.0)); + CHECK(accumulate(measurements | take(3)).value.imag() == doctest::Approx(30.0 / 6.0)); + CHECK(accumulate(measurements | take(3)).variance == doctest::Approx(1.0 / 10.0)); + }; + SUBCASE("UniformComplexRandVar") { + check.template operator()>(); + } + SUBCASE("IndependentComplexRandVar") { + check.template operator()>(); + } + } + + SUBCASE("UniformComplexRandVar") { + std::vector> const measurements{ + {.value = {RealValue{1.0, 2.0, -1.0}, RealValue{5.0, 6.0, 7.0}}, + .variance = 0.2}, + {.value = {RealValue{2.0, 4.0, 3.0}, RealValue{6.0, -7.0, 2.0}}, + .variance = 0.3}, + {.value = {RealValue{4.0, 5.0, 6.0}, RealValue{3.0, 1.0, 2.0}}, + .variance = 0.6}}; + + CHECK(accumulate(measurements | take(0)).value(0).real() == 0.0); + CHECK(accumulate(measurements | take(0)).value(1).real() == 0.0); + CHECK(accumulate(measurements | take(0)).value(2).real() == 0.0); + CHECK(accumulate(measurements | take(0)).value(0).imag() == 0.0); + CHECK(accumulate(measurements | take(0)).value(1).imag() == 0.0); + CHECK(accumulate(measurements | take(0)).value(2).imag() == 0.0); + CHECK(is_inf(accumulate(measurements | take(0)).variance)); + + CHECK(accumulate(measurements | take(1)).value(0).real() == measurements.front().value(0).real()); + CHECK(accumulate(measurements | take(1)).value(1).real() == measurements.front().value(1).real()); + CHECK(accumulate(measurements | take(1)).value(2).real() == measurements.front().value(2).real()); + CHECK(accumulate(measurements | take(1)).value(0).imag() == measurements.front().value(0).imag()); + CHECK(accumulate(measurements | take(1)).value(1).imag() == measurements.front().value(1).imag()); + CHECK(accumulate(measurements | take(1)).value(2).imag() == measurements.front().value(2).imag()); + CHECK(accumulate(measurements | take(1)).variance == measurements.front().variance); + + CHECK(accumulate(measurements | take(2)).value(0).real() == doctest::Approx(7.0 / 5.0)); + CHECK(accumulate(measurements | take(2)).value(1).real() == doctest::Approx(14.0 / 5.0)); + CHECK(accumulate(measurements | take(2)).value(2).real() == doctest::Approx(3.0 / 5.0)); + CHECK(accumulate(measurements | take(2)).value(0).imag() == doctest::Approx(27.0 / 5.0)); + CHECK(accumulate(measurements | take(2)).value(1).imag() == doctest::Approx(4.0 / 5.0)); + CHECK(accumulate(measurements | take(2)).value(2).imag() == doctest::Approx(25.0 / 5.0)); + CHECK(accumulate(measurements | take(2)).variance == doctest::Approx(3.0 / 25.0)); + + CHECK(accumulate(measurements | take(3)).value(0).real() == doctest::Approx(11.0 / 6.0)); + CHECK(accumulate(measurements | take(3)).value(1).real() == doctest::Approx(19.0 / 6.0)); + CHECK(accumulate(measurements | take(3)).value(2).real() == doctest::Approx(9.0 / 6.0)); + CHECK(accumulate(measurements | take(3)).value(0).imag() == doctest::Approx(30.0 / 6.0)); + CHECK(accumulate(measurements | take(3)).value(1).imag() == doctest::Approx(5.0 / 6.0)); + CHECK(accumulate(measurements | take(3)).value(2).imag() == doctest::Approx(27.0 / 6.0)); + CHECK(accumulate(measurements | take(3)).variance == doctest::Approx(1.0 / 10.0)); + } + + SUBCASE("IndependentComplexRandVar") { + std::vector> const measurements{ + {.value = {RealValue{1.0, 2.0, -1.0}, RealValue{5.0, 6.0, 7.0}}, + .variance = {0.2, 0.3, 0.4}}, + {.value = {RealValue{2.0, 4.0, 3.0}, RealValue{6.0, -7.0, 2.0}}, + .variance = {0.3, 0.4, 0.5}}, + {.value = {RealValue{4.0, 5.0, 6.0}, RealValue{3.0, 1.0, 2.0}}, + .variance = {0.6, 0.7, 0.8}}}; + + CHECK(accumulate(measurements | take(0)).value(0).real() == 0.0); + CHECK(accumulate(measurements | take(0)).value(1).real() == 0.0); + CHECK(accumulate(measurements | take(0)).value(2).real() == 0.0); + CHECK(accumulate(measurements | take(0)).value(0).imag() == 0.0); + CHECK(accumulate(measurements | take(0)).value(1).imag() == 0.0); + CHECK(accumulate(measurements | take(0)).value(2).imag() == 0.0); + CHECK(is_inf(accumulate(measurements | take(0)).variance(0))); + CHECK(is_inf(accumulate(measurements | take(0)).variance(1))); + CHECK(is_inf(accumulate(measurements | take(0)).variance(2))); + + CHECK(accumulate(measurements | take(1)).value(0).real() == measurements.front().value(0).real()); + CHECK(accumulate(measurements | take(1)).value(1).real() == measurements.front().value(1).real()); + CHECK(accumulate(measurements | take(1)).value(2).real() == measurements.front().value(2).real()); + CHECK(accumulate(measurements | take(1)).value(0).imag() == measurements.front().value(0).imag()); + CHECK(accumulate(measurements | take(1)).value(1).imag() == measurements.front().value(1).imag()); + CHECK(accumulate(measurements | take(1)).value(2).imag() == measurements.front().value(2).imag()); + CHECK(accumulate(measurements | take(1)).variance(0) == measurements.front().variance(0)); + CHECK(accumulate(measurements | take(1)).variance(1) == measurements.front().variance(1)); + CHECK(accumulate(measurements | take(1)).variance(2) == measurements.front().variance(2)); + + CHECK(accumulate(measurements | take(2)).value(0).real() == doctest::Approx(7.0 / 5.0)); + CHECK(accumulate(measurements | take(2)).value(1).real() == doctest::Approx(20.0 / 7.0)); + CHECK(accumulate(measurements | take(2)).value(2).real() == doctest::Approx(7.0 / 9.0)); + CHECK(accumulate(measurements | take(2)).value(0).imag() == doctest::Approx(27.0 / 5.0)); + CHECK(accumulate(measurements | take(2)).value(1).imag() == doctest::Approx(3.0 / 7.0)); + CHECK(accumulate(measurements | take(2)).value(2).imag() == doctest::Approx(43.0 / 9.0)); + CHECK(accumulate(measurements | take(2)).variance(0) == doctest::Approx(3.0 / 25.0)); + CHECK(accumulate(measurements | take(2)).variance(1) == doctest::Approx(6.0 / 35.0)); + CHECK(accumulate(measurements | take(2)).variance(2) == doctest::Approx(2.0 / 9.0)); + + CHECK(accumulate(measurements | take(3)).value(0).real() == doctest::Approx(11.0 / 6.0)); + CHECK(accumulate(measurements | take(3)).value(1).real() == doctest::Approx(200.0 / 61.0)); + CHECK(accumulate(measurements | take(3)).value(2).real() == doctest::Approx(44.0 / 23.0)); + CHECK(accumulate(measurements | take(3)).value(0).imag() == doctest::Approx(30.0 / 6.0)); + CHECK(accumulate(measurements | take(3)).value(1).imag() == doctest::Approx(33.0 / 61.0)); + CHECK(accumulate(measurements | take(3)).value(2).imag() == doctest::Approx(96.0 / 23.0)); + CHECK(accumulate(measurements | take(3)).variance(0) == doctest::Approx(1.0 / 10.0)); + CHECK(accumulate(measurements | take(3)).variance(1) == doctest::Approx(42.0 / 305.0)); + CHECK(accumulate(measurements | take(3)).variance(2) == doctest::Approx(4.0 / 23.0)); + } + + SUBCASE("DecomposedComplexRandVar") { + std::vector> const measurements{ + DecomposedComplexRandVar{.real_component = {.value = 1.0, .variance = 0.2}, + .imag_component = {.value = 5.0, .variance = 0.1}}, + DecomposedComplexRandVar{.real_component = {.value = 2.0, .variance = 0.3}, + .imag_component = {.value = 6.0, .variance = 0.2}}, + DecomposedComplexRandVar{.real_component = {.value = 4.0, .variance = 0.6}, + .imag_component = {.value = 3.0, .variance = 0.3}}}; + + CHECK(accumulate(measurements | take(0)).real_component.value == 0.0); + CHECK(accumulate(measurements | take(0)).imag_component.value == 0.0); + CHECK(is_inf(accumulate(measurements | take(0)).real_component.variance)); + CHECK(is_inf(accumulate(measurements | take(0)).imag_component.variance)); + + CHECK(accumulate(measurements | take(1)).real_component.value == measurements.front().real_component.value); + CHECK(accumulate(measurements | take(1)).imag_component.value == measurements.front().imag_component.value); + CHECK(accumulate(measurements | take(1)).real_component.variance == + measurements.front().real_component.variance); + CHECK(accumulate(measurements | take(1)).imag_component.variance == + measurements.front().imag_component.variance); + + CHECK(accumulate(measurements | take(2)).real_component.value == doctest::Approx(7.0 / 5.0)); + CHECK(accumulate(measurements | take(2)).imag_component.value == doctest::Approx(80.0 / 15.0)); + CHECK(accumulate(measurements | take(2)).real_component.variance == doctest::Approx(3.0 / 25.0)); + CHECK(accumulate(measurements | take(2)).imag_component.variance == doctest::Approx(1.0 / 15.0)); + + CHECK(accumulate(measurements | take(3)).real_component.value == doctest::Approx(11.0 / 6.0)); + CHECK(accumulate(measurements | take(3)).imag_component.value == doctest::Approx(270.0 / 55.0)); + CHECK(accumulate(measurements | take(3)).real_component.variance == doctest::Approx(1.0 / 10.0)); + CHECK(accumulate(measurements | take(3)).imag_component.variance == doctest::Approx(3.0 / 55.0)); + } + + SUBCASE("DecomposedComplexRandVar") { + std::vector> const measurements{ + DecomposedComplexRandVar{ + .real_component = {.value = RealValue{1.0, 2.0, -1.0}, .variance = {0.2, 0.3, 0.4}}, + .imag_component = {.value = RealValue{5.0, 6.0, 7.0}, .variance = {0.1, 0.2, 0.3}}}, + DecomposedComplexRandVar{ + .real_component = {.value = RealValue{2.0, 4.0, 3.0}, .variance = {0.3, 0.4, 0.5}}, + .imag_component = {.value = RealValue{6.0, -7.0, 2.0}, .variance = {0.2, 0.3, 0.4}}}, + DecomposedComplexRandVar{ + .real_component = {.value = RealValue{4.0, 5.0, 6.0}, .variance = {0.6, 0.7, 0.8}}, + .imag_component = {.value = RealValue{3.0, 1.0, 2.0}, .variance = {0.3, 0.4, 0.5}}}}; + + CHECK(accumulate(measurements | take(0)).real_component.value(0) == 0.0); + CHECK(accumulate(measurements | take(0)).real_component.value(2) == 0.0); + CHECK(accumulate(measurements | take(0)).real_component.value(1) == 0.0); + CHECK(accumulate(measurements | take(0)).imag_component.value(0) == 0.0); + CHECK(accumulate(measurements | take(0)).imag_component.value(2) == 0.0); + CHECK(accumulate(measurements | take(0)).imag_component.value(1) == 0.0); + CHECK(is_inf(accumulate(measurements | take(0)).real_component.variance(0))); + CHECK(is_inf(accumulate(measurements | take(0)).real_component.variance(1))); + CHECK(is_inf(accumulate(measurements | take(0)).real_component.variance(2))); + CHECK(is_inf(accumulate(measurements | take(0)).imag_component.variance(0))); + CHECK(is_inf(accumulate(measurements | take(0)).imag_component.variance(1))); + CHECK(is_inf(accumulate(measurements | take(0)).imag_component.variance(2))); + + CHECK(accumulate(measurements | take(1)).real_component.value(0) == + measurements.front().real_component.value(0)); + CHECK(accumulate(measurements | take(1)).real_component.value(1) == + measurements.front().real_component.value(1)); + CHECK(accumulate(measurements | take(1)).real_component.value(2) == + measurements.front().real_component.value(2)); + CHECK(accumulate(measurements | take(1)).imag_component.value(0) == + measurements.front().imag_component.value(0)); + CHECK(accumulate(measurements | take(1)).imag_component.value(1) == + measurements.front().imag_component.value(1)); + CHECK(accumulate(measurements | take(1)).imag_component.value(2) == + measurements.front().imag_component.value(2)); + CHECK(accumulate(measurements | take(1)).real_component.variance(0) == + measurements.front().real_component.variance(0)); + CHECK(accumulate(measurements | take(1)).real_component.variance(1) == + measurements.front().real_component.variance(1)); + CHECK(accumulate(measurements | take(1)).real_component.variance(2) == + measurements.front().real_component.variance(2)); + CHECK(accumulate(measurements | take(1)).imag_component.variance(0) == + measurements.front().imag_component.variance(0)); + CHECK(accumulate(measurements | take(1)).imag_component.variance(1) == + measurements.front().imag_component.variance(1)); + CHECK(accumulate(measurements | take(1)).imag_component.variance(2) == + measurements.front().imag_component.variance(2)); + + CHECK(accumulate(measurements | take(2)).real_component.value(0) == doctest::Approx(7.0 / 5.0)); + CHECK(accumulate(measurements | take(2)).real_component.value(1) == doctest::Approx(20.0 / 7.0)); + CHECK(accumulate(measurements | take(2)).real_component.value(2) == doctest::Approx(7.0 / 9.0)); + CHECK(accumulate(measurements | take(2)).imag_component.value(0) == doctest::Approx(80.0 / 15.0)); + CHECK(accumulate(measurements | take(2)).imag_component.value(1) == doctest::Approx(20 / 25.0)); + CHECK(accumulate(measurements | take(2)).imag_component.value(2) == doctest::Approx(170.0 / 35.0)); + CHECK(accumulate(measurements | take(2)).real_component.variance(0) == doctest::Approx(3.0 / 25.0)); + CHECK(accumulate(measurements | take(2)).real_component.variance(1) == doctest::Approx(6.0 / 35.0)); + CHECK(accumulate(measurements | take(2)).real_component.variance(2) == doctest::Approx(2.0 / 9.0)); + CHECK(accumulate(measurements | take(2)).imag_component.variance(0) == doctest::Approx(1.0 / 15.0)); + CHECK(accumulate(measurements | take(2)).imag_component.variance(1) == doctest::Approx(3.0 / 25.0)); + CHECK(accumulate(measurements | take(2)).imag_component.variance(2) == doctest::Approx(6.0 / 35.0)); + + CHECK(accumulate(measurements | take(3)).real_component.value(0) == doctest::Approx(11.0 / 6.0)); + CHECK(accumulate(measurements | take(3)).real_component.value(1) == doctest::Approx(200.0 / 61.0)); + CHECK(accumulate(measurements | take(3)).real_component.value(2) == doctest::Approx(44.0 / 23.0)); + CHECK(accumulate(measurements | take(3)).imag_component.value(0) == doctest::Approx(270.0 / 55.0)); + CHECK(accumulate(measurements | take(3)).imag_component.value(1) == doctest::Approx(55.0 / 65.0)); + CHECK(accumulate(measurements | take(3)).imag_component.value(2) == doctest::Approx(194.0 / 47.0)); + CHECK(accumulate(measurements | take(3)).real_component.variance(0) == doctest::Approx(1.0 / 10.0)); + CHECK(accumulate(measurements | take(3)).real_component.variance(1) == doctest::Approx(42.0 / 305.0)); + CHECK(accumulate(measurements | take(3)).real_component.variance(2) == doctest::Approx(4.0 / 23.0)); + CHECK(accumulate(measurements | take(3)).imag_component.variance(0) == doctest::Approx(3.0 / 55.0)); + CHECK(accumulate(measurements | take(3)).imag_component.variance(1) == doctest::Approx(6.0 / 65)); + CHECK(accumulate(measurements | take(3)).imag_component.variance(2) == doctest::Approx(6.0 / 47.0)); + } +} + TEST_SUITE_END(); } // namespace power_grid_model From 0dea4712f3dfde736474b77ff9e1fccf0edb5f72 Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Mon, 7 Apr 2025 15:46:13 +0200 Subject: [PATCH 09/25] measured values unit tests Signed-off-by: Martijn Govers --- .../iterative_linear_se_solver.hpp | 7 +- .../math_solver/measured_values.hpp | 24 ++- .../math_solver/newton_raphson_se_solver.hpp | 4 +- .../math_solver/observability.hpp | 3 +- tests/cpp_unit_tests/test_measured_values.cpp | 140 +++++++++++++++++- 5 files changed, 161 insertions(+), 17 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp index f99003a5d4..c10d31c1e5 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp @@ -142,7 +142,8 @@ template class IterativeLinearSESolver { private: // array selection function pointer - static constexpr std::array has_branch_{&MeasuredValues::has_branch_from, &MeasuredValues::has_branch_to}; + static constexpr std::array has_branch_power_{&MeasuredValues::has_branch_from_power, + &MeasuredValues::has_branch_to_power}; static constexpr std::array branch_power_{&MeasuredValues::branch_from_power, &MeasuredValues::branch_to_power}; @@ -210,7 +211,7 @@ template class IterativeLinearSESolver { // measured at from-side: 0, to-side: 1 for (IntS const measured_side : std::array{0, 1}) { // has measurement - if (std::invoke(has_branch_[measured_side], measured_value, obj)) { + if (std::invoke(has_branch_power_[measured_side], measured_value, obj)) { // G += Y{side, b0}^H * (variance^-1) * Y{side, b1} auto const& power = std::invoke(branch_power_[measured_side], measured_value, obj); block.g() += @@ -301,7 +302,7 @@ template class IterativeLinearSESolver { // measured at from-side: 0, to-side: 1 for (IntS const measured_side : std::array{0, 1}) { // has measurement - if (std::invoke(has_branch_[measured_side], measured_value, obj)) { + if (std::invoke(has_branch_power_[measured_side], measured_value, obj)) { PowerSensorCalcParam const& m = std::invoke(branch_power_[measured_side], measured_value, obj); // the current needs to be calculated with the voltage of the measured bus side diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp index 887baeafdb..da09a1ab27 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp @@ -75,12 +75,10 @@ template class MeasuredValues { constexpr bool has_voltage(Idx bus) const { return idx_voltage_[bus] >= 0; } constexpr bool has_angle_measurement(Idx bus) const { return !is_nan(imag(voltage(bus))); } constexpr bool has_bus_injection(Idx bus) const { return bus_injection_[bus].idx_bus_injection >= 0; } - constexpr bool has_branch_from(Idx branch) const { - return idx_branch_from_power_[branch] >= 0 || idx_branch_from_current_[branch] >= 0; - } - constexpr bool has_branch_to(Idx branch) const { - return idx_branch_to_power_[branch] >= 0 || idx_branch_to_current_[branch] >= 0; - } + constexpr bool has_branch_from_power(Idx branch) const { return idx_branch_from_power_[branch] >= 0; } + constexpr bool has_branch_to_power(Idx branch) const { return idx_branch_to_power_[branch] >= 0; } + constexpr bool has_branch_from_current(Idx branch) const { return idx_branch_from_current_[branch] >= 0; } + constexpr bool has_branch_to_current(Idx branch) const { return idx_branch_to_current_[branch] >= 0; } constexpr bool has_shunt(Idx shunt) const { return idx_shunt_power_[shunt] >= 0; } constexpr bool has_load_gen(Idx load_gen) const { return idx_load_gen_power_[load_gen] >= 0; } constexpr bool has_source(Idx source) const { return idx_source_power_[source] >= 0; } @@ -564,6 +562,16 @@ template class MeasuredValues { } } } + for (auto const& x : current_main_value_) { + auto const variance = x.measurement.real_component.variance + x.measurement.imag_component.variance; + if constexpr (is_symmetric_v) { + unconstrained_min(variance); + } else { + for (Idx const phase : {0, 1, 2}) { + unconstrained_min(variance[phase]); + } + } + } // scale auto const inv_norm_var = 1.0 / min_var; @@ -572,6 +580,10 @@ template class MeasuredValues { x.real_component.variance *= inv_norm_var; x.imag_component.variance *= inv_norm_var; }); + std::ranges::for_each(current_main_value_, [inv_norm_var](auto& x) { + x.measurement.real_component.variance *= inv_norm_var; + x.measurement.imag_component.variance *= inv_norm_var; + }); } void calculate_non_over_determined_injection(Idx n_unmeasured, IdxRange const& load_gens, IdxRange const& sources, diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp index 809734124a..faf69b16b3 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp @@ -304,14 +304,14 @@ template class NewtonRaphsonSESolver { [[fallthrough]]; case YBusElementType::btf: { auto const& y_branch = param.branch_param[obj]; - if (measured_values.has_branch_from(obj)) { + if (measured_values.has_branch_from_power(obj)) { auto const ij_voltage_order = (type == YBusElementType::bft) ? Order::row_major : Order::column_major; process_branch_measurement(block, diag_block, rhs_block, y_branch.yff(), y_branch.yft(), u_state, ij_voltage_order, measured_values.branch_from_power(obj)); } - if (measured_values.has_branch_to(obj)) { + if (measured_values.has_branch_to_power(obj)) { auto const ij_voltage_order = (type == YBusElementType::btf) ? Order::row_major : Order::column_major; process_branch_measurement(block, diag_block, rhs_block, y_branch.ytt(), y_branch.ytf(), diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp index 7e855f5d95..69c85bcbb6 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp @@ -63,7 +63,8 @@ std::pair, bool> count_flow_sensors(MeasuredValues cons // we only need one flow sensor, so the loop will break if (element.element_type != YBusElementType::shunt) { Idx const branch = element.idx; - if ((measured_values.has_branch_from(branch) || measured_values.has_branch_to(branch)) && + if ((measured_values.has_branch_from_power(branch) || + measured_values.has_branch_to_power(branch)) && topo.branch_bus_idx[branch][0] != -1 && topo.branch_bus_idx[branch][1] != -1) { flow_sensors[ybus_index] = 1; has_at_least_one_sensor = true; diff --git a/tests/cpp_unit_tests/test_measured_values.cpp b/tests/cpp_unit_tests/test_measured_values.cpp index d59e7e87e2..e265058c58 100644 --- a/tests/cpp_unit_tests/test_measured_values.cpp +++ b/tests/cpp_unit_tests/test_measured_values.cpp @@ -39,7 +39,7 @@ TEST_CASE("Measured Values") { MeasuredValues const values{std::make_shared(std::move(topo)), input}; - CHECK(values.has_bus_injection(0)); + REQUIRE(values.has_bus_injection(0)); auto const& injection = values.bus_injection(0); check_close(injection.value(), 1.0 + 0.1i); check_close(injection.real_component.variance, 0.75); @@ -61,7 +61,7 @@ TEST_CASE("Measured Values") { MeasuredValues const values{std::make_shared(std::move(topo)), input}; - CHECK(values.has_bus_injection(0)); + REQUIRE(values.has_bus_injection(0)); auto const& injection = values.bus_injection(0); check_close(injection.value(), ComplexValue{1.0 + 0.0i, 1.1 + 0.1i, 1.2 - 0.2i}); check_close(injection.real_component.variance, RealValue{0.75, 1.5, 1.625}); @@ -84,7 +84,7 @@ TEST_CASE("Measured Values") { MeasuredValues const values{std::make_shared(std::move(topo)), input}; - CHECK(values.has_bus_injection(0)); + REQUIRE(values.has_bus_injection(0)); auto const& injection = values.bus_injection(0); check_close(injection.value(), 2.0 + 1.0i); check_close(injection.real_component.variance, 16.0 / 61.0); @@ -107,7 +107,7 @@ TEST_CASE("Measured Values") { MeasuredValues const values{std::make_shared(std::move(topo)), input}; - CHECK(values.has_bus_injection(0)); + REQUIRE(values.has_bus_injection(0)); auto const& injection = values.bus_injection(0); check_close(injection.value(), 5.0 + 2.2i); check_close(injection.real_component.variance, 0.3); @@ -131,12 +131,142 @@ TEST_CASE("Measured Values") { MeasuredValues const values{std::make_shared(std::move(topo)), input}; - CHECK(values.has_bus_injection(0)); + REQUIRE(values.has_bus_injection(0)); auto const& injection = values.bus_injection(0); check_close(injection.value(), ComplexValue{5.0 + 2.2i, 0.1 + 0.2i, 0.8 + 3.0i}); check_close(injection.real_component.variance, RealValue{0.3, 0.2, 0.85}); check_close(injection.imag_component.variance, RealValue{0.7, 0.85, 0.2}); } + + SUBCASE("Accumulate branch flow sensors") { + using enum AngleMeasurementType; + + auto const get_single_branch_topology = [] { + auto topo = MathModelTopology{}; + topo.phase_shift = {0.0, 0.0}; + topo.branch_bus_idx = {{0, 1}}; + topo.shunts_per_bus = {from_dense, {}, 1}; + topo.load_gens_per_bus = {from_dense, {}, 1}; + topo.sources_per_bus = {from_dense, {}, 1}; + return topo; + }; + DecomposedComplexRandVar const measurement_a{.real_component = {.value = 1.0, .variance = 1.0}, + .imag_component = {.value = 1.5, .variance = 4.0}}; + DecomposedComplexRandVar const measurement_b{.real_component = {.value = 4.0, .variance = 2.0}, + .imag_component = {.value = 0.7, .variance = 3.0}}; + auto const check_accumulated = [](DecomposedComplexRandVar const& value) { + check_close(value.value(), 2.0 + (73.0i / 70.0)); + check_close(value.real_component.variance, 7.0 / 25.0); + check_close(value.imag_component.variance, 18.0 / 25.0); + }; + + SUBCASE("Accumulate power sensors on branch - sym") { + auto topo = get_single_branch_topology(); + topo.power_sensors_per_branch_from = {from_dense, {0, 0}, 1}; + + StateEstimationInput input{}; + input.measured_branch_from_power = {measurement_a, measurement_b}; + + MeasuredValues const values{std::make_shared(std::move(topo)), input}; + + REQUIRE(values.has_branch_from_power(0)); + CHECK_FALSE(values.has_branch_from_current(0)); + CHECK_FALSE(values.has_branch_to_power(0)); + CHECK_FALSE(values.has_branch_to_current(0)); + + check_accumulated(values.branch_from_power(0)); + } + + SUBCASE("Accumulate local angle current sensors on branch - sym") { + auto topo = get_single_branch_topology(); + topo.current_sensors_per_branch_from = {from_dense, {0, 0}, 1}; + + StateEstimationInput input{}; + input.measured_branch_from_current = { + {.angle_measurement_type = local_angle, .measurement = measurement_a}, + {.angle_measurement_type = local_angle, .measurement = measurement_b}}; + + MeasuredValues const values{std::make_shared(std::move(topo)), input}; + + REQUIRE(values.has_branch_from_current(0)); + CHECK_FALSE(values.has_branch_from_power(0)); + CHECK_FALSE(values.has_branch_to_power(0)); + CHECK_FALSE(values.has_branch_to_current(0)); + + CHECK(values.branch_from_current(0).angle_measurement_type == local_angle); + check_accumulated(values.branch_from_current(0).measurement); + } + + SUBCASE("Accumulate global angle current sensors on branch - sym") { + auto topo = get_single_branch_topology(); + topo.current_sensors_per_branch_from = {from_dense, {0, 0}, 1}; + topo.voltage_sensors_per_bus = {from_dense, {0}, 2}; + + StateEstimationInput input{}; + input.measured_voltage = {{.value = DoubleComplex{1.0, 1.0}, .variance = 10.0}}; + input.measured_branch_from_current = { + {.angle_measurement_type = global_angle, .measurement = measurement_a}, + {.angle_measurement_type = global_angle, .measurement = measurement_b}}; + + MeasuredValues const values{std::make_shared(std::move(topo)), input}; + + REQUIRE(values.has_branch_from_current(0)); + CHECK_FALSE(values.has_branch_from_power(0)); + CHECK_FALSE(values.has_branch_to_power(0)); + CHECK_FALSE(values.has_branch_to_current(0)); + + CHECK(values.branch_from_current(0).angle_measurement_type == global_angle); + check_accumulated(values.branch_from_current(0).measurement); + } + + SUBCASE("Cannot accumulate global angle current sensors if there is no voltage angle measurement") { + auto topo = get_single_branch_topology(); + topo.current_sensors_per_branch_from = {from_dense, {0, 0}, 1}; + + StateEstimationInput input{}; + input.measured_branch_from_current = { + {.angle_measurement_type = global_angle, .measurement = measurement_a}, + {.angle_measurement_type = global_angle, .measurement = measurement_b}}; + + auto const create_measured_values = [&topo, &input] { + return MeasuredValues{std::make_shared(topo), input}; + }; + + SUBCASE("No voltage measurement in topo") { + CHECK_THROWS_AS(create_measured_values(), ConflictingAngleMeasurementType); + } + SUBCASE("No voltage measurement values") { + topo.voltage_sensors_per_bus = {from_dense, {0}, 2}; + input.measured_voltage = {{.value = DoubleComplex{nan, nan}, .variance = nan}}; + CHECK_THROWS_AS(create_measured_values(), ConflictingAngleMeasurementType); + } + SUBCASE("No voltage angle measurement") { + topo.voltage_sensors_per_bus = {from_dense, {0}, 2}; + input.measured_voltage = {{.value = DoubleComplex{1.0, nan}, .variance = 0.1}}; + CHECK_THROWS_AS(create_measured_values(), ConflictingAngleMeasurementType); + } + SUBCASE("No voltage angle measurement") { + topo.voltage_sensors_per_bus = {from_dense, {0}, 2}; + input.measured_voltage = {{.value = DoubleComplex{1.0, 1.0}, .variance = 0.1}}; + CHECK_NOTHROW(create_measured_values()); + } + } + + SUBCASE("Cannot accumulate different angle measurement types on same terminal") { + auto topo = get_single_branch_topology(); + topo.current_sensors_per_branch_from = {from_dense, {0, 0}, 1}; + + StateEstimationInput input{}; + input.measured_branch_from_current = { + {.angle_measurement_type = AngleMeasurementType::local_angle, .measurement = measurement_a}, + {.angle_measurement_type = AngleMeasurementType::global_angle, .measurement = measurement_b}}; + + auto const create_measured_values = [topo_ = std::move(topo), &input] { + return MeasuredValues{std::make_shared(std::move(topo_)), input}; + }; + CHECK_THROWS_AS(create_measured_values(), ConflictingAngleMeasurementType); + } + } } } // namespace power_grid_model::math_solver \ No newline at end of file From 61781e847471590231a811a90aa5c67fa250236e Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Tue, 8 Apr 2025 09:03:10 +0200 Subject: [PATCH 10/25] make gcc happy Signed-off-by: Martijn Govers --- .../include/power_grid_model/common/statistics.hpp | 1 - 1 file changed, 1 deletion(-) 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 239c981a68..f255dca785 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 @@ -222,7 +222,6 @@ constexpr auto accumulate(std::ranges::view auto measurements) using RandVarType = std::ranges::range_value_t; using ValueType = decltype(RandVarType::value); using VarianceType = decltype(RandVarType::variance); - using sym = RandVarType::sym; VarianceType accumulated_inverse_variance{}; ValueType accumulated_value{}; From fb6606467d6d9d6729106703748cc55904ae6f34 Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Tue, 8 Apr 2025 16:02:48 +0200 Subject: [PATCH 11/25] well this is why you need testing Signed-off-by: Martijn Govers --- .../math_solver/observability.hpp | 5 +++-- tests/cpp_unit_tests/test_observability.cpp | 21 ++++++++++++------- 2 files changed, 17 insertions(+), 9 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp index 69c85bcbb6..2994387ed5 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp @@ -63,8 +63,9 @@ std::pair, bool> count_flow_sensors(MeasuredValues cons // we only need one flow sensor, so the loop will break if (element.element_type != YBusElementType::shunt) { Idx const branch = element.idx; - if ((measured_values.has_branch_from_power(branch) || - measured_values.has_branch_to_power(branch)) && + if ((measured_values.has_branch_from_power(branch) || measured_values.has_branch_to_power(branch) || + measured_values.has_branch_from_current(branch) || + measured_values.has_branch_to_current(branch)) && topo.branch_bus_idx[branch][0] != -1 && topo.branch_bus_idx[branch][1] != -1) { flow_sensors[ybus_index] = 1; has_at_least_one_sensor = true; diff --git a/tests/cpp_unit_tests/test_observability.cpp b/tests/cpp_unit_tests/test_observability.cpp index efebf3b3ed..26967d8002 100644 --- a/tests/cpp_unit_tests/test_observability.cpp +++ b/tests/cpp_unit_tests/test_observability.cpp @@ -119,18 +119,25 @@ TEST_CASE("Necessary observability check") { check_not_observable(topo, param, se_input); } SUBCASE("Power sensor and current sensor are equivalent") { + auto const branch_sensor_topo = topo.power_sensors_per_branch_from; + auto const branch_sensor_measurement = se_input.measured_branch_from_power; + check_observable(topo, param, se_input); + // remove the power sensor -> not observable - topo.power_sensors_per_branch_from = {from_sparse, {0, 0, 0, 0}}; + topo.power_sensors_per_branch_from = {from_dense, {}, 3}; se_input.measured_branch_from_power = {}; check_not_observable(topo, param, se_input); // add the current sensor -> observable - topo.current_sensors_per_branch_from = {from_sparse, {0, 1, 1, 1}}; - se_input.measured_branch_from_current = {{.angle_measurement_type = AngleMeasurementType::local_angle, - .measurement = {.real_component = {.value = 1.0, .variance = 1.0}, - .imag_component = {.value = 0.0, .variance = 1.0}}}}; - - check_observable(topo, param, se_input); + topo.current_sensors_per_branch_from = branch_sensor_topo; + for (auto angle_measurement_type : {AngleMeasurementType::local_angle, AngleMeasurementType::global_angle}) { + CAPTURE(angle_measurement_type); + REQUIRE(branch_sensor_measurement.size() == 1); + se_input.measured_branch_from_current = { + {.angle_measurement_type = angle_measurement_type, .measurement = branch_sensor_measurement.front()}}; + + check_observable(topo, param, se_input); + } } // namespace power_grid_model } From 3268e6abe69fc5b0f72c81952f443f46eafe9832 Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Tue, 8 Apr 2025 17:18:39 +0200 Subject: [PATCH 12/25] do not read out of bounds (thank you tests) Signed-off-by: Martijn Govers --- .../calculation_parameters.hpp | 4 ++++ .../power_grid_model/main_model_impl.hpp | 18 ++++++++++++++++++ 2 files changed, 22 insertions(+) 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 5929941044..904b64401a 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 @@ -182,6 +182,10 @@ struct MathModelTopology { Idx n_bus_power_sensor() const { return power_sensors_per_bus.element_size(); } + Idx n_branch_from_current_sensor() const { return current_sensors_per_branch_from.element_size(); } + + Idx n_branch_to_current_sensor() const { return current_sensors_per_branch_to.element_size(); } + Idx n_transformer_tap_regulator() const { return tap_regulators_per_branch.element_size(); } }; diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/main_model_impl.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/main_model_impl.hpp index cd54e8e7b8..4607e7f6ff 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/main_model_impl.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/main_model_impl.hpp @@ -1171,6 +1171,8 @@ class MainModelImpl, ComponentLis se_input[i].measured_branch_from_power.resize(state.math_topology[i]->n_branch_from_power_sensor()); se_input[i].measured_branch_to_power.resize(state.math_topology[i]->n_branch_to_power_sensor()); se_input[i].measured_bus_injection.resize(state.math_topology[i]->n_bus_power_sensor()); + se_input[i].measured_branch_from_current.resize(state.math_topology[i]->n_branch_from_current_sensor()); + se_input[i].measured_branch_to_current.resize(state.math_topology[i]->n_branch_to_current_sensor()); } prepare_input_status::shunt_status, Shunt>(state, state.topo_comp_coup->shunt, @@ -1217,6 +1219,22 @@ class MainModelImpl, ComponentLis state, state.topo_comp_coup->power_sensor, se_input, [&state](Idx i) { return state.comp_topo->power_sensor_terminal_type[i] == MeasuredTerminalType::node; }); + prepare_input, CurrentSensorCalcParam, + &StateEstimationInput::measured_branch_from_current, GenericCurrentSensor>( + state, state.topo_comp_coup->current_sensor, se_input, [&state](Idx i) { + using enum MeasuredTerminalType; + return state.comp_topo->current_sensor_terminal_type[i] == branch_from || + // all branch3 sensors are at from side in the mathematical model + state.comp_topo->current_sensor_terminal_type[i] == branch3_1 || + state.comp_topo->current_sensor_terminal_type[i] == branch3_2 || + state.comp_topo->current_sensor_terminal_type[i] == branch3_3; + }); + prepare_input, CurrentSensorCalcParam, + &StateEstimationInput::measured_branch_to_current, GenericCurrentSensor>( + state, state.topo_comp_coup->current_sensor, se_input, [&state](Idx i) { + return state.comp_topo->current_sensor_terminal_type[i] == MeasuredTerminalType::branch_to; + }); + return se_input; } From 8cc1bf2a11ad20ed83082addf2a1dbf20e43766d Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Tue, 8 Apr 2025 17:23:18 +0200 Subject: [PATCH 13/25] actually throw sparse matrix error Signed-off-by: Martijn Govers --- tests/cpp_unit_tests/test_observability.cpp | 2 +- tests/data/state_estimation/basic-current-sensor/params.json | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/cpp_unit_tests/test_observability.cpp b/tests/cpp_unit_tests/test_observability.cpp index 26967d8002..52f84a5024 100644 --- a/tests/cpp_unit_tests/test_observability.cpp +++ b/tests/cpp_unit_tests/test_observability.cpp @@ -138,7 +138,7 @@ TEST_CASE("Necessary observability check") { check_observable(topo, param, se_input); } - } // namespace power_grid_model + } } } // namespace power_grid_model diff --git a/tests/data/state_estimation/basic-current-sensor/params.json b/tests/data/state_estimation/basic-current-sensor/params.json index b8769a93ef..1646fba785 100644 --- a/tests/data/state_estimation/basic-current-sensor/params.json +++ b/tests/data/state_estimation/basic-current-sensor/params.json @@ -12,14 +12,14 @@ "iterative_linear": { "experimental_features": "enabled", "fail": { - "raises": "NotObservableError", + "raises": "SparseMatrixError", "reason": "Current sensors are not yet implemented for this calculation method" } }, "newton_raphson": { "experimental_features": "enabled", "fail": { - "raises": "NotObservableError", + "raises": "SparseMatrixError", "reason": "Current sensors are not yet implemented for this calculation method" } } From dd9233b40bf1517427f64f0df954e6a90424d61f Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Tue, 8 Apr 2025 17:26:13 +0200 Subject: [PATCH 14/25] fix python test Signed-off-by: Martijn Govers --- src/power_grid_model/enum.py | 13 +++++++++++++ tests/unit/test_error_handling.py | 19 +++++++++++++------ 2 files changed, 26 insertions(+), 6 deletions(-) diff --git a/src/power_grid_model/enum.py b/src/power_grid_model/enum.py index d6ae356007..7eefd773bc 100644 --- a/src/power_grid_model/enum.py +++ b/src/power_grid_model/enum.py @@ -137,6 +137,19 @@ class MeasuredTerminalType(IntEnum): """ +class AngleMeasurementType(IntEnum): + """The type of the angle measurement for current sensors.""" + + local_angle = 0 + """ + The angle is relative to the local voltage angle + """ + global_angle = 1 + """ + The angle is relative to the global voltage angle. + """ + + class FaultType(IntEnum): """The type of fault represented by a fault component""" diff --git a/tests/unit/test_error_handling.py b/tests/unit/test_error_handling.py index 92ed97a322..515fd9be9f 100644 --- a/tests/unit/test_error_handling.py +++ b/tests/unit/test_error_handling.py @@ -9,6 +9,7 @@ from power_grid_model import PowerGridModel from power_grid_model._core.power_grid_meta import initialize_array from power_grid_model.enum import ( + AngleMeasurementType, CalculationMethod, LoadGenType, MeasuredTerminalType, @@ -361,7 +362,7 @@ def test_conflicting_angle_measurement_type(): node_input["id"] = [0, 1] node_input["u_rated"] = [10e3, 10e3] - line_input = initialize_array("input", "transformer", 1) + line_input = initialize_array("input", "line", 1) line_input["id"] = [2] line_input["from_node"] = [0] line_input["to_node"] = [1] @@ -372,16 +373,22 @@ def test_conflicting_angle_measurement_type(): source_input["status"] = [1] source_input["u_ref"] = [1.0] - sym_current_sensor_input = initialize_array("input", "sym_current_sensor", 1) + sym_current_sensor_input = initialize_array("input", "sym_current_sensor", 2) sym_current_sensor_input["id"] = [4, 5] sym_current_sensor_input["measured_object"] = [2, 2] - sym_current_sensor_input["measured_terminal_type"] = [MeasuredTerminalType.branch_from, MeasuredTerminalType.branch_from] - sym_current_sensor_input["angle_measurement_type"] = [AngleMeasurementType.local, AngleMeasurementType.global] - sym_current_sensor_input["status"] = [1, 1] + sym_current_sensor_input["measured_terminal_type"] = [ + MeasuredTerminalType.branch_from, + MeasuredTerminalType.branch_from, + ] + sym_current_sensor_input["angle_measurement_type"] = [ + AngleMeasurementType.local_angle, + AngleMeasurementType.global_angle, + ] sym_current_sensor_input["i_sigma"] = [1.0, 1.0] sym_current_sensor_input["i_angle_sigma"] = [1.0, 1.0] sym_current_sensor_input["i_measured"] = [0.0, 0.0] sym_current_sensor_input["i_angle_measured"] = [0.0, 0.0] + repr(sym_current_sensor_input) model = PowerGridModel( input_data={ @@ -391,7 +398,7 @@ def test_conflicting_angle_measurement_type(): "sym_current_sensor": sym_current_sensor_input, } ) - + with pytest.raises(ConflictingAngleMeasurementType): model._calculate_state_estimation(decode_error=True, experimental_features=_ExperimentalFeatures.enabled) From 2f9dfe4b6618faa2ec3c081a1fff090edf5e4eb3 Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Wed, 9 Apr 2025 08:40:29 +0200 Subject: [PATCH 15/25] fix sonar cloud Signed-off-by: Martijn Govers --- .../power_grid_model/common/statistics.hpp | 41 ++++++++++--------- .../math_solver/measured_values.hpp | 8 ++-- tests/cpp_unit_tests/test_measured_values.cpp | 6 +-- 3 files changed, 28 insertions(+), 27 deletions(-) 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 f255dca785..99bcab7f9e 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 @@ -208,25 +208,25 @@ template struct PolarComplexRandVar { }; namespace statistics { -// combine multiple measurements of one quantity using Kalman filter -constexpr auto accumulate(std::ranges::view auto measurements) - requires std::same_as, - UniformRealRandVar::sym>> || - std::same_as, - IndependentRealRandVar::sym>> || - std::same_as, - UniformComplexRandVar::sym>> || - std::same_as, - IndependentComplexRandVar::sym>> -{ - using RandVarType = std::ranges::range_value_t; +// combine multiple random variables of one quantity using Kalman filter +template + requires std::same_as, + UniformRealRandVar::sym>> || + std::same_as, + IndependentRealRandVar::sym>> || + std::same_as, + UniformComplexRandVar::sym>> || + std::same_as, + IndependentComplexRandVar::sym>> +constexpr auto accumulate(RandVarsView rand_vars) { + using RandVarType = std::ranges::range_value_t; using ValueType = decltype(RandVarType::value); using VarianceType = decltype(RandVarType::variance); VarianceType accumulated_inverse_variance{}; ValueType accumulated_value{}; - std::ranges::for_each(measurements, [&accumulated_inverse_variance, &accumulated_value](auto const& measurement) { + std::ranges::for_each(rand_vars, [&accumulated_inverse_variance, &accumulated_value](auto const& measurement) { accumulated_inverse_variance += VarianceType{1.0} / measurement.variance; accumulated_value += measurement.value / measurement.variance; }); @@ -237,17 +237,18 @@ constexpr auto accumulate(std::ranges::view auto measurements) } return RandVarType{.value = accumulated_value, .variance = VarianceType{std::numeric_limits::infinity()}}; } -constexpr auto accumulate(std::ranges::view auto measurements) - requires std::same_as, - DecomposedComplexRandVar::sym>> -{ - using sym = std::ranges::range_value_t::sym; + +template + requires std::same_as, + DecomposedComplexRandVar::sym>> +constexpr auto accumulate(RandVarsView rand_vars) { + using sym = std::ranges::range_value_t::sym; DecomposedComplexRandVar result{ .real_component = - accumulate(measurements | std::views::transform([](auto const& x) -> auto& { return x.real_component; })), + accumulate(rand_vars | std::views::transform([](auto const& x) -> auto& { return x.real_component; })), .imag_component = - accumulate(measurements | std::views::transform([](auto const& x) -> auto& { return x.imag_component; }))}; + accumulate(rand_vars | std::views::transform([](auto const& x) -> auto& { return x.imag_component; }))}; if (!(is_normal(result.real_component.variance) && is_normal(result.imag_component.variance))) { result.real_component.variance = RealValue{std::numeric_limits::infinity()}; diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp index da09a1ab27..b93764486a 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp @@ -453,7 +453,7 @@ template class MeasuredValues { template static VoltageSensorCalcParam combine_measurements(std::vector> const& data, IdxRange const& sensors) { - auto complex_measurements = sensors | std::views::transform([&](Idx pos) -> auto& { return data[pos]; }); + auto complex_measurements = sensors | std::views::transform([&data](Idx pos) -> auto& { return data[pos]; }); if constexpr (only_magnitude) { auto const combined_magnitude_measurement = statistics::accumulate( complex_measurements | std::views::transform([](auto const& measurement) { @@ -476,13 +476,13 @@ template class MeasuredValues { requires(!only_magnitude) static PowerSensorCalcParam combine_measurements(std::vector> const& data, IdxRange const& sensors) { - return statistics::accumulate(sensors | std::views::transform([&](Idx pos) -> auto& { return data[pos]; })); + return statistics::accumulate(sensors | std::views::transform([&data](Idx pos) -> auto& { return data[pos]; })); } template requires(!only_magnitude) static CurrentSensorCalcParam combine_measurements(std::vector> const& data, IdxRange const& sensors) { - auto const params = sensors | std::views::transform([&](Idx pos) -> auto& { return data[pos]; }); + auto const params = sensors | std::views::transform([&data](Idx pos) -> auto& { return data[pos]; }); auto const angle_measurement_type = sensors.empty() ? AngleMeasurementType::local_angle // fallback : params.front().angle_measurement_type; if (std::ranges::any_of(params, [angle_measurement_type](auto const& params) { @@ -494,7 +494,7 @@ template class MeasuredValues { return {.angle_measurement_type = angle_measurement_type, .measurement = statistics::accumulate( - params | std::views::transform([&](auto const& params) -> auto& { return params.measurement; }))}; + params | std::views::transform([](auto const& param) -> auto& { return param.measurement; }))}; } template diff --git a/tests/cpp_unit_tests/test_measured_values.cpp b/tests/cpp_unit_tests/test_measured_values.cpp index e265058c58..6e9aced541 100644 --- a/tests/cpp_unit_tests/test_measured_values.cpp +++ b/tests/cpp_unit_tests/test_measured_values.cpp @@ -261,12 +261,12 @@ TEST_CASE("Measured Values") { {.angle_measurement_type = AngleMeasurementType::local_angle, .measurement = measurement_a}, {.angle_measurement_type = AngleMeasurementType::global_angle, .measurement = measurement_b}}; - auto const create_measured_values = [topo_ = std::move(topo), &input] { - return MeasuredValues{std::make_shared(std::move(topo_)), input}; + auto const create_measured_values = [&topo, &input] { + return MeasuredValues{std::make_shared(topo), input}; }; CHECK_THROWS_AS(create_measured_values(), ConflictingAngleMeasurementType); } } } -} // namespace power_grid_model::math_solver \ No newline at end of file +} // namespace power_grid_model::math_solver From 888271aa2af6de157e4a3b44c1c0fbb3c551db9a Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Thu, 10 Apr 2025 08:01:54 +0200 Subject: [PATCH 16/25] make missing voltage angle not observable Signed-off-by: Martijn Govers --- .../math_solver/measured_values.hpp | 10 ++-- .../math_solver/observability.hpp | 5 ++ tests/unit/test_error_handling.py | 60 ++++++++++++++++++- 3 files changed, 68 insertions(+), 7 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp index b93764486a..18e4a69cd6 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp @@ -71,6 +71,7 @@ template class MeasuredValues { constexpr bool has_angle() const { return n_voltage_angle_measurements_ > 0; } constexpr bool has_voltage_measurements() const { return n_voltage_measurements_ > 0; } + constexpr bool has_global_angle_current() const { return n_global_angle_current_measurements_ >= 0; } constexpr bool has_voltage(Idx bus) const { return idx_voltage_[bus] >= 0; } constexpr bool has_angle_measurement(Idx bus) const { return !is_nan(imag(voltage(bus))); } @@ -215,6 +216,7 @@ template class MeasuredValues { Idx n_voltage_measurements_{}; Idx n_voltage_angle_measurements_{}; + Idx n_global_angle_current_measurements_{}; // average angle shift of voltages with angle measurement // default is zero is no voltage has angle measurement @@ -437,12 +439,10 @@ template class MeasuredValues { process_one_object(branch, topo.current_sensors_per_branch_to, topo.branch_bus_idx, input.measured_branch_to_current, current_main_value_, branch_to_checker); - if ((!has_angle()) && std::ranges::any_of(current_main_value_, [](auto const& measurement) { + n_global_angle_current_measurements_ = + std::ranges::count_if(current_main_value_, [](auto const& measurement) { return measurement.angle_measurement_type == AngleMeasurementType::global_angle; - })) { - throw ConflictingAngleMeasurementType{ - "Global angle current measurements require a voltage angle measurement in the connected subgrid."}; - } + }); } } diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp index 2994387ed5..13d54fca8a 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp @@ -169,6 +169,11 @@ inline ObservabilityResult necessary_observability_check(MeasuredValues con result.is_sufficiently_observable = true; } + if (measured_values.has_global_angle_current() && n_voltage_phasor_sensor == 0) { + throw NotObservableError{ + "Global angle current sensors require at least one voltage angle measurement as a reference point.\n"}; + } + return result; } diff --git a/tests/unit/test_error_handling.py b/tests/unit/test_error_handling.py index 515fd9be9f..5a439b10b9 100644 --- a/tests/unit/test_error_handling.py +++ b/tests/unit/test_error_handling.py @@ -357,7 +357,7 @@ def test_automatic_tap_changing(): model.calculate_power_flow(tap_changing_strategy=TapChangingStrategy.min_voltage_tap) -def test_conflicting_angle_measurement_type(): +def test_conflicting_angle_measurement_type() -> None: node_input = initialize_array("input", "node", 2) node_input["id"] = [0, 1] node_input["u_rated"] = [10e3, 10e3] @@ -373,6 +373,14 @@ def test_conflicting_angle_measurement_type(): source_input["status"] = [1] source_input["u_ref"] = [1.0] + sym_voltage_sensor_input = initialize_array("input", "sym_voltage_sensor", 1) + sym_voltage_sensor_input["id"] = [6] + sym_voltage_sensor_input["measured_object"] = [0] + sym_voltage_sensor_input["u_sigma"] = [1.0] + sym_voltage_sensor_input["u_angle_sigma"] = [1.0] + sym_voltage_sensor_input["u_measured"] = [1.0] + sym_voltage_sensor_input["u_angle_measured"] = [0.0] + sym_current_sensor_input = initialize_array("input", "sym_current_sensor", 2) sym_current_sensor_input["id"] = [4, 5] sym_current_sensor_input["measured_object"] = [2, 2] @@ -388,13 +396,13 @@ def test_conflicting_angle_measurement_type(): sym_current_sensor_input["i_angle_sigma"] = [1.0, 1.0] sym_current_sensor_input["i_measured"] = [0.0, 0.0] sym_current_sensor_input["i_angle_measured"] = [0.0, 0.0] - repr(sym_current_sensor_input) model = PowerGridModel( input_data={ "node": node_input, "line": line_input, "source": source_input, + "sym_voltage_sensor": sym_voltage_sensor_input, "sym_current_sensor": sym_current_sensor_input, } ) @@ -403,6 +411,54 @@ def test_conflicting_angle_measurement_type(): model._calculate_state_estimation(decode_error=True, experimental_features=_ExperimentalFeatures.enabled) +def test_global_current_measurement_without_voltage_angle() -> None: + node_input = initialize_array("input", "node", 2) + node_input["id"] = [0, 1] + node_input["u_rated"] = [10e3, 10e3] + + line_input = initialize_array("input", "line", 1) + line_input["id"] = [2] + line_input["from_node"] = [0] + line_input["to_node"] = [1] + + source_input = initialize_array("input", "source", 1) + source_input["id"] = [3] + source_input["node"] = [0] + source_input["status"] = [1] + source_input["u_ref"] = [1.0] + + sym_voltage_sensor_input = initialize_array("input", "sym_voltage_sensor", 1) + sym_voltage_sensor_input["id"] = [4] + sym_voltage_sensor_input["measured_object"] = [0] + sym_voltage_sensor_input["u_sigma"] = [1.0] + sym_voltage_sensor_input["u_measured"] = [0.0] + + sym_current_sensor_input = initialize_array("input", "sym_current_sensor", 2) + sym_current_sensor_input["id"] = [5] + sym_current_sensor_input["measured_object"] = [2] + sym_current_sensor_input["measured_terminal_type"] = [ + MeasuredTerminalType.branch_from, + ] + sym_current_sensor_input["angle_measurement_type"] = [ + AngleMeasurementType.global_angle, + ] + sym_current_sensor_input["i_sigma"] = [1.0, 1.0] + sym_current_sensor_input["i_measured"] = [0.0, 0.0] + + model = PowerGridModel( + input_data={ + "node": node_input, + "line": line_input, + "source": source_input, + "sym_voltage_sensor": sym_voltage_sensor_input, + "sym_current_sensor": sym_current_sensor_input, + } + ) + + with pytest.raises(NotObservableError): + model._calculate_state_estimation(decode_error=True, experimental_features=_ExperimentalFeatures.enabled) + + @pytest.mark.skip(reason="TODO") def test_handle_power_grid_dataset_error(): pass From 0aae93c4b278ebebccb6569e084f9f913001410e Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Thu, 10 Apr 2025 08:29:18 +0200 Subject: [PATCH 17/25] move global angle current sensor with voltage angle to observability Signed-off-by: Martijn Govers --- .../math_solver/measured_values.hpp | 2 +- .../math_solver/observability.hpp | 15 +++-- tests/cpp_unit_tests/test_measured_values.cpp | 33 ---------- tests/cpp_unit_tests/test_observability.cpp | 64 ++++++++++++++----- 4 files changed, 56 insertions(+), 58 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp index 18e4a69cd6..d0a7ad67f7 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp @@ -71,7 +71,7 @@ template class MeasuredValues { constexpr bool has_angle() const { return n_voltage_angle_measurements_ > 0; } constexpr bool has_voltage_measurements() const { return n_voltage_measurements_ > 0; } - constexpr bool has_global_angle_current() const { return n_global_angle_current_measurements_ >= 0; } + constexpr bool has_global_angle_current() const { return n_global_angle_current_measurements_ > 0; } constexpr bool has_voltage(Idx bus) const { return idx_voltage_[bus] >= 0; } constexpr bool has_angle_measurement(Idx bus) const { return !is_nan(imag(voltage(bus))); } diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp index 13d54fca8a..1c9e1493ef 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp @@ -155,6 +155,11 @@ inline ObservabilityResult necessary_observability_check(MeasuredValues con throw NotObservableError{}; } + if (measured_values.has_global_angle_current() && n_voltage_phasor_sensor == 0) { + throw NotObservableError{ + "Global angle current sensors require at least one voltage angle measurement as a reference point.\n"}; + } + // for radial grid without phasor measurement, try to assign injection sensor to branch sensor // we can then check sufficient condition for observability if (topo.is_radial && n_voltage_phasor_sensor == 0) { @@ -163,17 +168,13 @@ inline ObservabilityResult necessary_observability_check(MeasuredValues con if (Idx const n_flow_sensor_new = std::reduce(flow_sensors.cbegin(), flow_sensors.cend(), Idx{}, std::plus{}); n_flow_sensor_new < n_bus - 1) { - throw NotObservableError{"The number of power sensors appears sufficient, but they are not independent " - "enough. The system is still not observable.\n"}; + throw NotObservableError{ + "The number of power/current sensors appears sufficient, but they are not independent " + "enough. The system is still not observable.\n"}; } result.is_sufficiently_observable = true; } - if (measured_values.has_global_angle_current() && n_voltage_phasor_sensor == 0) { - throw NotObservableError{ - "Global angle current sensors require at least one voltage angle measurement as a reference point.\n"}; - } - return result; } diff --git a/tests/cpp_unit_tests/test_measured_values.cpp b/tests/cpp_unit_tests/test_measured_values.cpp index 6e9aced541..83ffea6e5d 100644 --- a/tests/cpp_unit_tests/test_measured_values.cpp +++ b/tests/cpp_unit_tests/test_measured_values.cpp @@ -219,39 +219,6 @@ TEST_CASE("Measured Values") { check_accumulated(values.branch_from_current(0).measurement); } - SUBCASE("Cannot accumulate global angle current sensors if there is no voltage angle measurement") { - auto topo = get_single_branch_topology(); - topo.current_sensors_per_branch_from = {from_dense, {0, 0}, 1}; - - StateEstimationInput input{}; - input.measured_branch_from_current = { - {.angle_measurement_type = global_angle, .measurement = measurement_a}, - {.angle_measurement_type = global_angle, .measurement = measurement_b}}; - - auto const create_measured_values = [&topo, &input] { - return MeasuredValues{std::make_shared(topo), input}; - }; - - SUBCASE("No voltage measurement in topo") { - CHECK_THROWS_AS(create_measured_values(), ConflictingAngleMeasurementType); - } - SUBCASE("No voltage measurement values") { - topo.voltage_sensors_per_bus = {from_dense, {0}, 2}; - input.measured_voltage = {{.value = DoubleComplex{nan, nan}, .variance = nan}}; - CHECK_THROWS_AS(create_measured_values(), ConflictingAngleMeasurementType); - } - SUBCASE("No voltage angle measurement") { - topo.voltage_sensors_per_bus = {from_dense, {0}, 2}; - input.measured_voltage = {{.value = DoubleComplex{1.0, nan}, .variance = 0.1}}; - CHECK_THROWS_AS(create_measured_values(), ConflictingAngleMeasurementType); - } - SUBCASE("No voltage angle measurement") { - topo.voltage_sensors_per_bus = {from_dense, {0}, 2}; - input.measured_voltage = {{.value = DoubleComplex{1.0, 1.0}, .variance = 0.1}}; - CHECK_NOTHROW(create_measured_values()); - } - } - SUBCASE("Cannot accumulate different angle measurement types on same terminal") { auto topo = get_single_branch_topology(); topo.current_sensors_per_branch_from = {from_dense, {0, 0}, 1}; diff --git a/tests/cpp_unit_tests/test_observability.cpp b/tests/cpp_unit_tests/test_observability.cpp index 52f84a5024..4b099a0a8f 100644 --- a/tests/cpp_unit_tests/test_observability.cpp +++ b/tests/cpp_unit_tests/test_observability.cpp @@ -118,27 +118,57 @@ TEST_CASE("Necessary observability check") { // this will throw NotObservableError check_not_observable(topo, param, se_input); } - SUBCASE("Power sensor and current sensor are equivalent") { - auto const branch_sensor_topo = topo.power_sensors_per_branch_from; - auto const branch_sensor_measurement = se_input.measured_branch_from_power; - check_observable(topo, param, se_input); - - // remove the power sensor -> not observable + SUBCASE("Current sensors also measure branch flow") { topo.power_sensors_per_branch_from = {from_dense, {}, 3}; se_input.measured_branch_from_power = {}; - check_not_observable(topo, param, se_input); - - // add the current sensor -> observable - topo.current_sensors_per_branch_from = branch_sensor_topo; - for (auto angle_measurement_type : {AngleMeasurementType::local_angle, AngleMeasurementType::global_angle}) { - CAPTURE(angle_measurement_type); - REQUIRE(branch_sensor_measurement.size() == 1); - se_input.measured_branch_from_current = { - {.angle_measurement_type = angle_measurement_type, .measurement = branch_sensor_measurement.front()}}; + topo.current_sensors_per_branch_from = {from_dense, {2}, 3}; + + DecomposedComplexRandVar const current_measurement{ + .real_component = {.value = 1.0, .variance = 1.0}, .imag_component = {.value = 0.0, .variance = 1.0}}; + + SUBCASE("With voltage phasor measurement") { + SUBCASE("Local current sensor") { + se_input.measured_branch_from_current = { + {.angle_measurement_type = AngleMeasurementType::local_angle, .measurement = current_measurement}}; + check_observable(topo, param, se_input); + } + SUBCASE("Global angle current sensor") { + se_input.measured_branch_from_current = { + {.angle_measurement_type = AngleMeasurementType::global_angle, .measurement = current_measurement}}; + check_observable(topo, param, se_input); + } + } + SUBCASE("No voltage phasor measurement and single current sensor") { + se_input.measured_voltage = {{.value = {1.0, nan}, .variance = 1.0}}; - check_observable(topo, param, se_input); + SUBCASE("Local current sensor") { + se_input.measured_branch_from_current = { + {.angle_measurement_type = AngleMeasurementType::local_angle, .measurement = current_measurement}}; + check_not_observable(topo, param, se_input); + } + SUBCASE("Global angle current sensor") { + se_input.measured_branch_from_current = { + {.angle_measurement_type = AngleMeasurementType::global_angle, .measurement = current_measurement}}; + check_not_observable(topo, param, se_input); + } + } + SUBCASE("No voltage phasor measurement and two current sensors") { + se_input.measured_voltage = {{.value = {1.0, nan}, .variance = 1.0}}; + topo.current_sensors_per_branch_from = {from_dense, {0, 2}, 3}; + + SUBCASE("Local current sensor") { + se_input.measured_branch_from_current = { + {.angle_measurement_type = AngleMeasurementType::local_angle, .measurement = current_measurement}, + {.angle_measurement_type = AngleMeasurementType::local_angle, .measurement = current_measurement}}; + check_observable(topo, param, se_input); + } + SUBCASE("Global angle current sensor") { + se_input.measured_branch_from_current = { + {.angle_measurement_type = AngleMeasurementType::global_angle, .measurement = current_measurement}, + {.angle_measurement_type = AngleMeasurementType::global_angle, .measurement = current_measurement}}; + check_not_observable(topo, param, se_input); + } } } } - } // namespace power_grid_model From d65427d11362562aedba0da9f1e5331ad5b97263 Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Thu, 10 Apr 2025 10:08:13 +0200 Subject: [PATCH 18/25] fix test Signed-off-by: Martijn Govers --- tests/unit/test_error_handling.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/tests/unit/test_error_handling.py b/tests/unit/test_error_handling.py index 5a439b10b9..7e53009020 100644 --- a/tests/unit/test_error_handling.py +++ b/tests/unit/test_error_handling.py @@ -374,15 +374,14 @@ def test_conflicting_angle_measurement_type() -> None: source_input["u_ref"] = [1.0] sym_voltage_sensor_input = initialize_array("input", "sym_voltage_sensor", 1) - sym_voltage_sensor_input["id"] = [6] + sym_voltage_sensor_input["id"] = [4] sym_voltage_sensor_input["measured_object"] = [0] sym_voltage_sensor_input["u_sigma"] = [1.0] - sym_voltage_sensor_input["u_angle_sigma"] = [1.0] sym_voltage_sensor_input["u_measured"] = [1.0] sym_voltage_sensor_input["u_angle_measured"] = [0.0] sym_current_sensor_input = initialize_array("input", "sym_current_sensor", 2) - sym_current_sensor_input["id"] = [4, 5] + sym_current_sensor_input["id"] = [5, 6] sym_current_sensor_input["measured_object"] = [2, 2] sym_current_sensor_input["measured_terminal_type"] = [ MeasuredTerminalType.branch_from, @@ -433,7 +432,7 @@ def test_global_current_measurement_without_voltage_angle() -> None: sym_voltage_sensor_input["u_sigma"] = [1.0] sym_voltage_sensor_input["u_measured"] = [0.0] - sym_current_sensor_input = initialize_array("input", "sym_current_sensor", 2) + sym_current_sensor_input = initialize_array("input", "sym_current_sensor", 1) sym_current_sensor_input["id"] = [5] sym_current_sensor_input["measured_object"] = [2] sym_current_sensor_input["measured_terminal_type"] = [ @@ -442,8 +441,8 @@ def test_global_current_measurement_without_voltage_angle() -> None: sym_current_sensor_input["angle_measurement_type"] = [ AngleMeasurementType.global_angle, ] - sym_current_sensor_input["i_sigma"] = [1.0, 1.0] - sym_current_sensor_input["i_measured"] = [0.0, 0.0] + sym_current_sensor_input["i_sigma"] = [1.0] + sym_current_sensor_input["i_measured"] = [0.0] model = PowerGridModel( input_data={ From 2f3dae0701641b4fa76b03987cba219d836ea3f6 Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Thu, 10 Apr 2025 12:52:27 +0200 Subject: [PATCH 19/25] sonar cloud Signed-off-by: Martijn Govers --- .../math_solver/measured_values.hpp | 4 ++-- tests/cpp_unit_tests/test_observability.cpp | 18 ++++++++++-------- 2 files changed, 12 insertions(+), 10 deletions(-) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp index d0a7ad67f7..10721d9fa9 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp @@ -485,8 +485,8 @@ template class MeasuredValues { auto const params = sensors | std::views::transform([&data](Idx pos) -> auto& { return data[pos]; }); auto const angle_measurement_type = sensors.empty() ? AngleMeasurementType::local_angle // fallback : params.front().angle_measurement_type; - if (std::ranges::any_of(params, [angle_measurement_type](auto const& params) { - return params.angle_measurement_type != angle_measurement_type; + if (std::ranges::any_of(params, [angle_measurement_type](auto const& param) { + return param.angle_measurement_type != angle_measurement_type; })) { throw ConflictingAngleMeasurementType{ "Cannot mix local and global angle current measurements on the same terminal."}; diff --git a/tests/cpp_unit_tests/test_observability.cpp b/tests/cpp_unit_tests/test_observability.cpp index 4b099a0a8f..9bc377a64e 100644 --- a/tests/cpp_unit_tests/test_observability.cpp +++ b/tests/cpp_unit_tests/test_observability.cpp @@ -119,6 +119,8 @@ TEST_CASE("Necessary observability check") { check_not_observable(topo, param, se_input); } SUBCASE("Current sensors also measure branch flow") { + using enum AngleMeasurementType; + topo.power_sensors_per_branch_from = {from_dense, {}, 3}; se_input.measured_branch_from_power = {}; topo.current_sensors_per_branch_from = {from_dense, {2}, 3}; @@ -129,12 +131,12 @@ TEST_CASE("Necessary observability check") { SUBCASE("With voltage phasor measurement") { SUBCASE("Local current sensor") { se_input.measured_branch_from_current = { - {.angle_measurement_type = AngleMeasurementType::local_angle, .measurement = current_measurement}}; + {.angle_measurement_type = local_angle, .measurement = current_measurement}}; check_observable(topo, param, se_input); } SUBCASE("Global angle current sensor") { se_input.measured_branch_from_current = { - {.angle_measurement_type = AngleMeasurementType::global_angle, .measurement = current_measurement}}; + {.angle_measurement_type = global_angle, .measurement = current_measurement}}; check_observable(topo, param, se_input); } } @@ -143,12 +145,12 @@ TEST_CASE("Necessary observability check") { SUBCASE("Local current sensor") { se_input.measured_branch_from_current = { - {.angle_measurement_type = AngleMeasurementType::local_angle, .measurement = current_measurement}}; + {.angle_measurement_type = local_angle, .measurement = current_measurement}}; check_not_observable(topo, param, se_input); } SUBCASE("Global angle current sensor") { se_input.measured_branch_from_current = { - {.angle_measurement_type = AngleMeasurementType::global_angle, .measurement = current_measurement}}; + {.angle_measurement_type = global_angle, .measurement = current_measurement}}; check_not_observable(topo, param, se_input); } } @@ -158,14 +160,14 @@ TEST_CASE("Necessary observability check") { SUBCASE("Local current sensor") { se_input.measured_branch_from_current = { - {.angle_measurement_type = AngleMeasurementType::local_angle, .measurement = current_measurement}, - {.angle_measurement_type = AngleMeasurementType::local_angle, .measurement = current_measurement}}; + {.angle_measurement_type = local_angle, .measurement = current_measurement}, + {.angle_measurement_type = local_angle, .measurement = current_measurement}}; check_observable(topo, param, se_input); } SUBCASE("Global angle current sensor") { se_input.measured_branch_from_current = { - {.angle_measurement_type = AngleMeasurementType::global_angle, .measurement = current_measurement}, - {.angle_measurement_type = AngleMeasurementType::global_angle, .measurement = current_measurement}}; + {.angle_measurement_type = global_angle, .measurement = current_measurement}, + {.angle_measurement_type = global_angle, .measurement = current_measurement}}; check_not_observable(topo, param, se_input); } } From 7e5fc8e28a86e13d5373aaeb35dbb31878c9fad1 Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Thu, 10 Apr 2025 16:44:35 +0200 Subject: [PATCH 20/25] resolve some comments Signed-off-by: Martijn Govers --- .../power_grid_model/common/statistics.hpp | 20 ++++++++++--------- .../math_solver/measured_values.hpp | 16 +++++++-------- 2 files changed, 19 insertions(+), 17 deletions(-) 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 99bcab7f9e..7b395635ed 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 @@ -224,18 +224,20 @@ constexpr auto accumulate(RandVarsView rand_vars) { using VarianceType = decltype(RandVarType::variance); VarianceType accumulated_inverse_variance{}; - ValueType accumulated_value{}; + ValueType weighted_accumulated_value{}; - std::ranges::for_each(rand_vars, [&accumulated_inverse_variance, &accumulated_value](auto const& measurement) { - accumulated_inverse_variance += VarianceType{1.0} / measurement.variance; - accumulated_value += measurement.value / measurement.variance; - }); + std::ranges::for_each(rand_vars, + [&accumulated_inverse_variance, &weighted_accumulated_value](auto const& measurement) { + accumulated_inverse_variance += VarianceType{1.0} / measurement.variance; + weighted_accumulated_value += measurement.value / measurement.variance; + }); - if (is_normal(accumulated_inverse_variance)) { - return RandVarType{.value = accumulated_value / accumulated_inverse_variance, - .variance = VarianceType{1.0} / accumulated_inverse_variance}; + if (!is_normal(accumulated_inverse_variance)) { + return RandVarType{.value = weighted_accumulated_value, + .variance = VarianceType{std::numeric_limits::infinity()}}; } - return RandVarType{.value = accumulated_value, .variance = VarianceType{std::numeric_limits::infinity()}}; + return RandVarType{.value = weighted_accumulated_value / accumulated_inverse_variance, + .variance = VarianceType{1.0} / accumulated_inverse_variance}; } template diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp index 10721d9fa9..b45fb469f1 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp @@ -408,15 +408,15 @@ template class MeasuredValues { connected at that side. For each branch the checker checks if the from and to side are connected by checking if branch_bus_idx = disconnected. - If the branch_bus_idx = disconnected, idx_branch_(to|from)_(power|current)_ is set to disconnected. + If the branch_bus_idx = disconnected, idx_branch_(to/from)_(power/current)_ is set to disconnected. If the side is connected, but there are no measurements in this branch side - idx_branch_(to|from)_(power|current)_ is set to disconnected. - Else, idx_branch_(to|from)_(power|current)_ is set to the index of the aggregated data in + idx_branch_(to/from)_(power/current)_ is set to disconnected. + Else, idx_branch_(to/from)_(power/current)_ is set to the index of the aggregated data in power/current_main_value_. All measurement values for a single side of a branch are combined in a weighted average, which is appended to power/current_main_value_. The values in power/current_main_value_ can be found using - idx_branch_(to|from)_(power|current)_. + idx_branch_(to/from)_(power/current)_. */ MathModelTopology const& topo = math_topology(); static constexpr auto branch_from_checker = [](BranchIdx x) { return x[0] != -1; }; @@ -455,19 +455,19 @@ template class MeasuredValues { IdxRange const& sensors) { auto complex_measurements = sensors | std::views::transform([&data](Idx pos) -> auto& { return data[pos]; }); if constexpr (only_magnitude) { - auto const combined_magnitude_measurement = statistics::accumulate( + auto const weighted_average_magnitude_measurement = statistics::accumulate( complex_measurements | std::views::transform([](auto const& measurement) { return UniformRealRandVar{.value = detail::cabs_or_real(measurement.value), .variance = measurement.variance}; })); return UniformComplexRandVar{.value = - [&combined_magnitude_measurement]() { + [&weighted_average_magnitude_measurement]() { ComplexValue abs_value = piecewise_complex_value(DoubleComplex{0.0, nan}); - abs_value += combined_magnitude_measurement.value; + abs_value += weighted_average_magnitude_measurement.value; return abs_value; }(), - .variance = combined_magnitude_measurement.variance}; + .variance = weighted_average_magnitude_measurement.variance}; } else { return statistics::accumulate(complex_measurements); } From f6a400dac8c9d53ced3af922bd9ee753d77dfee3 Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Fri, 11 Apr 2025 08:41:37 +0200 Subject: [PATCH 21/25] add asym tests for measured values Signed-off-by: Martijn Govers --- tests/cpp_unit_tests/test_measured_values.cpp | 211 +++++++++++------- 1 file changed, 134 insertions(+), 77 deletions(-) diff --git a/tests/cpp_unit_tests/test_measured_values.cpp b/tests/cpp_unit_tests/test_measured_values.cpp index 83ffea6e5d..0c8183dcac 100644 --- a/tests/cpp_unit_tests/test_measured_values.cpp +++ b/tests/cpp_unit_tests/test_measured_values.cpp @@ -137,102 +137,159 @@ TEST_CASE("Measured Values") { check_close(injection.real_component.variance, RealValue{0.3, 0.2, 0.85}); check_close(injection.imag_component.variance, RealValue{0.7, 0.85, 0.2}); } +} - SUBCASE("Accumulate branch flow sensors") { - using enum AngleMeasurementType; - - auto const get_single_branch_topology = [] { - auto topo = MathModelTopology{}; - topo.phase_shift = {0.0, 0.0}; - topo.branch_bus_idx = {{0, 1}}; - topo.shunts_per_bus = {from_dense, {}, 1}; - topo.load_gens_per_bus = {from_dense, {}, 1}; - topo.sources_per_bus = {from_dense, {}, 1}; - return topo; - }; - DecomposedComplexRandVar const measurement_a{.real_component = {.value = 1.0, .variance = 1.0}, - .imag_component = {.value = 1.5, .variance = 4.0}}; - DecomposedComplexRandVar const measurement_b{.real_component = {.value = 4.0, .variance = 2.0}, - .imag_component = {.value = 0.7, .variance = 3.0}}; - auto const check_accumulated = [](DecomposedComplexRandVar const& value) { - check_close(value.value(), 2.0 + (73.0i / 70.0)); - check_close(value.real_component.variance, 7.0 / 25.0); - check_close(value.imag_component.variance, 18.0 / 25.0); - }; +TEST_CASE_TEMPLATE("Measured Values - Accumulate branch flow sensors", sym, symmetric_t, asymmetric_t) { + using enum AngleMeasurementType; - SUBCASE("Accumulate power sensors on branch - sym") { - auto topo = get_single_branch_topology(); - topo.power_sensors_per_branch_from = {from_dense, {0, 0}, 1}; + auto const get_single_branch_topology = [] { + auto topo = MathModelTopology{}; + topo.phase_shift = {0.0, 0.0}; + topo.branch_bus_idx = {{0, 1}}; + topo.shunts_per_bus = {from_dense, {}, 1}; + topo.load_gens_per_bus = {from_dense, {}, 1}; + topo.sources_per_bus = {from_dense, {}, 1}; + return topo; + }; + auto const measurement_a = [] { + if constexpr (is_symmetric_v) { + return DecomposedComplexRandVar{.real_component = {.value = 1.0, .variance = 1.0}, + .imag_component = {.value = 1.5, .variance = 4.0}}; + } else { + return DecomposedComplexRandVar{ + .real_component = {.value = RealValue{1.0, 1.0, 1.0}, + .variance = RealValue{1.0, 1.0, 1.0}}, + .imag_component = {.value = RealValue{1.5, 1.5, 1.5}, + .variance = RealValue{4.0, 4.0, 4.0}}}; + } + }(); + auto const measurement_b = [] { + if constexpr (is_symmetric_v) { + return DecomposedComplexRandVar{.real_component = {.value = 4.0, .variance = 2.0}, + .imag_component = {.value = 0.7, .variance = 3.0}}; + } else { + return DecomposedComplexRandVar{ + .real_component = {.value = RealValue{4.0, 4.0, 4.0}, + .variance = RealValue{2.0, 2.0, 2.0}}, + .imag_component = {.value = RealValue{0.7, 0.7, 0.7}, + .variance = RealValue{3.0, 3.0, 3.0}}}; + } + }(); + + auto const check_accumulated = [](DecomposedComplexRandVar const& value) { + auto const expected_value = 2.0 + (73.0i / 70.0); + auto const expected_real_variance = 7.0 / 25.0; + auto const expected_imag_variance = 18.0 / 25.0; + if constexpr (is_symmetric_v) { + check_close(value.value(), expected_value); + check_close(value.real_component.variance, expected_real_variance); + check_close(value.imag_component.variance, expected_imag_variance); + } else { + for (Idx phase : IdxRange(3)) { + check_close(value.value()(phase), expected_value); + check_close(value.real_component.variance(phase), expected_real_variance); + check_close(value.imag_component.variance(phase), expected_imag_variance); + } + } + }; - StateEstimationInput input{}; - input.measured_branch_from_power = {measurement_a, measurement_b}; + SUBCASE("Accumulate power sensors on branch") { + auto topo = get_single_branch_topology(); + topo.power_sensors_per_branch_from = {from_dense, {0, 0}, 1}; - MeasuredValues const values{std::make_shared(std::move(topo)), input}; + StateEstimationInput input{}; + input.measured_branch_from_power = {measurement_a, measurement_b}; - REQUIRE(values.has_branch_from_power(0)); - CHECK_FALSE(values.has_branch_from_current(0)); - CHECK_FALSE(values.has_branch_to_power(0)); - CHECK_FALSE(values.has_branch_to_current(0)); + MeasuredValues const values{std::make_shared(std::move(topo)), input}; - check_accumulated(values.branch_from_power(0)); - } + REQUIRE(values.has_branch_from_power(0)); + CHECK_FALSE(values.has_branch_from_current(0)); + CHECK_FALSE(values.has_branch_to_power(0)); + CHECK_FALSE(values.has_branch_to_current(0)); - SUBCASE("Accumulate local angle current sensors on branch - sym") { - auto topo = get_single_branch_topology(); - topo.current_sensors_per_branch_from = {from_dense, {0, 0}, 1}; + check_accumulated(values.branch_from_power(0)); + } - StateEstimationInput input{}; - input.measured_branch_from_current = { - {.angle_measurement_type = local_angle, .measurement = measurement_a}, - {.angle_measurement_type = local_angle, .measurement = measurement_b}}; + SUBCASE("Accumulate local angle current sensors on branch") { + auto topo = get_single_branch_topology(); + topo.current_sensors_per_branch_from = {from_dense, {0, 0}, 1}; - MeasuredValues const values{std::make_shared(std::move(topo)), input}; + StateEstimationInput input{}; + input.measured_branch_from_current = {{.angle_measurement_type = local_angle, .measurement = measurement_a}, + {.angle_measurement_type = local_angle, .measurement = measurement_b}}; - REQUIRE(values.has_branch_from_current(0)); - CHECK_FALSE(values.has_branch_from_power(0)); - CHECK_FALSE(values.has_branch_to_power(0)); - CHECK_FALSE(values.has_branch_to_current(0)); + MeasuredValues const values{std::make_shared(std::move(topo)), input}; - CHECK(values.branch_from_current(0).angle_measurement_type == local_angle); - check_accumulated(values.branch_from_current(0).measurement); - } + REQUIRE(values.has_branch_from_current(0)); + CHECK_FALSE(values.has_branch_from_power(0)); + CHECK_FALSE(values.has_branch_to_power(0)); + CHECK_FALSE(values.has_branch_to_current(0)); - SUBCASE("Accumulate global angle current sensors on branch - sym") { - auto topo = get_single_branch_topology(); - topo.current_sensors_per_branch_from = {from_dense, {0, 0}, 1}; - topo.voltage_sensors_per_bus = {from_dense, {0}, 2}; + CHECK(values.branch_from_current(0).angle_measurement_type == local_angle); + check_accumulated(values.branch_from_current(0).measurement); + } - StateEstimationInput input{}; - input.measured_voltage = {{.value = DoubleComplex{1.0, 1.0}, .variance = 10.0}}; - input.measured_branch_from_current = { - {.angle_measurement_type = global_angle, .measurement = measurement_a}, - {.angle_measurement_type = global_angle, .measurement = measurement_b}}; + SUBCASE("Accumulate global angle current sensors on branch") { + auto topo = get_single_branch_topology(); + topo.current_sensors_per_branch_from = {from_dense, {0, 0}, 1}; + topo.voltage_sensors_per_bus = {from_dense, {0}, 2}; + + StateEstimationInput input{}; + ComplexValue const measured_voltage = [] { + if constexpr (is_symmetric_v) { + return ComplexValue{1.0, 1.0}; + } else { + return ComplexValue{RealValue{1.0, 1.0, 1.0}, + RealValue{1.0, 1.0, 1.0}}; + } + }(); + input.measured_voltage = {{.value = measured_voltage, .variance = 10.0}}; + input.measured_branch_from_current = {{.angle_measurement_type = global_angle, .measurement = measurement_a}, + {.angle_measurement_type = global_angle, .measurement = measurement_b}}; + + MeasuredValues const values{std::make_shared(std::move(topo)), input}; + + REQUIRE(values.has_branch_from_current(0)); + CHECK_FALSE(values.has_branch_from_power(0)); + CHECK_FALSE(values.has_branch_to_power(0)); + CHECK_FALSE(values.has_branch_to_current(0)); + + CHECK(values.branch_from_current(0).angle_measurement_type == global_angle); + check_accumulated(values.branch_from_current(0).measurement); + } - MeasuredValues const values{std::make_shared(std::move(topo)), input}; + SUBCASE("Accumulate local angle current sensors on branch") { + auto topo = get_single_branch_topology(); + topo.current_sensors_per_branch_from = {from_dense, {0, 0}, 1}; - REQUIRE(values.has_branch_from_current(0)); - CHECK_FALSE(values.has_branch_from_power(0)); - CHECK_FALSE(values.has_branch_to_power(0)); - CHECK_FALSE(values.has_branch_to_current(0)); + StateEstimationInput input{}; + input.measured_branch_from_current = {{.angle_measurement_type = local_angle, .measurement = measurement_a}, + {.angle_measurement_type = local_angle, .measurement = measurement_b}}; - CHECK(values.branch_from_current(0).angle_measurement_type == global_angle); - check_accumulated(values.branch_from_current(0).measurement); - } + MeasuredValues const values{std::make_shared(std::move(topo)), input}; - SUBCASE("Cannot accumulate different angle measurement types on same terminal") { - auto topo = get_single_branch_topology(); - topo.current_sensors_per_branch_from = {from_dense, {0, 0}, 1}; + REQUIRE(values.has_branch_from_current(0)); + CHECK_FALSE(values.has_branch_from_power(0)); + CHECK_FALSE(values.has_branch_to_power(0)); + CHECK_FALSE(values.has_branch_to_current(0)); - StateEstimationInput input{}; - input.measured_branch_from_current = { - {.angle_measurement_type = AngleMeasurementType::local_angle, .measurement = measurement_a}, - {.angle_measurement_type = AngleMeasurementType::global_angle, .measurement = measurement_b}}; + CHECK(values.branch_from_current(0).angle_measurement_type == local_angle); + check_accumulated(values.branch_from_current(0).measurement); + } - auto const create_measured_values = [&topo, &input] { - return MeasuredValues{std::make_shared(topo), input}; - }; - CHECK_THROWS_AS(create_measured_values(), ConflictingAngleMeasurementType); - } + SUBCASE("Cannot accumulate different angle measurement types on same terminal") { + auto topo = get_single_branch_topology(); + topo.current_sensors_per_branch_from = {from_dense, {0, 0}, 1}; + + StateEstimationInput input{}; + input.measured_branch_from_current = { + {.angle_measurement_type = AngleMeasurementType::local_angle, .measurement = measurement_a}, + {.angle_measurement_type = AngleMeasurementType::global_angle, .measurement = measurement_b}}; + + auto const create_measured_values = [&topo, &input] { + return MeasuredValues{std::make_shared(topo), input}; + }; + CHECK_THROWS_AS(create_measured_values(), ConflictingAngleMeasurementType); } } From 1953c856b9357b88c679551065b17b9e43b7eecb Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Mon, 14 Apr 2025 08:42:46 +0200 Subject: [PATCH 22/25] resolve comments Signed-off-by: Martijn Govers --- .../power_grid_model/common/statistics.hpp | 8 +- .../math_solver/measured_values.hpp | 8 +- tests/cpp_unit_tests/test_measured_values.cpp | 1 + tests/cpp_unit_tests/test_statistics.cpp | 420 +++++++++--------- 4 files changed, 215 insertions(+), 222 deletions(-) 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 7b395635ed..c10ddd65fd 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 @@ -218,7 +218,7 @@ template UniformComplexRandVar::sym>> || std::same_as, IndependentComplexRandVar::sym>> -constexpr auto accumulate(RandVarsView rand_vars) { +constexpr auto combine(RandVarsView rand_vars) { using RandVarType = std::ranges::range_value_t; using ValueType = decltype(RandVarType::value); using VarianceType = decltype(RandVarType::variance); @@ -243,14 +243,14 @@ constexpr auto accumulate(RandVarsView rand_vars) { template requires std::same_as, DecomposedComplexRandVar::sym>> -constexpr auto accumulate(RandVarsView rand_vars) { +constexpr auto combine(RandVarsView rand_vars) { using sym = std::ranges::range_value_t::sym; DecomposedComplexRandVar result{ .real_component = - accumulate(rand_vars | std::views::transform([](auto const& x) -> auto& { return x.real_component; })), + combine(rand_vars | std::views::transform([](auto const& x) -> auto& { return x.real_component; })), .imag_component = - accumulate(rand_vars | std::views::transform([](auto const& x) -> auto& { return x.imag_component; }))}; + combine(rand_vars | std::views::transform([](auto const& x) -> auto& { return x.imag_component; }))}; if (!(is_normal(result.real_component.variance) && is_normal(result.imag_component.variance))) { result.real_component.variance = RealValue{std::numeric_limits::infinity()}; diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp index b45fb469f1..12d5f73102 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp @@ -455,7 +455,7 @@ template class MeasuredValues { IdxRange const& sensors) { auto complex_measurements = sensors | std::views::transform([&data](Idx pos) -> auto& { return data[pos]; }); if constexpr (only_magnitude) { - auto const weighted_average_magnitude_measurement = statistics::accumulate( + auto const weighted_average_magnitude_measurement = statistics::combine( complex_measurements | std::views::transform([](auto const& measurement) { return UniformRealRandVar{.value = detail::cabs_or_real(measurement.value), .variance = measurement.variance}; @@ -469,14 +469,14 @@ template class MeasuredValues { }(), .variance = weighted_average_magnitude_measurement.variance}; } else { - return statistics::accumulate(complex_measurements); + return statistics::combine(complex_measurements); } } template requires(!only_magnitude) static PowerSensorCalcParam combine_measurements(std::vector> const& data, IdxRange const& sensors) { - return statistics::accumulate(sensors | std::views::transform([&data](Idx pos) -> auto& { return data[pos]; })); + return statistics::combine(sensors | std::views::transform([&data](Idx pos) -> auto& { return data[pos]; })); } template requires(!only_magnitude) @@ -493,7 +493,7 @@ template class MeasuredValues { } return {.angle_measurement_type = angle_measurement_type, - .measurement = statistics::accumulate( + .measurement = statistics::combine( params | std::views::transform([](auto const& param) -> auto& { return param.measurement; }))}; } diff --git a/tests/cpp_unit_tests/test_measured_values.cpp b/tests/cpp_unit_tests/test_measured_values.cpp index 0c8183dcac..f63592acb6 100644 --- a/tests/cpp_unit_tests/test_measured_values.cpp +++ b/tests/cpp_unit_tests/test_measured_values.cpp @@ -186,6 +186,7 @@ TEST_CASE_TEMPLATE("Measured Values - Accumulate branch flow sensors", sym, symm check_close(value.imag_component.variance, expected_imag_variance); } else { for (Idx phase : IdxRange(3)) { + // eigen index-based element access check_close(value.value()(phase), expected_value); check_close(value.real_component.variance(phase), expected_real_variance); check_close(value.imag_component.variance(phase), expected_imag_variance); diff --git a/tests/cpp_unit_tests/test_statistics.cpp b/tests/cpp_unit_tests/test_statistics.cpp index 750a5f7e4e..a7ac5e1f2b 100644 --- a/tests/cpp_unit_tests/test_statistics.cpp +++ b/tests/cpp_unit_tests/test_statistics.cpp @@ -816,7 +816,7 @@ TEST_CASE("Test statistics") { } TEST_CASE("Test statistics - accumulate") { - using statistics::accumulate; + using statistics::combine; using std::views::take; SUBCASE("UniformRealRandVar | IndependentRealRandVar") { @@ -825,17 +825,17 @@ TEST_CASE("Test statistics - accumulate") { std::vector const measurements{ {.value = 1.0, .variance = 0.2}, {.value = 2.0, .variance = 0.3}, {.value = 3.0, .variance = 0.6}}; - CHECK(accumulate(measurements | take(0)).value == 0.0); - CHECK(is_inf(accumulate(measurements | take(0)).variance)); + CHECK(combine(measurements | take(0)).value == 0.0); + CHECK(is_inf(combine(measurements | take(0)).variance)); - CHECK(accumulate(measurements | take(1)).value == measurements.front().value); - CHECK(accumulate(measurements | take(1)).variance == measurements.front().variance); + CHECK(combine(measurements | take(1)).value == measurements.front().value); + CHECK(combine(measurements | take(1)).variance == measurements.front().variance); - CHECK(accumulate(measurements | take(2)).value == doctest::Approx(7.0 / 5.0)); - CHECK(accumulate(measurements | take(2)).variance == doctest::Approx(3.0 / 25.0)); + CHECK(combine(measurements | take(2)).value == doctest::Approx(7.0 / 5.0)); + CHECK(combine(measurements | take(2)).variance == doctest::Approx(3.0 / 25.0)); - CHECK(accumulate(measurements | take(3)).value == doctest::Approx(5.0 / 3.0)); - CHECK(accumulate(measurements | take(3)).variance == doctest::Approx(1.0 / 10.0)); + CHECK(combine(measurements | take(3)).value == doctest::Approx(5.0 / 3.0)); + CHECK(combine(measurements | take(3)).variance == doctest::Approx(1.0 / 10.0)); }; SUBCASE("UniformRealRandVar") { check.template operator()>(); } SUBCASE("IndependentRealRandVar") { @@ -848,25 +848,25 @@ TEST_CASE("Test statistics - accumulate") { {.value = {2.0, 4.0, 3.0}, .variance = 0.3}, {.value = {4.0, 5.0, 6.0}, .variance = 0.6}}; - CHECK(accumulate(measurements | take(0)).value(0) == 0.0); - CHECK(accumulate(measurements | take(0)).value(1) == 0.0); - CHECK(accumulate(measurements | take(0)).value(2) == 0.0); - CHECK(is_inf(accumulate(measurements | take(0)).variance)); - - CHECK(accumulate(measurements | take(1)).value(0) == measurements.front().value(0)); - CHECK(accumulate(measurements | take(1)).value(1) == measurements.front().value(1)); - CHECK(accumulate(measurements | take(1)).value(2) == measurements.front().value(2)); - CHECK(accumulate(measurements | take(1)).variance == measurements.front().variance); - - CHECK(accumulate(measurements | take(2)).value(0) == doctest::Approx(7.0 / 5.0)); - CHECK(accumulate(measurements | take(2)).value(1) == doctest::Approx(14.0 / 5.0)); - CHECK(accumulate(measurements | take(2)).value(2) == doctest::Approx(3.0 / 5.0)); - CHECK(accumulate(measurements | take(2)).variance == doctest::Approx(3.0 / 25.0)); - - CHECK(accumulate(measurements | take(3)).value(0) == doctest::Approx(11.0 / 6.0)); - CHECK(accumulate(measurements | take(3)).value(1) == doctest::Approx(19.0 / 6.0)); - CHECK(accumulate(measurements | take(3)).value(2) == doctest::Approx(3.0 / 2.0)); - CHECK(accumulate(measurements | take(3)).variance == doctest::Approx(1.0 / 10.0)); + CHECK(combine(measurements | take(0)).value(0) == 0.0); + CHECK(combine(measurements | take(0)).value(1) == 0.0); + CHECK(combine(measurements | take(0)).value(2) == 0.0); + CHECK(is_inf(combine(measurements | take(0)).variance)); + + CHECK(combine(measurements | take(1)).value(0) == measurements.front().value(0)); + CHECK(combine(measurements | take(1)).value(1) == measurements.front().value(1)); + CHECK(combine(measurements | take(1)).value(2) == measurements.front().value(2)); + CHECK(combine(measurements | take(1)).variance == measurements.front().variance); + + CHECK(combine(measurements | take(2)).value(0) == doctest::Approx(7.0 / 5.0)); + CHECK(combine(measurements | take(2)).value(1) == doctest::Approx(14.0 / 5.0)); + CHECK(combine(measurements | take(2)).value(2) == doctest::Approx(3.0 / 5.0)); + CHECK(combine(measurements | take(2)).variance == doctest::Approx(3.0 / 25.0)); + + CHECK(combine(measurements | take(3)).value(0) == doctest::Approx(11.0 / 6.0)); + CHECK(combine(measurements | take(3)).value(1) == doctest::Approx(19.0 / 6.0)); + CHECK(combine(measurements | take(3)).value(2) == doctest::Approx(3.0 / 2.0)); + CHECK(combine(measurements | take(3)).variance == doctest::Approx(1.0 / 10.0)); } SUBCASE("IndependentRealRandVar") { @@ -875,33 +875,33 @@ TEST_CASE("Test statistics - accumulate") { {.value = {2.0, 4.0, 3.0}, .variance = {0.3, 0.4, 0.5}}, {.value = {4.0, 5.0, 6.0}, .variance = {0.6, 0.7, 0.8}}}; - CHECK(accumulate(measurements | take(0)).value(0) == 0.0); - CHECK(accumulate(measurements | take(0)).value(1) == 0.0); - CHECK(accumulate(measurements | take(0)).value(2) == 0.0); - CHECK(is_inf(accumulate(measurements | take(0)).variance(0))); - CHECK(is_inf(accumulate(measurements | take(0)).variance(1))); - CHECK(is_inf(accumulate(measurements | take(0)).variance(2))); - - CHECK(accumulate(measurements | take(1)).value(0) == measurements.front().value(0)); - CHECK(accumulate(measurements | take(1)).value(1) == measurements.front().value(1)); - CHECK(accumulate(measurements | take(1)).value(2) == measurements.front().value(2)); - CHECK(accumulate(measurements | take(1)).variance(0) == measurements.front().variance(0)); - CHECK(accumulate(measurements | take(1)).variance(1) == measurements.front().variance(1)); - CHECK(accumulate(measurements | take(1)).variance(2) == measurements.front().variance(2)); - - CHECK(accumulate(measurements | take(2)).value(0) == doctest::Approx(7.0 / 5.0)); - CHECK(accumulate(measurements | take(2)).value(1) == doctest::Approx(20.0 / 7.0)); - CHECK(accumulate(measurements | take(2)).value(2) == doctest::Approx(7.0 / 9.0)); - CHECK(accumulate(measurements | take(2)).variance(0) == doctest::Approx(3.0 / 25.0)); - CHECK(accumulate(measurements | take(2)).variance(1) == doctest::Approx(6.0 / 35.0)); - CHECK(accumulate(measurements | take(2)).variance(2) == doctest::Approx(2.0 / 9.0)); - - CHECK(accumulate(measurements | take(3)).value(0) == doctest::Approx(11.0 / 6.0)); - CHECK(accumulate(measurements | take(3)).value(1) == doctest::Approx(200.0 / 61.0)); - CHECK(accumulate(measurements | take(3)).value(2) == doctest::Approx(44.0 / 23.0)); - CHECK(accumulate(measurements | take(3)).variance(0) == doctest::Approx(1.0 / 10.0)); - CHECK(accumulate(measurements | take(3)).variance(1) == doctest::Approx(42.0 / 305.0)); - CHECK(accumulate(measurements | take(3)).variance(2) == doctest::Approx(4.0 / 23.0)); + CHECK(combine(measurements | take(0)).value(0) == 0.0); + CHECK(combine(measurements | take(0)).value(1) == 0.0); + CHECK(combine(measurements | take(0)).value(2) == 0.0); + CHECK(is_inf(combine(measurements | take(0)).variance(0))); + CHECK(is_inf(combine(measurements | take(0)).variance(1))); + CHECK(is_inf(combine(measurements | take(0)).variance(2))); + + CHECK(combine(measurements | take(1)).value(0) == measurements.front().value(0)); + CHECK(combine(measurements | take(1)).value(1) == measurements.front().value(1)); + CHECK(combine(measurements | take(1)).value(2) == measurements.front().value(2)); + CHECK(combine(measurements | take(1)).variance(0) == measurements.front().variance(0)); + CHECK(combine(measurements | take(1)).variance(1) == measurements.front().variance(1)); + CHECK(combine(measurements | take(1)).variance(2) == measurements.front().variance(2)); + + CHECK(combine(measurements | take(2)).value(0) == doctest::Approx(7.0 / 5.0)); + CHECK(combine(measurements | take(2)).value(1) == doctest::Approx(20.0 / 7.0)); + CHECK(combine(measurements | take(2)).value(2) == doctest::Approx(7.0 / 9.0)); + CHECK(combine(measurements | take(2)).variance(0) == doctest::Approx(3.0 / 25.0)); + CHECK(combine(measurements | take(2)).variance(1) == doctest::Approx(6.0 / 35.0)); + CHECK(combine(measurements | take(2)).variance(2) == doctest::Approx(2.0 / 9.0)); + + CHECK(combine(measurements | take(3)).value(0) == doctest::Approx(11.0 / 6.0)); + CHECK(combine(measurements | take(3)).value(1) == doctest::Approx(200.0 / 61.0)); + CHECK(combine(measurements | take(3)).value(2) == doctest::Approx(44.0 / 23.0)); + CHECK(combine(measurements | take(3)).variance(0) == doctest::Approx(1.0 / 10.0)); + CHECK(combine(measurements | take(3)).variance(1) == doctest::Approx(42.0 / 305.0)); + CHECK(combine(measurements | take(3)).variance(2) == doctest::Approx(4.0 / 23.0)); } SUBCASE("UniformComplexRandVar | IndependentComplexRandVar") { @@ -911,21 +911,21 @@ TEST_CASE("Test statistics - accumulate") { T{.value = 2.0 + 6.0i, .variance = 0.3}, T{.value = 4.0 + 3.0i, .variance = 0.6}}; - CHECK(accumulate(measurements | take(0)).value.real() == 0.0); - CHECK(accumulate(measurements | take(0)).value.imag() == 0.0); - CHECK(is_inf(accumulate(measurements | take(0)).variance)); + CHECK(combine(measurements | take(0)).value.real() == 0.0); + CHECK(combine(measurements | take(0)).value.imag() == 0.0); + CHECK(is_inf(combine(measurements | take(0)).variance)); - CHECK(accumulate(measurements | take(1)).value.real() == measurements.front().value.real()); - CHECK(accumulate(measurements | take(1)).value.imag() == measurements.front().value.imag()); - CHECK(accumulate(measurements | take(1)).variance == measurements.front().variance); + CHECK(combine(measurements | take(1)).value.real() == measurements.front().value.real()); + CHECK(combine(measurements | take(1)).value.imag() == measurements.front().value.imag()); + CHECK(combine(measurements | take(1)).variance == measurements.front().variance); - CHECK(accumulate(measurements | take(2)).value.real() == doctest::Approx(7.0 / 5.0)); - CHECK(accumulate(measurements | take(2)).value.imag() == doctest::Approx(27.0 / 5.0)); - CHECK(accumulate(measurements | take(2)).variance == doctest::Approx(3.0 / 25.0)); + CHECK(combine(measurements | take(2)).value.real() == doctest::Approx(7.0 / 5.0)); + CHECK(combine(measurements | take(2)).value.imag() == doctest::Approx(27.0 / 5.0)); + CHECK(combine(measurements | take(2)).variance == doctest::Approx(3.0 / 25.0)); - CHECK(accumulate(measurements | take(3)).value.real() == doctest::Approx(11.0 / 6.0)); - CHECK(accumulate(measurements | take(3)).value.imag() == doctest::Approx(30.0 / 6.0)); - CHECK(accumulate(measurements | take(3)).variance == doctest::Approx(1.0 / 10.0)); + CHECK(combine(measurements | take(3)).value.real() == doctest::Approx(11.0 / 6.0)); + CHECK(combine(measurements | take(3)).value.imag() == doctest::Approx(30.0 / 6.0)); + CHECK(combine(measurements | take(3)).variance == doctest::Approx(1.0 / 10.0)); }; SUBCASE("UniformComplexRandVar") { check.template operator()>(); @@ -944,37 +944,37 @@ TEST_CASE("Test statistics - accumulate") { {.value = {RealValue{4.0, 5.0, 6.0}, RealValue{3.0, 1.0, 2.0}}, .variance = 0.6}}; - CHECK(accumulate(measurements | take(0)).value(0).real() == 0.0); - CHECK(accumulate(measurements | take(0)).value(1).real() == 0.0); - CHECK(accumulate(measurements | take(0)).value(2).real() == 0.0); - CHECK(accumulate(measurements | take(0)).value(0).imag() == 0.0); - CHECK(accumulate(measurements | take(0)).value(1).imag() == 0.0); - CHECK(accumulate(measurements | take(0)).value(2).imag() == 0.0); - CHECK(is_inf(accumulate(measurements | take(0)).variance)); - - CHECK(accumulate(measurements | take(1)).value(0).real() == measurements.front().value(0).real()); - CHECK(accumulate(measurements | take(1)).value(1).real() == measurements.front().value(1).real()); - CHECK(accumulate(measurements | take(1)).value(2).real() == measurements.front().value(2).real()); - CHECK(accumulate(measurements | take(1)).value(0).imag() == measurements.front().value(0).imag()); - CHECK(accumulate(measurements | take(1)).value(1).imag() == measurements.front().value(1).imag()); - CHECK(accumulate(measurements | take(1)).value(2).imag() == measurements.front().value(2).imag()); - CHECK(accumulate(measurements | take(1)).variance == measurements.front().variance); - - CHECK(accumulate(measurements | take(2)).value(0).real() == doctest::Approx(7.0 / 5.0)); - CHECK(accumulate(measurements | take(2)).value(1).real() == doctest::Approx(14.0 / 5.0)); - CHECK(accumulate(measurements | take(2)).value(2).real() == doctest::Approx(3.0 / 5.0)); - CHECK(accumulate(measurements | take(2)).value(0).imag() == doctest::Approx(27.0 / 5.0)); - CHECK(accumulate(measurements | take(2)).value(1).imag() == doctest::Approx(4.0 / 5.0)); - CHECK(accumulate(measurements | take(2)).value(2).imag() == doctest::Approx(25.0 / 5.0)); - CHECK(accumulate(measurements | take(2)).variance == doctest::Approx(3.0 / 25.0)); - - CHECK(accumulate(measurements | take(3)).value(0).real() == doctest::Approx(11.0 / 6.0)); - CHECK(accumulate(measurements | take(3)).value(1).real() == doctest::Approx(19.0 / 6.0)); - CHECK(accumulate(measurements | take(3)).value(2).real() == doctest::Approx(9.0 / 6.0)); - CHECK(accumulate(measurements | take(3)).value(0).imag() == doctest::Approx(30.0 / 6.0)); - CHECK(accumulate(measurements | take(3)).value(1).imag() == doctest::Approx(5.0 / 6.0)); - CHECK(accumulate(measurements | take(3)).value(2).imag() == doctest::Approx(27.0 / 6.0)); - CHECK(accumulate(measurements | take(3)).variance == doctest::Approx(1.0 / 10.0)); + CHECK(combine(measurements | take(0)).value(0).real() == 0.0); + CHECK(combine(measurements | take(0)).value(1).real() == 0.0); + CHECK(combine(measurements | take(0)).value(2).real() == 0.0); + CHECK(combine(measurements | take(0)).value(0).imag() == 0.0); + CHECK(combine(measurements | take(0)).value(1).imag() == 0.0); + CHECK(combine(measurements | take(0)).value(2).imag() == 0.0); + CHECK(is_inf(combine(measurements | take(0)).variance)); + + CHECK(combine(measurements | take(1)).value(0).real() == measurements.front().value(0).real()); + CHECK(combine(measurements | take(1)).value(1).real() == measurements.front().value(1).real()); + CHECK(combine(measurements | take(1)).value(2).real() == measurements.front().value(2).real()); + CHECK(combine(measurements | take(1)).value(0).imag() == measurements.front().value(0).imag()); + CHECK(combine(measurements | take(1)).value(1).imag() == measurements.front().value(1).imag()); + CHECK(combine(measurements | take(1)).value(2).imag() == measurements.front().value(2).imag()); + CHECK(combine(measurements | take(1)).variance == measurements.front().variance); + + CHECK(combine(measurements | take(2)).value(0).real() == doctest::Approx(7.0 / 5.0)); + CHECK(combine(measurements | take(2)).value(1).real() == doctest::Approx(14.0 / 5.0)); + CHECK(combine(measurements | take(2)).value(2).real() == doctest::Approx(3.0 / 5.0)); + CHECK(combine(measurements | take(2)).value(0).imag() == doctest::Approx(27.0 / 5.0)); + CHECK(combine(measurements | take(2)).value(1).imag() == doctest::Approx(4.0 / 5.0)); + CHECK(combine(measurements | take(2)).value(2).imag() == doctest::Approx(25.0 / 5.0)); + CHECK(combine(measurements | take(2)).variance == doctest::Approx(3.0 / 25.0)); + + CHECK(combine(measurements | take(3)).value(0).real() == doctest::Approx(11.0 / 6.0)); + CHECK(combine(measurements | take(3)).value(1).real() == doctest::Approx(19.0 / 6.0)); + CHECK(combine(measurements | take(3)).value(2).real() == doctest::Approx(9.0 / 6.0)); + CHECK(combine(measurements | take(3)).value(0).imag() == doctest::Approx(30.0 / 6.0)); + CHECK(combine(measurements | take(3)).value(1).imag() == doctest::Approx(5.0 / 6.0)); + CHECK(combine(measurements | take(3)).value(2).imag() == doctest::Approx(27.0 / 6.0)); + CHECK(combine(measurements | take(3)).variance == doctest::Approx(1.0 / 10.0)); } SUBCASE("IndependentComplexRandVar") { @@ -986,45 +986,45 @@ TEST_CASE("Test statistics - accumulate") { {.value = {RealValue{4.0, 5.0, 6.0}, RealValue{3.0, 1.0, 2.0}}, .variance = {0.6, 0.7, 0.8}}}; - CHECK(accumulate(measurements | take(0)).value(0).real() == 0.0); - CHECK(accumulate(measurements | take(0)).value(1).real() == 0.0); - CHECK(accumulate(measurements | take(0)).value(2).real() == 0.0); - CHECK(accumulate(measurements | take(0)).value(0).imag() == 0.0); - CHECK(accumulate(measurements | take(0)).value(1).imag() == 0.0); - CHECK(accumulate(measurements | take(0)).value(2).imag() == 0.0); - CHECK(is_inf(accumulate(measurements | take(0)).variance(0))); - CHECK(is_inf(accumulate(measurements | take(0)).variance(1))); - CHECK(is_inf(accumulate(measurements | take(0)).variance(2))); - - CHECK(accumulate(measurements | take(1)).value(0).real() == measurements.front().value(0).real()); - CHECK(accumulate(measurements | take(1)).value(1).real() == measurements.front().value(1).real()); - CHECK(accumulate(measurements | take(1)).value(2).real() == measurements.front().value(2).real()); - CHECK(accumulate(measurements | take(1)).value(0).imag() == measurements.front().value(0).imag()); - CHECK(accumulate(measurements | take(1)).value(1).imag() == measurements.front().value(1).imag()); - CHECK(accumulate(measurements | take(1)).value(2).imag() == measurements.front().value(2).imag()); - CHECK(accumulate(measurements | take(1)).variance(0) == measurements.front().variance(0)); - CHECK(accumulate(measurements | take(1)).variance(1) == measurements.front().variance(1)); - CHECK(accumulate(measurements | take(1)).variance(2) == measurements.front().variance(2)); - - CHECK(accumulate(measurements | take(2)).value(0).real() == doctest::Approx(7.0 / 5.0)); - CHECK(accumulate(measurements | take(2)).value(1).real() == doctest::Approx(20.0 / 7.0)); - CHECK(accumulate(measurements | take(2)).value(2).real() == doctest::Approx(7.0 / 9.0)); - CHECK(accumulate(measurements | take(2)).value(0).imag() == doctest::Approx(27.0 / 5.0)); - CHECK(accumulate(measurements | take(2)).value(1).imag() == doctest::Approx(3.0 / 7.0)); - CHECK(accumulate(measurements | take(2)).value(2).imag() == doctest::Approx(43.0 / 9.0)); - CHECK(accumulate(measurements | take(2)).variance(0) == doctest::Approx(3.0 / 25.0)); - CHECK(accumulate(measurements | take(2)).variance(1) == doctest::Approx(6.0 / 35.0)); - CHECK(accumulate(measurements | take(2)).variance(2) == doctest::Approx(2.0 / 9.0)); - - CHECK(accumulate(measurements | take(3)).value(0).real() == doctest::Approx(11.0 / 6.0)); - CHECK(accumulate(measurements | take(3)).value(1).real() == doctest::Approx(200.0 / 61.0)); - CHECK(accumulate(measurements | take(3)).value(2).real() == doctest::Approx(44.0 / 23.0)); - CHECK(accumulate(measurements | take(3)).value(0).imag() == doctest::Approx(30.0 / 6.0)); - CHECK(accumulate(measurements | take(3)).value(1).imag() == doctest::Approx(33.0 / 61.0)); - CHECK(accumulate(measurements | take(3)).value(2).imag() == doctest::Approx(96.0 / 23.0)); - CHECK(accumulate(measurements | take(3)).variance(0) == doctest::Approx(1.0 / 10.0)); - CHECK(accumulate(measurements | take(3)).variance(1) == doctest::Approx(42.0 / 305.0)); - CHECK(accumulate(measurements | take(3)).variance(2) == doctest::Approx(4.0 / 23.0)); + CHECK(combine(measurements | take(0)).value(0).real() == 0.0); + CHECK(combine(measurements | take(0)).value(1).real() == 0.0); + CHECK(combine(measurements | take(0)).value(2).real() == 0.0); + CHECK(combine(measurements | take(0)).value(0).imag() == 0.0); + CHECK(combine(measurements | take(0)).value(1).imag() == 0.0); + CHECK(combine(measurements | take(0)).value(2).imag() == 0.0); + CHECK(is_inf(combine(measurements | take(0)).variance(0))); + CHECK(is_inf(combine(measurements | take(0)).variance(1))); + CHECK(is_inf(combine(measurements | take(0)).variance(2))); + + CHECK(combine(measurements | take(1)).value(0).real() == measurements.front().value(0).real()); + CHECK(combine(measurements | take(1)).value(1).real() == measurements.front().value(1).real()); + CHECK(combine(measurements | take(1)).value(2).real() == measurements.front().value(2).real()); + CHECK(combine(measurements | take(1)).value(0).imag() == measurements.front().value(0).imag()); + CHECK(combine(measurements | take(1)).value(1).imag() == measurements.front().value(1).imag()); + CHECK(combine(measurements | take(1)).value(2).imag() == measurements.front().value(2).imag()); + CHECK(combine(measurements | take(1)).variance(0) == measurements.front().variance(0)); + CHECK(combine(measurements | take(1)).variance(1) == measurements.front().variance(1)); + CHECK(combine(measurements | take(1)).variance(2) == measurements.front().variance(2)); + + CHECK(combine(measurements | take(2)).value(0).real() == doctest::Approx(7.0 / 5.0)); + CHECK(combine(measurements | take(2)).value(1).real() == doctest::Approx(20.0 / 7.0)); + CHECK(combine(measurements | take(2)).value(2).real() == doctest::Approx(7.0 / 9.0)); + CHECK(combine(measurements | take(2)).value(0).imag() == doctest::Approx(27.0 / 5.0)); + CHECK(combine(measurements | take(2)).value(1).imag() == doctest::Approx(3.0 / 7.0)); + CHECK(combine(measurements | take(2)).value(2).imag() == doctest::Approx(43.0 / 9.0)); + CHECK(combine(measurements | take(2)).variance(0) == doctest::Approx(3.0 / 25.0)); + CHECK(combine(measurements | take(2)).variance(1) == doctest::Approx(6.0 / 35.0)); + CHECK(combine(measurements | take(2)).variance(2) == doctest::Approx(2.0 / 9.0)); + + CHECK(combine(measurements | take(3)).value(0).real() == doctest::Approx(11.0 / 6.0)); + CHECK(combine(measurements | take(3)).value(1).real() == doctest::Approx(200.0 / 61.0)); + CHECK(combine(measurements | take(3)).value(2).real() == doctest::Approx(44.0 / 23.0)); + CHECK(combine(measurements | take(3)).value(0).imag() == doctest::Approx(30.0 / 6.0)); + CHECK(combine(measurements | take(3)).value(1).imag() == doctest::Approx(33.0 / 61.0)); + CHECK(combine(measurements | take(3)).value(2).imag() == doctest::Approx(96.0 / 23.0)); + CHECK(combine(measurements | take(3)).variance(0) == doctest::Approx(1.0 / 10.0)); + CHECK(combine(measurements | take(3)).variance(1) == doctest::Approx(42.0 / 305.0)); + CHECK(combine(measurements | take(3)).variance(2) == doctest::Approx(4.0 / 23.0)); } SUBCASE("DecomposedComplexRandVar") { @@ -1036,27 +1036,25 @@ TEST_CASE("Test statistics - accumulate") { DecomposedComplexRandVar{.real_component = {.value = 4.0, .variance = 0.6}, .imag_component = {.value = 3.0, .variance = 0.3}}}; - CHECK(accumulate(measurements | take(0)).real_component.value == 0.0); - CHECK(accumulate(measurements | take(0)).imag_component.value == 0.0); - CHECK(is_inf(accumulate(measurements | take(0)).real_component.variance)); - CHECK(is_inf(accumulate(measurements | take(0)).imag_component.variance)); - - CHECK(accumulate(measurements | take(1)).real_component.value == measurements.front().real_component.value); - CHECK(accumulate(measurements | take(1)).imag_component.value == measurements.front().imag_component.value); - CHECK(accumulate(measurements | take(1)).real_component.variance == - measurements.front().real_component.variance); - CHECK(accumulate(measurements | take(1)).imag_component.variance == - measurements.front().imag_component.variance); - - CHECK(accumulate(measurements | take(2)).real_component.value == doctest::Approx(7.0 / 5.0)); - CHECK(accumulate(measurements | take(2)).imag_component.value == doctest::Approx(80.0 / 15.0)); - CHECK(accumulate(measurements | take(2)).real_component.variance == doctest::Approx(3.0 / 25.0)); - CHECK(accumulate(measurements | take(2)).imag_component.variance == doctest::Approx(1.0 / 15.0)); - - CHECK(accumulate(measurements | take(3)).real_component.value == doctest::Approx(11.0 / 6.0)); - CHECK(accumulate(measurements | take(3)).imag_component.value == doctest::Approx(270.0 / 55.0)); - CHECK(accumulate(measurements | take(3)).real_component.variance == doctest::Approx(1.0 / 10.0)); - CHECK(accumulate(measurements | take(3)).imag_component.variance == doctest::Approx(3.0 / 55.0)); + CHECK(combine(measurements | take(0)).real_component.value == 0.0); + CHECK(combine(measurements | take(0)).imag_component.value == 0.0); + CHECK(is_inf(combine(measurements | take(0)).real_component.variance)); + CHECK(is_inf(combine(measurements | take(0)).imag_component.variance)); + + CHECK(combine(measurements | take(1)).real_component.value == measurements.front().real_component.value); + CHECK(combine(measurements | take(1)).imag_component.value == measurements.front().imag_component.value); + CHECK(combine(measurements | take(1)).real_component.variance == measurements.front().real_component.variance); + CHECK(combine(measurements | take(1)).imag_component.variance == measurements.front().imag_component.variance); + + CHECK(combine(measurements | take(2)).real_component.value == doctest::Approx(7.0 / 5.0)); + CHECK(combine(measurements | take(2)).imag_component.value == doctest::Approx(80.0 / 15.0)); + CHECK(combine(measurements | take(2)).real_component.variance == doctest::Approx(3.0 / 25.0)); + CHECK(combine(measurements | take(2)).imag_component.variance == doctest::Approx(1.0 / 15.0)); + + CHECK(combine(measurements | take(3)).real_component.value == doctest::Approx(11.0 / 6.0)); + CHECK(combine(measurements | take(3)).imag_component.value == doctest::Approx(270.0 / 55.0)); + CHECK(combine(measurements | take(3)).real_component.variance == doctest::Approx(1.0 / 10.0)); + CHECK(combine(measurements | take(3)).imag_component.variance == doctest::Approx(3.0 / 55.0)); } SUBCASE("DecomposedComplexRandVar") { @@ -1071,69 +1069,63 @@ TEST_CASE("Test statistics - accumulate") { .real_component = {.value = RealValue{4.0, 5.0, 6.0}, .variance = {0.6, 0.7, 0.8}}, .imag_component = {.value = RealValue{3.0, 1.0, 2.0}, .variance = {0.3, 0.4, 0.5}}}}; - CHECK(accumulate(measurements | take(0)).real_component.value(0) == 0.0); - CHECK(accumulate(measurements | take(0)).real_component.value(2) == 0.0); - CHECK(accumulate(measurements | take(0)).real_component.value(1) == 0.0); - CHECK(accumulate(measurements | take(0)).imag_component.value(0) == 0.0); - CHECK(accumulate(measurements | take(0)).imag_component.value(2) == 0.0); - CHECK(accumulate(measurements | take(0)).imag_component.value(1) == 0.0); - CHECK(is_inf(accumulate(measurements | take(0)).real_component.variance(0))); - CHECK(is_inf(accumulate(measurements | take(0)).real_component.variance(1))); - CHECK(is_inf(accumulate(measurements | take(0)).real_component.variance(2))); - CHECK(is_inf(accumulate(measurements | take(0)).imag_component.variance(0))); - CHECK(is_inf(accumulate(measurements | take(0)).imag_component.variance(1))); - CHECK(is_inf(accumulate(measurements | take(0)).imag_component.variance(2))); - - CHECK(accumulate(measurements | take(1)).real_component.value(0) == - measurements.front().real_component.value(0)); - CHECK(accumulate(measurements | take(1)).real_component.value(1) == - measurements.front().real_component.value(1)); - CHECK(accumulate(measurements | take(1)).real_component.value(2) == - measurements.front().real_component.value(2)); - CHECK(accumulate(measurements | take(1)).imag_component.value(0) == - measurements.front().imag_component.value(0)); - CHECK(accumulate(measurements | take(1)).imag_component.value(1) == - measurements.front().imag_component.value(1)); - CHECK(accumulate(measurements | take(1)).imag_component.value(2) == - measurements.front().imag_component.value(2)); - CHECK(accumulate(measurements | take(1)).real_component.variance(0) == + CHECK(combine(measurements | take(0)).real_component.value(0) == 0.0); + CHECK(combine(measurements | take(0)).real_component.value(2) == 0.0); + CHECK(combine(measurements | take(0)).real_component.value(1) == 0.0); + CHECK(combine(measurements | take(0)).imag_component.value(0) == 0.0); + CHECK(combine(measurements | take(0)).imag_component.value(2) == 0.0); + CHECK(combine(measurements | take(0)).imag_component.value(1) == 0.0); + CHECK(is_inf(combine(measurements | take(0)).real_component.variance(0))); + CHECK(is_inf(combine(measurements | take(0)).real_component.variance(1))); + CHECK(is_inf(combine(measurements | take(0)).real_component.variance(2))); + CHECK(is_inf(combine(measurements | take(0)).imag_component.variance(0))); + CHECK(is_inf(combine(measurements | take(0)).imag_component.variance(1))); + CHECK(is_inf(combine(measurements | take(0)).imag_component.variance(2))); + + CHECK(combine(measurements | take(1)).real_component.value(0) == measurements.front().real_component.value(0)); + CHECK(combine(measurements | take(1)).real_component.value(1) == measurements.front().real_component.value(1)); + CHECK(combine(measurements | take(1)).real_component.value(2) == measurements.front().real_component.value(2)); + CHECK(combine(measurements | take(1)).imag_component.value(0) == measurements.front().imag_component.value(0)); + CHECK(combine(measurements | take(1)).imag_component.value(1) == measurements.front().imag_component.value(1)); + CHECK(combine(measurements | take(1)).imag_component.value(2) == measurements.front().imag_component.value(2)); + CHECK(combine(measurements | take(1)).real_component.variance(0) == measurements.front().real_component.variance(0)); - CHECK(accumulate(measurements | take(1)).real_component.variance(1) == + CHECK(combine(measurements | take(1)).real_component.variance(1) == measurements.front().real_component.variance(1)); - CHECK(accumulate(measurements | take(1)).real_component.variance(2) == + CHECK(combine(measurements | take(1)).real_component.variance(2) == measurements.front().real_component.variance(2)); - CHECK(accumulate(measurements | take(1)).imag_component.variance(0) == + CHECK(combine(measurements | take(1)).imag_component.variance(0) == measurements.front().imag_component.variance(0)); - CHECK(accumulate(measurements | take(1)).imag_component.variance(1) == + CHECK(combine(measurements | take(1)).imag_component.variance(1) == measurements.front().imag_component.variance(1)); - CHECK(accumulate(measurements | take(1)).imag_component.variance(2) == + CHECK(combine(measurements | take(1)).imag_component.variance(2) == measurements.front().imag_component.variance(2)); - CHECK(accumulate(measurements | take(2)).real_component.value(0) == doctest::Approx(7.0 / 5.0)); - CHECK(accumulate(measurements | take(2)).real_component.value(1) == doctest::Approx(20.0 / 7.0)); - CHECK(accumulate(measurements | take(2)).real_component.value(2) == doctest::Approx(7.0 / 9.0)); - CHECK(accumulate(measurements | take(2)).imag_component.value(0) == doctest::Approx(80.0 / 15.0)); - CHECK(accumulate(measurements | take(2)).imag_component.value(1) == doctest::Approx(20 / 25.0)); - CHECK(accumulate(measurements | take(2)).imag_component.value(2) == doctest::Approx(170.0 / 35.0)); - CHECK(accumulate(measurements | take(2)).real_component.variance(0) == doctest::Approx(3.0 / 25.0)); - CHECK(accumulate(measurements | take(2)).real_component.variance(1) == doctest::Approx(6.0 / 35.0)); - CHECK(accumulate(measurements | take(2)).real_component.variance(2) == doctest::Approx(2.0 / 9.0)); - CHECK(accumulate(measurements | take(2)).imag_component.variance(0) == doctest::Approx(1.0 / 15.0)); - CHECK(accumulate(measurements | take(2)).imag_component.variance(1) == doctest::Approx(3.0 / 25.0)); - CHECK(accumulate(measurements | take(2)).imag_component.variance(2) == doctest::Approx(6.0 / 35.0)); - - CHECK(accumulate(measurements | take(3)).real_component.value(0) == doctest::Approx(11.0 / 6.0)); - CHECK(accumulate(measurements | take(3)).real_component.value(1) == doctest::Approx(200.0 / 61.0)); - CHECK(accumulate(measurements | take(3)).real_component.value(2) == doctest::Approx(44.0 / 23.0)); - CHECK(accumulate(measurements | take(3)).imag_component.value(0) == doctest::Approx(270.0 / 55.0)); - CHECK(accumulate(measurements | take(3)).imag_component.value(1) == doctest::Approx(55.0 / 65.0)); - CHECK(accumulate(measurements | take(3)).imag_component.value(2) == doctest::Approx(194.0 / 47.0)); - CHECK(accumulate(measurements | take(3)).real_component.variance(0) == doctest::Approx(1.0 / 10.0)); - CHECK(accumulate(measurements | take(3)).real_component.variance(1) == doctest::Approx(42.0 / 305.0)); - CHECK(accumulate(measurements | take(3)).real_component.variance(2) == doctest::Approx(4.0 / 23.0)); - CHECK(accumulate(measurements | take(3)).imag_component.variance(0) == doctest::Approx(3.0 / 55.0)); - CHECK(accumulate(measurements | take(3)).imag_component.variance(1) == doctest::Approx(6.0 / 65)); - CHECK(accumulate(measurements | take(3)).imag_component.variance(2) == doctest::Approx(6.0 / 47.0)); + CHECK(combine(measurements | take(2)).real_component.value(0) == doctest::Approx(7.0 / 5.0)); + CHECK(combine(measurements | take(2)).real_component.value(1) == doctest::Approx(20.0 / 7.0)); + CHECK(combine(measurements | take(2)).real_component.value(2) == doctest::Approx(7.0 / 9.0)); + CHECK(combine(measurements | take(2)).imag_component.value(0) == doctest::Approx(80.0 / 15.0)); + CHECK(combine(measurements | take(2)).imag_component.value(1) == doctest::Approx(20 / 25.0)); + CHECK(combine(measurements | take(2)).imag_component.value(2) == doctest::Approx(170.0 / 35.0)); + CHECK(combine(measurements | take(2)).real_component.variance(0) == doctest::Approx(3.0 / 25.0)); + CHECK(combine(measurements | take(2)).real_component.variance(1) == doctest::Approx(6.0 / 35.0)); + CHECK(combine(measurements | take(2)).real_component.variance(2) == doctest::Approx(2.0 / 9.0)); + CHECK(combine(measurements | take(2)).imag_component.variance(0) == doctest::Approx(1.0 / 15.0)); + CHECK(combine(measurements | take(2)).imag_component.variance(1) == doctest::Approx(3.0 / 25.0)); + CHECK(combine(measurements | take(2)).imag_component.variance(2) == doctest::Approx(6.0 / 35.0)); + + CHECK(combine(measurements | take(3)).real_component.value(0) == doctest::Approx(11.0 / 6.0)); + CHECK(combine(measurements | take(3)).real_component.value(1) == doctest::Approx(200.0 / 61.0)); + CHECK(combine(measurements | take(3)).real_component.value(2) == doctest::Approx(44.0 / 23.0)); + CHECK(combine(measurements | take(3)).imag_component.value(0) == doctest::Approx(270.0 / 55.0)); + CHECK(combine(measurements | take(3)).imag_component.value(1) == doctest::Approx(55.0 / 65.0)); + CHECK(combine(measurements | take(3)).imag_component.value(2) == doctest::Approx(194.0 / 47.0)); + CHECK(combine(measurements | take(3)).real_component.variance(0) == doctest::Approx(1.0 / 10.0)); + CHECK(combine(measurements | take(3)).real_component.variance(1) == doctest::Approx(42.0 / 305.0)); + CHECK(combine(measurements | take(3)).real_component.variance(2) == doctest::Approx(4.0 / 23.0)); + CHECK(combine(measurements | take(3)).imag_component.variance(0) == doctest::Approx(3.0 / 55.0)); + CHECK(combine(measurements | take(3)).imag_component.variance(1) == doctest::Approx(6.0 / 65)); + CHECK(combine(measurements | take(3)).imag_component.variance(2) == doctest::Approx(6.0 / 47.0)); } } From 04e1273e5bbbca2643dcec8f93ac1a02d4d5cd2c Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Mon, 14 Apr 2025 10:10:22 +0200 Subject: [PATCH 23/25] also migrate combine_magnitude to statistics Signed-off-by: Martijn Govers --- .../power_grid_model/common/statistics.hpp | 46 +++++++++++ .../math_solver/measured_values.hpp | 24 +----- .../math_solver/newton_raphson_se_solver.hpp | 8 +- tests/cpp_unit_tests/test_statistics.cpp | 77 ++++++++++++++++++- 4 files changed, 129 insertions(+), 26 deletions(-) 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 c10ddd65fd..17e1678823 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 @@ -258,5 +258,51 @@ constexpr auto combine(RandVarsView rand_vars) { } return result; } + +namespace detail { +template inline RealValue cabs_or_real(ComplexValue const& value) { + if (is_nan(imag(value))) { + return real(value); // only keep real part + } + return cabs(value); // get abs of the value +} +} // namespace detail + +template + requires std::same_as, + UniformComplexRandVar::sym>> +constexpr auto combine_magnitude(RandVarsView rand_vars) { + using sym = std::ranges::range_value_t::sym; + + assert(std::ranges::all_of(rand_vars, [](auto const& measurement) { + auto const check = [](auto const& value) { + return is_nan(imag(measurement.value)) || (real(measurement.value) > 0.0); + }; + if constexpr (is_symmetric_v) { + return check(value); + } else { + for (Idx const phase : IdxRange(3)) { + if (!check(measurement.value(phase))) { + return false; + } + } + return false; + } + })); + + auto const weighted_average_magnitude_measurement = + statistics::combine(rand_vars | std::views::transform([](auto const& measurement) { + return UniformRealRandVar{.value = detail::cabs_or_real(measurement.value), + .variance = measurement.variance}; + })); + return UniformComplexRandVar{.value = + [&weighted_average_magnitude_measurement]() { + ComplexValue abs_value = + piecewise_complex_value(DoubleComplex{0.0, nan}); + abs_value += weighted_average_magnitude_measurement.value; + return abs_value; + }(), + .variance = weighted_average_magnitude_measurement.variance}; +} } // namespace statistics } // namespace power_grid_model diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp index 12d5f73102..01704ec7ee 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp @@ -15,16 +15,6 @@ Collect all measured Values #include namespace power_grid_model::math_solver { - -namespace detail { -template inline RealValue cabs_or_real(ComplexValue const& value) { - if (is_nan(imag(value))) { - return real(value); // only keep real part - } - return cabs(value); // get abs of the value -} -} // namespace detail - // processed measurement struct // combined all measurement of the same quantity // accumulate for bus injection measurement @@ -455,19 +445,7 @@ template class MeasuredValues { IdxRange const& sensors) { auto complex_measurements = sensors | std::views::transform([&data](Idx pos) -> auto& { return data[pos]; }); if constexpr (only_magnitude) { - auto const weighted_average_magnitude_measurement = statistics::combine( - complex_measurements | std::views::transform([](auto const& measurement) { - return UniformRealRandVar{.value = detail::cabs_or_real(measurement.value), - .variance = measurement.variance}; - })); - return UniformComplexRandVar{.value = - [&weighted_average_magnitude_measurement]() { - ComplexValue abs_value = - piecewise_complex_value(DoubleComplex{0.0, nan}); - abs_value += weighted_average_magnitude_measurement.value; - return abs_value; - }(), - .variance = weighted_average_magnitude_measurement.variance}; + return statistics::combine_magnitude(complex_measurements); } else { return statistics::combine(complex_measurements); } diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp index faf69b16b3..5878bb3b17 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp @@ -212,6 +212,8 @@ template class NewtonRaphsonSESolver { typename SparseLUSolver, NRSERhs, NRSEUnknown>::BlockPermArray perm_; void initialize_unknown(ComplexValueVector& initial_u, MeasuredValues const& measured_values) { + using statistics::detail::cabs_or_real; + reset_unknown(); RealValue const mean_angle_shift = measured_values.mean_angle_shift(); for (Idx bus = 0; bus != n_bus_; ++bus) { @@ -222,7 +224,7 @@ template class NewtonRaphsonSESolver { if (measured_values.has_angle_measurement(bus)) { estimated_result.theta() = arg(measured_values.voltage(bus)); } - estimated_result.v() = detail::cabs_or_real(measured_values.voltage(bus)); + estimated_result.v() = cabs_or_real(measured_values.voltage(bus)); } initial_u[bus] = estimated_result.v() * exp(1.0i * estimated_result.theta()); } @@ -530,6 +532,8 @@ template class NewtonRaphsonSESolver { /// @param bus bus with voltage measurement void process_voltage_measurements(NRSEGainBlock& block, NRSERhs& rhs_block, MeasuredValues const& measured_values, Idx const& bus) { + using statistics::detail::cabs_or_real; + if (!measured_values.has_voltage(bus)) { return; } @@ -537,7 +541,7 @@ template class NewtonRaphsonSESolver { // G += 1.0 / variance // for 3x3 tensor, fill diagonal auto const w_v = RealTensor{1.0 / measured_values.voltage_var(bus)}; - auto const abs_measured_v = detail::cabs_or_real(measured_values.voltage(bus)); + auto const abs_measured_v = cabs_or_real(measured_values.voltage(bus)); auto const delta_v = abs_measured_v - x_[bus].v(); auto const virtual_angle_measurement_bus = measured_values.has_voltage(math_topo_->slack_bus) diff --git a/tests/cpp_unit_tests/test_statistics.cpp b/tests/cpp_unit_tests/test_statistics.cpp index a7ac5e1f2b..ce18de5c65 100644 --- a/tests/cpp_unit_tests/test_statistics.cpp +++ b/tests/cpp_unit_tests/test_statistics.cpp @@ -815,7 +815,7 @@ TEST_CASE("Test statistics") { } } -TEST_CASE("Test statistics - accumulate") { +TEST_CASE("Test statistics - combine") { using statistics::combine; using std::views::take; @@ -1129,6 +1129,81 @@ TEST_CASE("Test statistics - accumulate") { } } +TEST_CASE("Test statistics - combine_magnitude") { + using statistics::combine_magnitude; + using std::views::take; + + SUBCASE("UniformComplexRandVar") { + // using a template lambda to avoid code duplication and to avoid having to create a separate test case + std::vector> const measurements{ + {.value = ComplexValue{1.0, 5.0}, .variance = 0.2}, + {.value = ComplexValue{2.0, nan}, .variance = 0.3}, + {.value = ComplexValue{4.0, nan}, .variance = 0.6}}; + + CHECK(combine_magnitude(measurements | take(0)).value.real() == 0.0); + CHECK(is_nan(combine_magnitude(measurements | take(0)).value.imag())); + CHECK(is_inf(combine_magnitude(measurements | take(0)).variance)); + + CHECK(combine_magnitude(measurements | take(1)).value.real() == cabs(measurements.front().value)); + CHECK(is_nan(combine_magnitude(measurements | take(1)).value.imag())); + CHECK(combine_magnitude(measurements | take(1)).variance == measurements.front().variance); + + CHECK(combine_magnitude(measurements | take(2)).value.real() == + doctest::Approx((3.0 * std::sqrt(26.0) + 4.0) / 5.0)); + CHECK(is_nan(combine_magnitude(measurements | take(2)).value.imag())); + CHECK(combine_magnitude(measurements | take(2)).variance == doctest::Approx(3.0 / 25.0)); + + CHECK(combine_magnitude(measurements | take(3)).value.real() == + doctest::Approx((3.0 * std::sqrt(26.0) + 5.0) / 6.0)); + CHECK(is_nan(combine_magnitude(measurements | take(3)).value.imag())); + CHECK(combine_magnitude(measurements | take(3)).variance == doctest::Approx(1.0 / 10.0)); + } + + SUBCASE("UniformComplexRandVar") { + std::vector> const measurements{ + {.value = {RealValue{1.0, 2.0, -1.0}, RealValue{5.0, nan, 7.0}}, + .variance = 0.2}, + {.value = {RealValue{2.0, 4.0, 3.0}, RealValue{nan}}, .variance = 0.3}, + {.value = {RealValue{4.0, 5.0, 6.0}, RealValue{nan}}, .variance = 0.6}}; + + CHECK(combine_magnitude(measurements | take(0)).value(0).real() == 0.0); + CHECK(combine_magnitude(measurements | take(0)).value(1).real() == 0.0); + CHECK(combine_magnitude(measurements | take(0)).value(2).real() == 0.0); + CHECK(is_nan(combine_magnitude(measurements | take(0)).value(0).imag())); + CHECK(is_nan(combine_magnitude(measurements | take(0)).value(0).imag())); + CHECK(is_nan(combine_magnitude(measurements | take(0)).value(0).imag())); + CHECK(is_inf(combine_magnitude(measurements | take(0)).variance)); + + CHECK(combine_magnitude(measurements | take(1)).value(0).real() == cabs(measurements.front().value(0))); + CHECK(combine_magnitude(measurements | take(1)).value(1).real() == cabs(measurements.front().value(1))); + CHECK(combine_magnitude(measurements | take(1)).value(2).real() == cabs(measurements.front().value(2))); + CHECK(is_nan(combine_magnitude(measurements | take(1)).value(0).imag())); + CHECK(is_nan(combine_magnitude(measurements | take(1)).value(1).imag())); + CHECK(is_nan(combine_magnitude(measurements | take(1)).value(2).imag())); + CHECK(combine_magnitude(measurements | take(1)).variance == measurements.front().variance); + + CHECK(combine_magnitude(measurements | take(2)).value(0).real() == + doctest::Approx((3.0 * std::sqrt(26.0) + 4.0) / 5.0)); + CHECK(combine_magnitude(measurements | take(2)).value(1).real() == doctest::Approx(14.0 / 5.0)); + CHECK(combine_magnitude(measurements | take(2)).value(2).real() == + doctest::Approx((6.0 + 15.0 * std::sqrt(2.0)) / 5.0)); + CHECK(is_nan(combine_magnitude(measurements | take(2)).value(0).imag())); + CHECK(is_nan(combine_magnitude(measurements | take(2)).value(1).imag())); + CHECK(is_nan(combine_magnitude(measurements | take(2)).value(2).imag())); + CHECK(combine_magnitude(measurements | take(2)).variance == doctest::Approx(3.0 / 25.0)); + + CHECK(combine_magnitude(measurements | take(3)).value(0).real() == + doctest::Approx((3.0 * std::sqrt(26.0) + 5.0) / 6.0)); + CHECK(combine_magnitude(measurements | take(3)).value(1).real() == doctest::Approx(19.0 / 6.0)); + CHECK(combine_magnitude(measurements | take(3)).value(2).real() == + doctest::Approx((4.0 + 5.0 * std::sqrt(2.0)) / 2.0)); + CHECK(is_nan(combine_magnitude(measurements | take(3)).value(0).imag())); + CHECK(is_nan(combine_magnitude(measurements | take(3)).value(1).imag())); + CHECK(is_nan(combine_magnitude(measurements | take(3)).value(2).imag())); + CHECK(combine_magnitude(measurements | take(3)).variance == doctest::Approx(1.0 / 10.0)); + } +} + TEST_SUITE_END(); } // namespace power_grid_model From 51d0e8e4bbd3e6b579ce0cf7e053b1ca4855f75b Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Mon, 14 Apr 2025 12:56:58 +0200 Subject: [PATCH 24/25] fix tests Signed-off-by: Martijn Govers --- .../power_grid_model/common/statistics.hpp | 12 ++--------- tests/cpp_unit_tests/test_statistics.cpp | 21 ++++++++++++------- 2 files changed, 15 insertions(+), 18 deletions(-) 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 17e1678823..2c09d935f7 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 @@ -275,18 +275,10 @@ constexpr auto combine_magnitude(RandVarsView rand_vars) { using sym = std::ranges::range_value_t::sym; assert(std::ranges::all_of(rand_vars, [](auto const& measurement) { - auto const check = [](auto const& value) { - return is_nan(imag(measurement.value)) || (real(measurement.value) > 0.0); - }; if constexpr (is_symmetric_v) { - return check(value); + return (!is_nan(imag(measurement.value))) || real(measurement.value) >= 0.0; } else { - for (Idx const phase : IdxRange(3)) { - if (!check(measurement.value(phase))) { - return false; - } - } - return false; + return (!is_nan(imag(measurement.value))) || (real(measurement.value) >= 0.0).all(); } })); diff --git a/tests/cpp_unit_tests/test_statistics.cpp b/tests/cpp_unit_tests/test_statistics.cpp index ce18de5c65..4df368497b 100644 --- a/tests/cpp_unit_tests/test_statistics.cpp +++ b/tests/cpp_unit_tests/test_statistics.cpp @@ -1154,14 +1154,14 @@ TEST_CASE("Test statistics - combine_magnitude") { CHECK(combine_magnitude(measurements | take(2)).variance == doctest::Approx(3.0 / 25.0)); CHECK(combine_magnitude(measurements | take(3)).value.real() == - doctest::Approx((3.0 * std::sqrt(26.0) + 5.0) / 6.0)); + doctest::Approx((8.0 + 3.0 * std::sqrt(26.0)) / 6.0)); CHECK(is_nan(combine_magnitude(measurements | take(3)).value.imag())); CHECK(combine_magnitude(measurements | take(3)).variance == doctest::Approx(1.0 / 10.0)); } SUBCASE("UniformComplexRandVar") { std::vector> const measurements{ - {.value = {RealValue{1.0, 2.0, -1.0}, RealValue{5.0, nan, 7.0}}, + {.value = {RealValue{1.0, 2.0, -1.0}, RealValue{5.0, 6.0, 7.0}}, .variance = 0.2}, {.value = {RealValue{2.0, 4.0, 3.0}, RealValue{nan}}, .variance = 0.3}, {.value = {RealValue{4.0, 5.0, 6.0}, RealValue{nan}}, .variance = 0.6}}; @@ -1174,9 +1174,12 @@ TEST_CASE("Test statistics - combine_magnitude") { CHECK(is_nan(combine_magnitude(measurements | take(0)).value(0).imag())); CHECK(is_inf(combine_magnitude(measurements | take(0)).variance)); - CHECK(combine_magnitude(measurements | take(1)).value(0).real() == cabs(measurements.front().value(0))); - CHECK(combine_magnitude(measurements | take(1)).value(1).real() == cabs(measurements.front().value(1))); - CHECK(combine_magnitude(measurements | take(1)).value(2).real() == cabs(measurements.front().value(2))); + CHECK(combine_magnitude(measurements | take(1)).value(0).real() == + doctest::Approx(cabs(measurements.front().value(0)))); + CHECK(combine_magnitude(measurements | take(1)).value(1).real() == + doctest::Approx(cabs(measurements.front().value(1)))); + CHECK(combine_magnitude(measurements | take(1)).value(2).real() == + doctest::Approx(cabs(measurements.front().value(2)))); CHECK(is_nan(combine_magnitude(measurements | take(1)).value(0).imag())); CHECK(is_nan(combine_magnitude(measurements | take(1)).value(1).imag())); CHECK(is_nan(combine_magnitude(measurements | take(1)).value(2).imag())); @@ -1184,7 +1187,8 @@ TEST_CASE("Test statistics - combine_magnitude") { CHECK(combine_magnitude(measurements | take(2)).value(0).real() == doctest::Approx((3.0 * std::sqrt(26.0) + 4.0) / 5.0)); - CHECK(combine_magnitude(measurements | take(2)).value(1).real() == doctest::Approx(14.0 / 5.0)); + CHECK(combine_magnitude(measurements | take(2)).value(1).real() == + doctest::Approx((8.0 + 6.0 * std::sqrt(10.0)) / 5.0)); CHECK(combine_magnitude(measurements | take(2)).value(2).real() == doctest::Approx((6.0 + 15.0 * std::sqrt(2.0)) / 5.0)); CHECK(is_nan(combine_magnitude(measurements | take(2)).value(0).imag())); @@ -1193,8 +1197,9 @@ TEST_CASE("Test statistics - combine_magnitude") { CHECK(combine_magnitude(measurements | take(2)).variance == doctest::Approx(3.0 / 25.0)); CHECK(combine_magnitude(measurements | take(3)).value(0).real() == - doctest::Approx((3.0 * std::sqrt(26.0) + 5.0) / 6.0)); - CHECK(combine_magnitude(measurements | take(3)).value(1).real() == doctest::Approx(19.0 / 6.0)); + doctest::Approx((8.0 + 3.0 * std::sqrt(26.0)) / 6.0)); + CHECK(combine_magnitude(measurements | take(3)).value(1).real() == + doctest::Approx((13.0 + 6.0 * std::sqrt(10.0)) / 6.0)); CHECK(combine_magnitude(measurements | take(3)).value(2).real() == doctest::Approx((4.0 + 5.0 * std::sqrt(2.0)) / 2.0)); CHECK(is_nan(combine_magnitude(measurements | take(3)).value(0).imag())); From 5b6113fe0bc49828120a932de9a6f0f0840ff5ed Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Mon, 14 Apr 2025 14:08:04 +0200 Subject: [PATCH 25/25] assertion in combine magnitude causes crashes in first batch validation trial run Signed-off-by: Martijn Govers --- .../include/power_grid_model/common/statistics.hpp | 8 -------- 1 file changed, 8 deletions(-) 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 2c09d935f7..247a9b6764 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 @@ -274,14 +274,6 @@ template constexpr auto combine_magnitude(RandVarsView rand_vars) { using sym = std::ranges::range_value_t::sym; - assert(std::ranges::all_of(rand_vars, [](auto const& measurement) { - if constexpr (is_symmetric_v) { - return (!is_nan(imag(measurement.value))) || real(measurement.value) >= 0.0; - } else { - return (!is_nan(imag(measurement.value))) || (real(measurement.value) >= 0.0).all(); - } - })); - auto const weighted_average_magnitude_measurement = statistics::combine(rand_vars | std::views::transform([](auto const& measurement) { return UniformRealRandVar{.value = detail::cabs_or_real(measurement.value),