Skip to content

Commit 4ab0e9c

Browse files
mgoversmgovers
authored
Merge pull request #954 from PowerGridModel/feature/current-sensor-math-solver
current sensor: math solver implementation
2 parents 5d38475 + e02333d commit 4ab0e9c

File tree

5 files changed

+138
-23
lines changed

5 files changed

+138
-23
lines changed

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

Lines changed: 56 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -66,6 +66,8 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
6666
using sym = sym_type;
6767

6868
static constexpr auto is_iterative = true;
69+
static constexpr auto has_current_sensor_implemented =
70+
true; // TODO(mgovers): for testing purposes; remove after NRSE has current sensor implemented
6971

7072
private:
7173
// block size 2 for symmetric, 6 for asym
@@ -146,6 +148,10 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
146148
&MeasuredValues<sym>::has_branch_to_power};
147149
static constexpr std::array branch_power_{&MeasuredValues<sym>::branch_from_power,
148150
&MeasuredValues<sym>::branch_to_power};
151+
static constexpr std::array has_branch_current_{&MeasuredValues<sym>::has_branch_from_current,
152+
&MeasuredValues<sym>::has_branch_to_current};
153+
static constexpr std::array branch_current_{&MeasuredValues<sym>::branch_from_current,
154+
&MeasuredValues<sym>::branch_to_current};
149155

