Skip to content

Commit 04e1273

Browse files
mgoversmgovers
committed
also migrate combine_magnitude to statistics
Signed-off-by: Martijn Govers <Martijn.Govers@Alliander.com>
1 parent 1953c85 commit 04e1273

File tree

4 files changed

+129
-26
lines changed

4 files changed

+129
-26
lines changed

‎power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp

Lines changed: 46 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -258,5 +258,51 @@ constexpr auto combine(RandVarsView rand_vars) {
258258
}
259259
return result;
260260
}
261+
262+
namespace detail {
263+
template <symmetry_tag sym> inline RealValue<sym> cabs_or_real(ComplexValue<sym> const& value) {
264+
if (is_nan(imag(value))) {
265+
return real(value); // only keep real part
266+
}
267+
return cabs(value); // get abs of the value
268+
}
269+
} // namespace detail
270+
271+
template <std::ranges::view RandVarsView>
272+
requires std::same_as<std::ranges::range_value_t<RandVarsView>,
273+
UniformComplexRandVar<typename std::ranges::range_value_t<RandVarsView>::sym>>
274+
constexpr auto combine_magnitude(RandVarsView rand_vars) {
275+
using sym = std::ranges::range_value_t<RandVarsView>::sym;
276+
277+
assert(std::ranges::all_of(rand_vars, [](auto const& measurement) {
278+
auto const check = [](auto const& value) {
279+
return is_nan(imag(measurement.value)) || (real(measurement.value) > 0.0);
280+
};
281+
if constexpr (is_symmetric_v<sym>) {
282+
return check(value);
283+
} else {
284+
for (Idx const phase : IdxRange(3)) {
285+
if (!check(measurement.value(phase))) {
286+
return false;
287+
}
288+
}
289+
return false;
290+
}
291+
}));
292+
293+
auto const weighted_average_magnitude_measurement =
294+
statistics::combine(rand_vars | std::views::transform([](auto const& measurement) {
295+
return UniformRealRandVar<sym>{.value = detail::cabs_or_real<sym>(measurement.value),
296+
.variance = measurement.variance};
297+
}));
298+
return UniformComplexRandVar<sym>{.value =
299+
[&weighted_average_magnitude_measurement]() {
300+
ComplexValue<sym> abs_value =
301+
piecewise_complex_value<sym>(DoubleComplex{0.0, nan});
302+
abs_value += weighted_average_magnitude_measurement.value;
303+
return abs_value;
304+
}(),
305+
.variance = weighted_average_magnitude_measurement.variance};
306+
}
261307
} // namespace statistics
262308
} // namespace power_grid_model

‎power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp

Lines changed: 1 addition & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -15,16 +15,6 @@ Collect all measured Values
1515
#include <memory>
1616

