From bd448ebbebefc981deaba72e9494d823daff019c Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Wed, 16 Apr 2025 09:17:31 +0200 Subject: [PATCH 1/4] add conjugation to statistics Signed-off-by: Martijn Govers --- .../power_grid_model/common/statistics.hpp | 17 ++++ tests/cpp_unit_tests/test_statistics.cpp | 98 +++++++++++++++++++ 2 files changed, 115 insertions(+) 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 3686156d2e..7d03d9249a 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 @@ -206,6 +206,23 @@ template struct PolarComplexRandVar { 9.0}}; } }; +template typename RandVarType> + requires std::same_as, UniformComplexRandVar> || + std::same_as, IndependentComplexRandVar> +inline auto conj(RandVarType var) { + var.value = conj(var.value); + return var; +} + +template inline auto conj(DecomposedComplexRandVar var) { + var.imag_component.value *= -1; + return var; +} + +template inline auto conj(PolarComplexRandVar var) { + var.angle.value = -var.angle.value; + return var; +} namespace statistics { // Var(s x) ≈ Var(x) * ||s||² diff --git a/tests/cpp_unit_tests/test_statistics.cpp b/tests/cpp_unit_tests/test_statistics.cpp index 504433959a..46ac8652a8 100644 --- a/tests/cpp_unit_tests/test_statistics.cpp +++ b/tests/cpp_unit_tests/test_statistics.cpp @@ -815,6 +815,104 @@ TEST_CASE("Test statistics") { } } +TEST_CASE("Test statistics - conjugate") { + SUBCASE("UniformComplexRandVar | IndependentComplexRandVar") { + auto const check = [](auto const& var) { + auto const conjugated = conj(var); + CHECK(real(conjugated.value) == real(conj(var.value))); + CHECK(imag(conjugated.value) == imag(conj(var.value))); + CHECK(conjugated.variance == var.variance); + }; + + SUBCASE("UniformComplexRandVar") { + check(UniformComplexRandVar{.value = {1.0, 2.0}, .variance = 2.0}); + } + SUBCASE("IndependentComplexRandVar") { + check(IndependentComplexRandVar{.value = {1.0, 2.0}, .variance = 2.0}); + } + } + SUBCASE("UniformComplexRandVar") { + auto const var = UniformComplexRandVar{ + .value = {RealValue{1.0, 2.0, 3.0}, RealValue{4.0, 5.0, 6.0}}, .variance = 2.0}; + auto const conjugated = conj(var); + CHECK(real(conjugated.value(0)) == real(conj(var.value(0)))); + CHECK(imag(conjugated.value(0)) == imag(conj(var.value(0)))); + CHECK(real(conjugated.value(1)) == real(conj(var.value(1)))); + CHECK(imag(conjugated.value(1)) == imag(conj(var.value(1)))); + CHECK(real(conjugated.value(2)) == real(conj(var.value(2)))); + CHECK(imag(conjugated.value(2)) == imag(conj(var.value(2)))); + CHECK(conjugated.variance == var.variance); + } + SUBCASE("IndependentComplexRandVar") { + auto const var = IndependentComplexRandVar{ + .value = {RealValue{1.0, 2.0, 3.0}, RealValue{4.0, 5.0, 6.0}}, + .variance = {2.0, 3.0, 4.0}}; + auto const conjugated = conj(var); + CHECK(real(conjugated.value(0)) == real(conj(var.value(0)))); + CHECK(imag(conjugated.value(0)) == imag(conj(var.value(0)))); + CHECK(real(conjugated.value(1)) == real(conj(var.value(1)))); + CHECK(imag(conjugated.value(1)) == imag(conj(var.value(1)))); + CHECK(real(conjugated.value(2)) == real(conj(var.value(2)))); + CHECK(imag(conjugated.value(2)) == imag(conj(var.value(2)))); + CHECK(conjugated.variance(0) == var.variance(0)); + CHECK(conjugated.variance(1) == var.variance(1)); + CHECK(conjugated.variance(2) == var.variance(2)); + } + SUBCASE("DecomposedComplexRandVar") { + auto const var = DecomposedComplexRandVar{.real_component = {.value = 1.0, .variance = 2.0}, + .imag_component = {.value = 3.0, .variance = 4.0}}; + auto const conjugated = conj(var); + CHECK(real(conjugated.value()) == real(conj(var.value()))); + CHECK(imag(conjugated.value()) == imag(conj(var.value()))); + CHECK(conjugated.real_component.variance == var.real_component.variance); + CHECK(conjugated.imag_component.variance == var.imag_component.variance); + } + SUBCASE("DecomposedComplexRandVar") { + auto const var = DecomposedComplexRandVar{ + .real_component = {.value = RealValue{1.0, 2.0, 3.0}, + .variance = RealValue{2.0, 3.0, 4.0}}, + .imag_component = {.value = RealValue{4.0, 5.0, 6.0}, + .variance = RealValue{3.0, 4.0, 5.0}}}; + auto const conjugated = conj(var); + CHECK(real(conjugated.value()(0)) == real(conj(var.value())(0))); + CHECK(imag(conjugated.value()(0)) == imag(conj(var.value())(0))); + CHECK(real(conjugated.value()(1)) == real(conj(var.value())(1))); + CHECK(imag(conjugated.value()(1)) == imag(conj(var.value())(1))); + CHECK(real(conjugated.value()(2)) == real(conj(var.value())(2))); + CHECK(imag(conjugated.value()(2)) == imag(conj(var.value())(2))); + CHECK(conjugated.real_component.variance(0) == var.real_component.variance(0)); + CHECK(conjugated.imag_component.variance(0) == var.imag_component.variance(0)); + CHECK(conjugated.real_component.variance(1) == var.real_component.variance(1)); + CHECK(conjugated.imag_component.variance(1) == var.imag_component.variance(1)); + CHECK(conjugated.real_component.variance(2) == var.real_component.variance(2)); + CHECK(conjugated.imag_component.variance(2) == var.imag_component.variance(2)); + } + + SUBCASE("PolarComplexRandVar") { + auto const var = PolarComplexRandVar{.magnitude = {.value = 1.0, .variance = 2.0}, + .angle = {.value = 3.0, .variance = 4.0}}; + auto const conjugated = conj(var); + CHECK(real(conjugated.value()) == real(conj(var.value()))); + CHECK(imag(conjugated.value()) == imag(conj(var.value()))); + CHECK(conjugated.magnitude.variance == var.magnitude.variance); + CHECK(conjugated.angle.variance == var.angle.variance); + } + SUBCASE("PolarComplexRandVar") { + auto const var = PolarComplexRandVar{ + .magnitude = {.value = RealValue{1.0, 2.0, 3.0}, .variance = 2.0}, + .angle = {.value = RealValue{4.0, 5.0, 6.0}, .variance = 4.0}}; + auto const conjugated = conj(var); + CHECK(real(conjugated.value()(0)) == real(conj(var.value()(0)))); + CHECK(imag(conjugated.value()(0)) == imag(conj(var.value()(0)))); + CHECK(real(conjugated.value()(1)) == real(conj(var.value()(1)))); + CHECK(imag(conjugated.value()(1)) == imag(conj(var.value()(1)))); + CHECK(real(conjugated.value()(2)) == real(conj(var.value()(2)))); + CHECK(imag(conjugated.value()(2)) == imag(conj(var.value()(2)))); + CHECK(conjugated.magnitude.variance == var.magnitude.variance); + CHECK(conjugated.angle.variance == var.angle.variance); + } +} + TEST_CASE("Test statistics - scale") { using statistics::scale; From d66e476003a7e8b49e4df4538f8c3aba1c5c215f Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Wed, 16 Apr 2025 09:24:53 +0200 Subject: [PATCH 2/4] add conjugation to statistics Signed-off-by: Martijn Govers --- .../include/power_grid_model/common/statistics.hpp | 2 +- 1 file changed, 1 insertion(+), 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 7d03d9249a..2760947e3c 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 @@ -215,7 +215,7 @@ inline auto conj(RandVarType var) { } template inline auto conj(DecomposedComplexRandVar var) { - var.imag_component.value *= -1; + var.imag_component.value = -var.imag_component.value; return var; } From 7d04eaf9d91e46a6ca9f5fc1552344e3df1098b7 Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Wed, 16 Apr 2025 11:55:27 +0200 Subject: [PATCH 3/4] more clearly reflect maths behind ILSE Signed-off-by: Martijn Govers --- .../power_grid_model/common/statistics.hpp | 1 + .../iterative_linear_se_solver.hpp | 125 ++++++++++++------ 2 files changed, 84 insertions(+), 42 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 2760947e3c..c8c01ad500 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 @@ -206,6 +206,7 @@ template struct PolarComplexRandVar { 9.0}}; } }; + template typename RandVarType> requires std::same_as, UniformComplexRandVar> || std::same_as, IndependentComplexRandVar> 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 7ab95e281a..1a8806cda8 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 @@ -203,40 +203,39 @@ template class IterativeLinearSESolver { if (measured_value.has_shunt(obj)) { // G += (-Ys)^H * (variance^-1) * (-Ys) auto const& shunt_power = measured_value.shunt_power(obj); - block.g() += - dot(hermitian_transpose(param.shunt_param[obj]), - diagonal_inverse(static_cast>(shunt_power).variance), - param.shunt_param[obj]); + auto const current = power_to_global_current_measurement(shunt_power); + block.g() += dot(hermitian_transpose(param.shunt_param[obj]), + diagonal_inverse(current.variance), param.shunt_param[obj]); } } // branch else { - auto const add_branch_measurement = [&block, ¶m, obj, - type](IntS measured_side, - RealValue const& branch_current_variance) { + auto const add_branch_measurement = [&block, ¶m, obj, type]( + IntS measured_side, + IndependentComplexRandVar const& branch_current) { // branch from- and to-side index at 0, and 1 position IntS const b0 = static_cast(type) / 2; IntS const b1 = static_cast(type) % 2; + + // G += Y{side, b0}^H * (variance^-1) * Y{side, b1} block.g() += dot(hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b0]), - diagonal_inverse(branch_current_variance), + diagonal_inverse(branch_current.variance), param.branch_param[obj].value[measured_side * 2 + b1]); }; // measured at from-side: 0, to-side: 1 for (IntS const measured_side : std::array{0, 1}) { // has measurement 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); - // assume voltage ~ 1 p.u. => current variance = power variance * 1² = power variance + auto const& branch_power = + std::invoke(branch_power_[measured_side], measured_value, obj); add_branch_measurement(measured_side, - static_cast>(power).variance); + power_to_global_current_measurement(branch_power)); } if (std::invoke(has_branch_current_[measured_side], measured_value, obj)) { - // G += (variance^-1) - auto const& current = - std::invoke(branch_current_[measured_side], measured_value, obj).measurement; + auto const& branch_current = + std::invoke(branch_current_[measured_side], measured_value, obj); add_branch_measurement(measured_side, - static_cast>(current).variance); + current_to_global_current_measurement(branch_current)); } } } @@ -250,8 +249,8 @@ template class IterativeLinearSESolver { if (row == col) { // assign variance to diagonal of 3x3 tensor, for asym auto const& injection = measured_value.bus_injection(row); - block.r() = ComplexTensor{static_cast>( - -(static_cast>(injection).variance))}; + block.r() = ComplexTensor{ + static_cast>(-power_to_global_current_measurement(injection).variance)}; } } // injection measurement not exist @@ -305,26 +304,24 @@ template class IterativeLinearSESolver { // shunt if (type == YBusElementType::shunt) { if (measured_value.has_shunt(obj)) { - PowerSensorCalcParam const& m = measured_value.shunt_power(obj); + PowerSensorCalcParam const& shunt_power = measured_value.shunt_power(obj); + auto const current = power_to_global_current_measurement(shunt_power, u[bus]); // eta += (-Ys)^H * (variance^-1) * i_shunt - rhs_block.eta() -= - dot(hermitian_transpose(param.shunt_param[obj]), - diagonal_inverse(static_cast>(m).variance), - conj(m.value() / u[bus])); + rhs_block.eta() -= dot(hermitian_transpose(param.shunt_param[obj]), + diagonal_inverse(current.variance), current.value); } } // branch else { auto const add_branch_measurement = [&rhs_block, ¶m, obj, type](IntS measured_side, - IndependentComplexRandVar const& current) { + IndependentComplexRandVar const& branch_current) { // branch is either ff or tt IntS const b = static_cast(type) / 2; // eta += Y{side, b}^H * (variance^-1) * i_branch_{f, t} rhs_block.eta() += dot(hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b]), - - diagonal_inverse(current.variance), current.value); + diagonal_inverse(branch_current.variance), branch_current.value); }; // measured at from-side: 0, to-side: 1 for (IntS const measured_side : std::array{0, 1}) { @@ -334,31 +331,22 @@ template class IterativeLinearSESolver { // has measurement if (std::invoke(has_branch_power_[measured_side], measured_value, obj)) { - PowerSensorCalcParam const& m = - std::invoke(branch_power_[measured_side], measured_value, obj); - auto const apparent_power = static_cast>(m); + auto const& branch_power = std::invoke(branch_power_[measured_side], measured_value, obj); add_branch_measurement(measured_side, - {.value = conj(apparent_power.value / u[measured_bus]), - .variance = apparent_power.variance}); + power_to_global_current_measurement(branch_power, u[measured_bus])); } if (std::invoke(has_branch_current_[measured_side], measured_value, obj)) { - CurrentSensorCalcParam const m = + auto const& branch_current = std::invoke(branch_current_[measured_side], measured_value, obj); - // for local angle current sensors, the current needs to be offset with the phase offset - // of the measured bus side. NOTE: not the bus that is currently being processed! - auto measured_current = static_cast>(m.measurement); - if (m.angle_measurement_type == AngleMeasurementType::local_angle) { - measured_current.value = conj(measured_current.value) * - phase_shift(u[measured_bus]); // offset with the phase shift - } - add_branch_measurement(measured_side, measured_current); + add_branch_measurement( + measured_side, current_to_global_current_measurement(branch_current, u[measured_bus])); } } } } // fill block with injection measurement, need to convert to current if (measured_value.has_bus_injection(bus)) { - rhs_block.tau() = conj(measured_value.bus_injection(bus).value() / u[bus]); + rhs_block.tau() = power_to_global_current_measurement(measured_value.bus_injection(bus), u[bus]).value; } } } @@ -394,9 +382,62 @@ template class IterativeLinearSESolver { return max_dev; } - auto linearize_measurements(ComplexValueVector const& current_u, MeasuredValues const& measured_values) { + auto linearize_measurements(ComplexValueVector const& current_u, + MeasuredValues const& measured_values) const { return measured_values.combine_voltage_iteration_with_measurements(current_u); } + + // The variance is not scaled as an approximation under the assumptions of: + // - linearization of obtaining the current from power measurements + // - voltages are ~1pu + // - power sensor variances are often an approximation dominated by heuristics in the first place + // See also https://github.com/PowerGridModel/power-grid-model/pull/951#issuecomment-2805154436 + IndependentComplexRandVar + power_to_global_current_measurement(PowerSensorCalcParam const& power_measurement, + ComplexValue const& voltage) const { + auto measurement = static_cast>(power_measurement); + measurement.value = conj(measurement.value / voltage); + return measurement; + } + + // Overload when the voltage is not present: the value can't be determined, but the variance assumption still holds. + // The variance is not scaled as an approximation under the assumptions of: + // - linearization of obtaining the current from power measurements + // - voltages are ~1pu + // - power sensor variances are often an approximation dominated by heuristics in the first place + // See also https://github.com/PowerGridModel/power-grid-model/pull/951#issuecomment-2805154436 + IndependentComplexRandVar + power_to_global_current_measurement(PowerSensorCalcParam const& power_measurement) const { + auto measurement = static_cast>(power_measurement); + measurement.value = ComplexValue{nan}; + return measurement; + } + + IndependentComplexRandVar + current_to_global_current_measurement(CurrentSensorCalcParam const& current_measurement, + ComplexValue const& voltage) const { + using statistics::scale; + + auto const measurement = static_cast>(current_measurement.measurement); + + switch (current_measurement.angle_measurement_type) { + case AngleMeasurementType::global_angle: + return measurement; // no offset + case AngleMeasurementType::local_angle: + return scale(conj(measurement), + phase_shift(voltage)); // offset with the phase shift + default: + throw MissingCaseForEnumError{"AngleMeasurementType", current_measurement.angle_measurement_type}; + } + } + + // Overload when the voltage is not present: the value can't be determined, but the variance assumption still holds. + IndependentComplexRandVar + current_to_global_current_measurement(CurrentSensorCalcParam const& current_measurement) const { + auto measurement = static_cast>(current_measurement.measurement); + measurement.value = ComplexValue{nan}; + return measurement; + } }; } // namespace iterative_linear_se From 841a8effd5717aab1fd0629640bfd60a1918dc52 Mon Sep 17 00:00:00 2001 From: Martijn Govers Date: Wed, 16 Apr 2025 12:46:26 +0200 Subject: [PATCH 4/4] update readme to get rid of duplicate CI badges Signed-off-by: Martijn Govers --- README.md | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/README.md b/README.md index 8cc086a85a..7f5674b2cb 100644 --- a/README.md +++ b/README.md @@ -10,10 +10,7 @@ SPDX-License-Identifier: MPL-2.0 [![Downloads](https://static.pepy.tech/badge/power-grid-model)](https://pepy.tech/project/power-grid-model) [![Downloads](https://static.pepy.tech/badge/power-grid-model/month)](https://pepy.tech/project/power-grid-model) -[![Build and Test C++ and Python](https://github.com/PowerGridModel/power-grid-model/actions/workflows/ci.yml/badge.svg)](https://github.com/PowerGridModel/power-grid-model/actions/workflows/build-test-release.yml) -[![Check Code Quality](https://github.com/PowerGridModel/power-grid-model/actions/workflows/ci.yml/badge.svg)](https://github.com/PowerGridModel/power-grid-model/actions/workflows/check-code-quality.yml) -[![Clang Tidy](https://github.com/PowerGridModel/power-grid-model/actions/workflows/ci.yml/badge.svg)](https://github.com/PowerGridModel/power-grid-model/actions/workflows/clang-tidy.yml) -[![REUSE Compliance Check](https://github.com/PowerGridModel/power-grid-model/actions/workflows/ci.yml/badge.svg)](https://github.com/PowerGridModel/power-grid-model/actions/workflows/reuse-compliance.yml) +[![CI Build](https://github.com/PowerGridModel/power-grid-model/actions/workflows/ci.yml/badge.svg)](https://github.com/PowerGridModel/power-grid-model/actions/workflows/ci.yml) [![docs](https://readthedocs.org/projects/power-grid-model/badge/)](https://power-grid-model.readthedocs.io/en/stable/) [![Quality Gate Status](https://sonarcloud.io/api/project_badges/measure?project=PowerGridModel_power-grid-model&metric=alert_status)](https://sonarcloud.io/summary/new_code?id=PowerGridModel_power-grid-model)