150156
Idx n_bus_;
151157
// shared topo data
@@ -205,19 +211,32 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
205211
}
206212
// branch
207213
else {
208-
// branch from- and to-side index at 0, and 1 position
209-
IntS const b0 = static_cast<IntS>(type) / 2;
210-
IntS const b1 = static_cast<IntS>(type) % 2;
214+
auto const add_branch_measurement = [&block, &param, obj,
215+
type](IntS measured_side,
216+
RealValue<sym> const& branch_current_variance) {
217+
// branch from- and to-side index at 0, and 1 position
218+
IntS const b0 = static_cast<IntS>(type) / 2;
219+
IntS const b1 = static_cast<IntS>(type) % 2;
220+
block.g() += dot(hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b0]),
221+
diagonal_inverse(branch_current_variance),
222+
param.branch_param[obj].value[measured_side * 2 + b1]);
223+
};
211224
// measured at from-side: 0, to-side: 1
212225
for (IntS const measured_side : std::array<IntS, 2>{0, 1}) {
213226
// has measurement
214227
if (std::invoke(has_branch_power_[measured_side], measured_value, obj)) {
215228
// G += Y{side, b0}^H * (variance^-1) * Y{side, b1}
216229
auto const& power = std::invoke(branch_power_[measured_side], measured_value, obj);
217-
block.g() +=
218-
dot(hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b0]),
219-
diagonal_inverse(static_cast<IndependentComplexRandVar<sym>>(power).variance),
220-
param.branch_param[obj].value[measured_side * 2 + b1]);
230+
// assume voltage ~ 1 p.u. => current variance = power variance * 1² = power variance
231+
add_branch_measurement(measured_side,
232+
static_cast<IndependentComplexRandVar<sym>>(power).variance);
233+
}
234+
if (std::invoke(has_branch_current_[measured_side], measured_value, obj)) {
235+
// G += (variance^-1)
236+
auto const& current =
237+
std::invoke(branch_current_[measured_side], measured_value, obj).measurement;
238+
add_branch_measurement(measured_side,
239+
static_cast<IndependentComplexRandVar<sym>>(current).variance);
221240
}
222241
}
223242
}
@@ -296,23 +315,42 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
296315
}
297316
// branch
298317
else {
299-
// branch is either ff or tt
300-
IntS const b = static_cast<IntS>(type) / 2;
301-
assert(b == static_cast<IntS>(type) % 2);
318+
auto const add_branch_measurement = [&rhs_block, &param, obj,
319+
type](IntS measured_side,
320+
IndependentComplexRandVar<sym> const& current) {
321+
// branch is either ff or tt
322+
IntS const b = static_cast<IntS>(type) / 2;
323+
// eta += Y{side, b}^H * (variance^-1) * i_branch_{f, t}
324+
rhs_block.eta() +=
325+
dot(hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b]),
326+
327+
diagonal_inverse(current.variance), current.value);
328+
};
302329
// measured at from-side: 0, to-side: 1
303330
for (IntS const measured_side : std::array<IntS, 2>{0, 1}) {
331+
// the current needs to be calculated with the voltage of the measured bus side
332+
// NOTE: not the bus that is currently being processed!
333+
Idx const measured_bus = branch_bus_idx[obj][measured_side];
334+
304335
// has measurement
305336
if (std::invoke(has_branch_power_[measured_side], measured_value, obj)) {
306337
PowerSensorCalcParam<sym> const& m =
307338
std::invoke(branch_power_[measured_side], measured_value, obj);
308-
// the current needs to be calculated with the voltage of the measured bus side
309-
// NOTE: not the current bus!
310-
Idx const measured_bus = branch_bus_idx[obj][measured_side];
311-
// eta += Y{side, b}^H * (variance^-1) * i_branch_{f, t}
312-
rhs_block.eta() +=
313-
dot(hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b]),
314-
diagonal_inverse(static_cast<IndependentComplexRandVar<sym>>(m).variance),
315-
conj(m.value() / u[measured_bus]));
339+
auto const apparent_power = static_cast<IndependentComplexRandVar<sym>>(m);
340+
add_branch_measurement(measured_side,
341+
{.value = conj(apparent_power.value / u[measured_bus]),
342+
.variance = apparent_power.variance});
343+
}
344+
if (std::invoke(has_branch_current_[measured_side], measured_value, obj)) {
345+
CurrentSensorCalcParam<sym> const m =
346+
std::invoke(branch_current_[measured_side], measured_value, obj);
347+
// for local angle current sensors, the current needs to be offset with the phase offset
348+
// of the measured bus side. NOTE: not the bus that is currently being processed!
349+
auto measured_current = static_cast<IndependentComplexRandVar<sym>>(m.measurement);
350+
if (m.angle_measurement_type == AngleMeasurementType::local_angle) {
351+
measured_current.value *= phase_shift(u[measured_bus]); // offset with the phase shift
352+
}
353+
add_branch_measurement(measured_side, measured_current);
316354
}
317355
}
318356
}

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

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -91,6 +91,8 @@ template <symmetry_tag sym_type> class NewtonRaphsonSESolver {
9191
using sym = sym_type;
9292

9393
static constexpr auto is_iterative = true;
94+
static constexpr auto has_current_sensor_implemented =
95+
false; // TODO(mgovers): for testing purposes; remove after NRSE has current sensor implemented
9496

9597
private:
9698
enum class Order : IntS { row_major = 0, column_major = 1 };

‎tests/cpp_unit_tests/test_math_solver_se.hpp

Lines changed: 74 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -344,6 +344,80 @@ TEST_CASE_TEMPLATE_DEFINE("Test math solver - SE, measurements", SolverType, tes
344344
CHECK(real(output.branch[0].s_f) == doctest::Approx(1.95));
345345
}
346346

347+
SUBCASE("Source and branch current") {
348+
/*
349+
network, v means voltage measured, c means current measured
350+
351+
bus_0(v) -(c)-branch_0-- bus_1
352+
| |
353+
source_0(p) load_0
354+
355+
*/
356+
// to make the distinction between global and local angle current measurement
357+
auto const global_shift = phase_shift(1.0 + 1.0i);
358+
359+
topo.power_sensors_per_source = {from_sparse, {0, 1}};
360+
topo.current_sensors_per_branch_from = {from_sparse, {0, 1}};
361+
se_input.measured_voltage = {{.value = 1.0 * global_shift, .variance = 0.1}};
362+
363+
auto param_ptr = std::make_shared<MathModelParam<symmetric_t> const>(param);
364+
auto topo_ptr = std::make_shared<MathModelTopology const>(topo);
365+
YBus<symmetric_t> const y_bus_sym{topo_ptr, param_ptr};
366+
367+
SolverType solver{y_bus_sym, topo_ptr};
368+
369+
SUBCASE("Local angle current sensor") {
370+
se_input.measured_source_power = {{.real_component = {.value = 1.93, .variance = 0.05},
371+
.imag_component = {.value = 0.0, .variance = 0.05}}};
372+
se_input.measured_branch_from_current = {
373+
{.angle_measurement_type = AngleMeasurementType::local_angle,
374+
.measurement = {.real_component = {.value = 1.97, .variance = 0.05},
375+
.imag_component = {.value = 0.0, .variance = 0.05}}}};
376+
377+
output = run_state_estimation(solver, y_bus_sym, se_input, error_tolerance, num_iter, info);
378+
379+
if (SolverType::has_current_sensor_implemented) { // TODO(mgovers): for testing purposes; remove if
380+
// statement after NRSE has current sensor implemented
381+
CHECK(real(output.bus_injection[0]) == doctest::Approx(1.95));
382+
CHECK(real(output.source[0].s) == doctest::Approx(1.95));
383+
CHECK(real(output.branch[0].s_f) == doctest::Approx(1.95));
384+
CHECK(real(output.branch[0].i_f) == doctest::Approx(real(1.95 * global_shift)));
385+
CHECK(imag(output.branch[0].i_f) == doctest::Approx(imag(1.95 * global_shift)));
386+
} else {
387+
CHECK_FALSE(real(output.bus_injection[0]) == doctest::Approx(1.95));
388+
CHECK_FALSE(real(output.source[0].s) == doctest::Approx(1.95));
389+
CHECK_FALSE(real(output.branch[0].s_f) == doctest::Approx(1.95));
390+
CHECK_FALSE(real(output.branch[0].i_f) == doctest::Approx(real(1.95 * global_shift)));
391+
CHECK_FALSE(imag(output.branch[0].i_f) == doctest::Approx(imag(1.95 * global_shift)));
392+
}
393+
}
394+
SUBCASE("Global angle current sensor") {
395+
se_input.measured_source_power = {{.real_component = {.value = 1.93, .variance = 0.05},
396+
.imag_component = {.value = 0.0, .variance = 0.05}}};
397+
se_input.measured_branch_from_current = {
398+
{.angle_measurement_type = AngleMeasurementType::global_angle,
399+
.measurement = {.real_component = {.value = real(1.97 * global_shift), .variance = 0.05},
400+
.imag_component = {.value = imag(1.97 * global_shift), .variance = 0.05}}}};
401+
402+
output = run_state_estimation(solver, y_bus_sym, se_input, error_tolerance, num_iter, info);
403+
404+
if (SolverType::has_current_sensor_implemented) { // TODO(mgovers): for testing purposes; remove if
405+
// statement after NRSE has current sensor implemented
406+
CHECK(real(output.bus_injection[0]) == doctest::Approx(1.95));
407+
CHECK(real(output.source[0].s) == doctest::Approx(1.95));
408+
CHECK(real(output.branch[0].s_f) == doctest::Approx(1.95));
409+
CHECK(real(output.branch[0].i_f) == doctest::Approx(real(1.95 * global_shift)));
410+
CHECK(imag(output.branch[0].i_f) == doctest::Approx(imag(1.95 * global_shift)));
411+
} else {
412+
CHECK_FALSE(real(output.bus_injection[0]) == doctest::Approx(1.95));
413+
CHECK_FALSE(real(output.source[0].s) == doctest::Approx(1.95));
414+
CHECK_FALSE(real(output.branch[0].s_f) == doctest::Approx(1.95));
415+
CHECK_FALSE(real(output.branch[0].i_f) == doctest::Approx(real(1.95 * global_shift)));
416+
CHECK_FALSE(imag(output.branch[0].i_f) == doctest::Approx(imag(1.95 * global_shift)));
417+
}
418+
}
419+
}
420+
347421
SUBCASE("Load and branch") {
348422
/*
349423
network, v means voltage measured, p means power measured

‎tests/cpp_validation_tests/test_validation.cpp

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -295,6 +295,8 @@ std::map<std::string, PGM_TapChangingStrategy, std::less<>> const optimizer_stra
295295
{"min_voltage_tap", PGM_tap_changing_strategy_min_voltage_tap},
296296
{"max_voltage_tap", PGM_tap_changing_strategy_max_voltage_tap},
297297
{"fast_any_tap", PGM_tap_changing_strategy_fast_any_tap}};
298+
std::map<std::string, PGM_ExperimentalFeatures, std::less<>> const experimental_features_mapping = {
299+
{"disabled", PGM_experimental_features_disabled}, {"enabled", PGM_experimental_features_enabled}};
298300

299301
// case parameters
300302
struct CaseParam {
@@ -304,6 +306,7 @@ struct CaseParam {
304306
std::string calculation_method;
305307
std::string short_circuit_voltage_scaling;
306308
std::string tap_changing_strategy;
309+
std::string experimental_features;
307310
double err_tol = 1e-8;
308311
Idx max_iter = 20;
309312
bool sym{};
@@ -329,6 +332,7 @@ Options get_options(CaseParam const& param, Idx threading = -1) {
329332
options.set_threading(threading);
330333
options.set_short_circuit_voltage_scaling(sc_voltage_scaling_mapping.at(param.short_circuit_voltage_scaling));
331334
options.set_tap_changing_strategy(optimizer_strategy_mapping.at(param.tap_changing_strategy));
335+
options.set_experimental_features(experimental_features_mapping.at(param.experimental_features));
332336
return options;
333337
}
334338

@@ -390,6 +394,7 @@ std::optional<CaseParam> construct_case(std::filesystem::path const& case_dir, j
390394
}
391395

392396
param.tap_changing_strategy = calculation_method_params.value("tap_changing_strategy", "disabled");
397+
param.experimental_features = calculation_method_params.value("experimental_features", "disabled");
393398
param.case_name += sym ? "-sym"s : "-asym"s;
394399
param.case_name += "-"s + param.calculation_method;
395400
param.case_name += is_batch ? "_batch"s : ""s;

‎tests/data/state_estimation/basic-current-sensor/params.json

Lines changed: 1 addition & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -10,11 +10,7 @@
1010
},
1111
"extra_params": {
1212
"iterative_linear": {
13-
"experimental_features": "enabled",
14-
"fail": {
15-
"raises": "SparseMatrixError",
16-
"reason": "Current sensors are not yet implemented for this calculation method"
17-
}
13+
"experimental_features": "enabled"
1814
},
1915
"newton_raphson": {
2016
"experimental_features": "enabled",

0 commit comments

Comments
 (0)