1717
namespace power_grid_model::math_solver {
18-
19-
namespace detail {
20-
template <symmetry_tag sym> inline RealValue<sym> cabs_or_real(ComplexValue<sym> const& value) {
21-
if (is_nan(imag(value))) {
22-
return real(value); // only keep real part
23-
}
24-
return cabs(value); // get abs of the value
25-
}
26-
} // namespace detail
27-
2818
// processed measurement struct
2919
// combined all measurement of the same quantity
3020
// accumulate for bus injection measurement
@@ -455,19 +445,7 @@ template <symmetry_tag sym> class MeasuredValues {
455445
IdxRange const& sensors) {
456446
auto complex_measurements = sensors | std::views::transform([&data](Idx pos) -> auto& { return data[pos]; });
457447
if constexpr (only_magnitude) {
458-
auto const weighted_average_magnitude_measurement = statistics::combine(
459-
complex_measurements | std::views::transform([](auto const& measurement) {
460-
return UniformRealRandVar<sym>{.value = detail::cabs_or_real<sym>(measurement.value),
461-
.variance = measurement.variance};
462-
}));
463-
return UniformComplexRandVar<sym>{.value =
464-
[&weighted_average_magnitude_measurement]() {
465-
ComplexValue<sym> abs_value =
466-
piecewise_complex_value<sym>(DoubleComplex{0.0, nan});
467-
abs_value += weighted_average_magnitude_measurement.value;
468-
return abs_value;
469-
}(),
470-
.variance = weighted_average_magnitude_measurement.variance};
448+
return statistics::combine_magnitude(complex_measurements);
471449
} else {
472450
return statistics::combine(complex_measurements);
473451
}

‎power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -212,6 +212,8 @@ template <symmetry_tag sym_type> class NewtonRaphsonSESolver {
212212
typename SparseLUSolver<NRSEGainBlock<sym>, NRSERhs<sym>, NRSEUnknown<sym>>::BlockPermArray perm_;
213213

214214
void initialize_unknown(ComplexValueVector<sym>& initial_u, MeasuredValues<sym> const& measured_values) {
215+
using statistics::detail::cabs_or_real;
216+
215217
reset_unknown();
216218
RealValue<sym> const mean_angle_shift = measured_values.mean_angle_shift();
217219
for (Idx bus = 0; bus != n_bus_; ++bus) {
@@ -222,7 +224,7 @@ template <symmetry_tag sym_type> class NewtonRaphsonSESolver {
222224
if (measured_values.has_angle_measurement(bus)) {
223225
estimated_result.theta() = arg(measured_values.voltage(bus));
224226
}
225-
estimated_result.v() = detail::cabs_or_real<sym>(measured_values.voltage(bus));
227+
estimated_result.v() = cabs_or_real<sym>(measured_values.voltage(bus));
226228
}
227229
initial_u[bus] = estimated_result.v() * exp(1.0i * estimated_result.theta());
228230
}
@@ -530,14 +532,16 @@ template <symmetry_tag sym_type> class NewtonRaphsonSESolver {
530532
/// @param bus bus with voltage measurement
531533
void process_voltage_measurements(NRSEGainBlock<sym>& block, NRSERhs<sym>& rhs_block,
532534
MeasuredValues<sym> const& measured_values, Idx const& bus) {
535+
using statistics::detail::cabs_or_real;
536+
533537
if (!measured_values.has_voltage(bus)) {
534538
return;
535539
}
536540

537541
// G += 1.0 / variance
538542
// for 3x3 tensor, fill diagonal
539543
auto const w_v = RealTensor<sym>{1.0 / measured_values.voltage_var(bus)};
540-
auto const abs_measured_v = detail::cabs_or_real<sym>(measured_values.voltage(bus));
544+
auto const abs_measured_v = cabs_or_real<sym>(measured_values.voltage(bus));
541545
auto const delta_v = abs_measured_v - x_[bus].v();
542546

543547
auto const virtual_angle_measurement_bus = measured_values.has_voltage(math_topo_->slack_bus)

‎tests/cpp_unit_tests/test_statistics.cpp

Lines changed: 76 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -815,7 +815,7 @@ TEST_CASE("Test statistics") {
815815
}
816816
}
817817

818-
TEST_CASE("Test statistics - accumulate") {
818+
TEST_CASE("Test statistics - combine") {
819819
using statistics::combine;
820820
using std::views::take;
821821

@@ -1129,6 +1129,81 @@ TEST_CASE("Test statistics - accumulate") {
11291129
}
11301130
}
11311131

1132+
TEST_CASE("Test statistics - combine_magnitude") {
1133+
using statistics::combine_magnitude;
1134+
using std::views::take;
1135+
1136+
SUBCASE("UniformComplexRandVar<symmetric_t>") {
1137+
// using a template lambda to avoid code duplication and to avoid having to create a separate test case
1138+
std::vector<UniformComplexRandVar<symmetric_t>> const measurements{
1139+
{.value = ComplexValue<symmetric_t>{1.0, 5.0}, .variance = 0.2},
1140+
{.value = ComplexValue<symmetric_t>{2.0, nan}, .variance = 0.3},
1141+
{.value = ComplexValue<symmetric_t>{4.0, nan}, .variance = 0.6}};
1142+
1143+
CHECK(combine_magnitude(measurements | take(0)).value.real() == 0.0);
1144+
CHECK(is_nan(combine_magnitude(measurements | take(0)).value.imag()));
1145+
CHECK(is_inf(combine_magnitude(measurements | take(0)).variance));
1146+
1147+
CHECK(combine_magnitude(measurements | take(1)).value.real() == cabs(measurements.front().value));
1148+
CHECK(is_nan(combine_magnitude(measurements | take(1)).value.imag()));
1149+
CHECK(combine_magnitude(measurements | take(1)).variance == measurements.front().variance);
1150+
1151+
CHECK(combine_magnitude(measurements | take(2)).value.real() ==
1152+
doctest::Approx((3.0 * std::sqrt(26.0) + 4.0) / 5.0));
1153+
CHECK(is_nan(combine_magnitude(measurements | take(2)).value.imag()));
1154+
CHECK(combine_magnitude(measurements | take(2)).variance == doctest::Approx(3.0 / 25.0));
1155+
1156+
CHECK(combine_magnitude(measurements | take(3)).value.real() ==
1157+
doctest::Approx((3.0 * std::sqrt(26.0) + 5.0) / 6.0));
1158+
CHECK(is_nan(combine_magnitude(measurements | take(3)).value.imag()));
1159+
CHECK(combine_magnitude(measurements | take(3)).variance == doctest::Approx(1.0 / 10.0));
1160+
}
1161+
1162+
SUBCASE("UniformComplexRandVar<asymmetric_t>") {
1163+
std::vector<UniformComplexRandVar<asymmetric_t>> const measurements{
1164+
{.value = {RealValue<asymmetric_t>{1.0, 2.0, -1.0}, RealValue<asymmetric_t>{5.0, nan, 7.0}},
1165+
.variance = 0.2},
1166+
{.value = {RealValue<asymmetric_t>{2.0, 4.0, 3.0}, RealValue<asymmetric_t>{nan}}, .variance = 0.3},
1167+
{.value = {RealValue<asymmetric_t>{4.0, 5.0, 6.0}, RealValue<asymmetric_t>{nan}}, .variance = 0.6}};
1168+
1169+
CHECK(combine_magnitude(measurements | take(0)).value(0).real() == 0.0);
1170+
CHECK(combine_magnitude(measurements | take(0)).value(1).real() == 0.0);
1171+
CHECK(combine_magnitude(measurements | take(0)).value(2).real() == 0.0);
1172+
CHECK(is_nan(combine_magnitude(measurements | take(0)).value(0).imag()));
1173+
CHECK(is_nan(combine_magnitude(measurements | take(0)).value(0).imag()));
1174+
CHECK(is_nan(combine_magnitude(measurements | take(0)).value(0).imag()));
1175+
CHECK(is_inf(combine_magnitude(measurements | take(0)).variance));
1176+
1177+
CHECK(combine_magnitude(measurements | take(1)).value(0).real() == cabs(measurements.front().value(0)));
1178+
CHECK(combine_magnitude(measurements | take(1)).value(1).real() == cabs(measurements.front().value(1)));
1179+
CHECK(combine_magnitude(measurements | take(1)).value(2).real() == cabs(measurements.front().value(2)));
1180+
CHECK(is_nan(combine_magnitude(measurements | take(1)).value(0).imag()));
1181+
CHECK(is_nan(combine_magnitude(measurements | take(1)).value(1).imag()));
1182+
CHECK(is_nan(combine_magnitude(measurements | take(1)).value(2).imag()));
1183+
CHECK(combine_magnitude(measurements | take(1)).variance == measurements.front().variance);
1184+
1185+
CHECK(combine_magnitude(measurements | take(2)).value(0).real() ==
1186+
doctest::Approx((3.0 * std::sqrt(26.0) + 4.0) / 5.0));
1187+
CHECK(combine_magnitude(measurements | take(2)).value(1).real() == doctest::Approx(14.0 / 5.0));
1188+
CHECK(combine_magnitude(measurements | take(2)).value(2).real() ==
1189+
doctest::Approx((6.0 + 15.0 * std::sqrt(2.0)) / 5.0));
1190+
CHECK(is_nan(combine_magnitude(measurements | take(2)).value(0).imag()));
1191+
CHECK(is_nan(combine_magnitude(measurements | take(2)).value(1).imag()));
1192+
CHECK(is_nan(combine_magnitude(measurements | take(2)).value(2).imag()));
1193+
CHECK(combine_magnitude(measurements | take(2)).variance == doctest::Approx(3.0 / 25.0));
1194+
1195+
CHECK(combine_magnitude(measurements | take(3)).value(0).real() ==
1196+
doctest::Approx((3.0 * std::sqrt(26.0) + 5.0) / 6.0));
1197+
CHECK(combine_magnitude(measurements | take(3)).value(1).real() == doctest::Approx(19.0 / 6.0));
1198+
CHECK(combine_magnitude(measurements | take(3)).value(2).real() ==
1199+
doctest::Approx((4.0 + 5.0 * std::sqrt(2.0)) / 2.0));
1200+
CHECK(is_nan(combine_magnitude(measurements | take(3)).value(0).imag()));
1201+
CHECK(is_nan(combine_magnitude(measurements | take(3)).value(1).imag()));
1202+
CHECK(is_nan(combine_magnitude(measurements | take(3)).value(2).imag()));
1203+
CHECK(combine_magnitude(measurements | take(3)).variance == doctest::Approx(1.0 / 10.0));
1204+
}
1205+
}
1206+
11321207
TEST_SUITE_END();
11331208

11341209
} // namespace power_grid_model

0 commit comments

Comments
 (0)