Skip to content

Commit 33c4176

Browse files
TonyXiang8787TonyXiang8787
authored
Merge pull request #958 from PowerGridModel/feature/minor-refactor-ilse
Current sensor: minor refactor ilse math solver
2 parents c7e3730 + 841a8ef commit 33c4176

File tree

4 files changed

+200
-46
lines changed

4 files changed

+200
-46
lines changed

‎README.md

Lines changed: 1 addition & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -10,10 +10,7 @@ SPDX-License-Identifier: MPL-2.0
1010
[![Downloads](https://static.pepy.tech/badge/power-grid-model)](https://pepy.tech/project/power-grid-model)
1111
[![Downloads](https://static.pepy.tech/badge/power-grid-model/month)](https://pepy.tech/project/power-grid-model)
1212

13-
[![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)
14-
[![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)
15-
[![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)
16-
[![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)
13+
[![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)
1714
[![docs](https://readthedocs.org/projects/power-grid-model/badge/)](https://power-grid-model.readthedocs.io/en/stable/)
1815

1916
[![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)

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

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -207,6 +207,24 @@ template <symmetry_tag sym_type> struct PolarComplexRandVar {
207207
}
208208
};
209209

210+
template <symmetry_tag sym, template <symmetry_tag> typename RandVarType>
211+
requires std::same_as<RandVarType<sym>, UniformComplexRandVar<sym>> ||
212+
std::same_as<RandVarType<sym>, IndependentComplexRandVar<sym>>
213+
inline auto conj(RandVarType<sym> var) {
214+
var.value = conj(var.value);
215+
return var;
216+
}
217+
218+
template <symmetry_tag sym> inline auto conj(DecomposedComplexRandVar<sym> var) {
219+
var.imag_component.value = -var.imag_component.value;
220+
return var;
221+
}
222+
223+
template <symmetry_tag sym> inline auto conj(PolarComplexRandVar<sym> var) {
224+
var.angle.value = -var.angle.value;
225+
return var;
226+
}
227+
210228
namespace statistics {
211229
// Var(s x) ≈ Var(x) * ||s||²
212230
template <symmetry_tag sym, template <symmetry_tag> typename RandVarType>

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

Lines changed: 83 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -203,40 +203,39 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
203203
if (measured_value.has_shunt(obj)) {
204204
// G += (-Ys)^H * (variance^-1) * (-Ys)
205205
auto const& shunt_power = measured_value.shunt_power(obj);
206-
block.g() +=
207-
dot(hermitian_transpose(param.shunt_param[obj]),
208-
diagonal_inverse(static_cast<IndependentComplexRandVar<sym>>(shunt_power).variance),
209-
param.shunt_param[obj]);
206+
auto const current = power_to_global_current_measurement(shunt_power);
207+
block.g() += dot(hermitian_transpose(param.shunt_param[obj]),
208+
diagonal_inverse(current.variance), param.shunt_param[obj]);
210209
}
211210
}
212211
// branch
213212
else {
214-
auto const add_branch_measurement = [&block, &param, obj,
215-
type](IntS measured_side,
216-
RealValue<sym> const& branch_current_variance) {
213+
auto const add_branch_measurement = [&block, &param, obj, type](
214+
IntS measured_side,
215+
IndependentComplexRandVar<sym> const& branch_current) {
217216
// branch from- and to-side index at 0, and 1 position
218217
IntS const b0 = static_cast<IntS>(type) / 2;
219218
IntS const b1 = static_cast<IntS>(type) % 2;
219+
220+
// G += Y{side, b0}^H * (variance^-1) * Y{side, b1}
220221
block.g() += dot(hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b0]),
221-
diagonal_inverse(branch_current_variance),
222+
diagonal_inverse(branch_current.variance),
222223
param.branch_param[obj].value[measured_side * 2 + b1]);
223224
};
224225
// measured at from-side: 0, to-side: 1
225226
for (IntS const measured_side : std::array<IntS, 2>{0, 1}) {
226227
// has measurement
227228
if (std::invoke(has_branch_power_[measured_side], measured_value, obj)) {
228-
// G += Y{side, b0}^H * (variance^-1) * Y{side, b1}
229-
auto const& power = std::invoke(branch_power_[measured_side], measured_value, obj);
230-
// assume voltage ~ 1 p.u. => current variance = power variance * 1² = power variance
229+
auto const& branch_power =
230+
std::invoke(branch_power_[measured_side], measured_value, obj);
231231
add_branch_measurement(measured_side,
232-
static_cast<IndependentComplexRandVar<sym>>(power).variance);
232+
power_to_global_current_measurement(branch_power));
233233
}
234234
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;
235+
auto const& branch_current =
236+
std::invoke(branch_current_[measured_side], measured_value, obj);
238237
add_branch_measurement(measured_side,
239-
static_cast<IndependentComplexRandVar<sym>>(current).variance);
238+
current_to_global_current_measurement(branch_current));
240239
}
241240
}
242241
}
@@ -250,8 +249,8 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
250249
if (row == col) {
251250
// assign variance to diagonal of 3x3 tensor, for asym
252251
auto const& injection = measured_value.bus_injection(row);
253-
block.r() = ComplexTensor<sym>{static_cast<ComplexValue<sym>>(
254-
-(static_cast<IndependentComplexRandVar<sym>>(injection).variance))};
252+
block.r() = ComplexTensor<sym>{
253+
static_cast<ComplexValue<sym>>(-power_to_global_current_measurement(injection).variance)};
255254
}
256255
}
257256
// injection measurement not exist
@@ -305,26 +304,24 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
305304
// shunt
306305
if (type == YBusElementType::shunt) {
307306
if (measured_value.has_shunt(obj)) {
308-
PowerSensorCalcParam<sym> const& m = measured_value.shunt_power(obj);
307+
PowerSensorCalcParam<sym> const& shunt_power = measured_value.shunt_power(obj);
308+
auto const current = power_to_global_current_measurement(shunt_power, u[bus]);
309309
// eta += (-Ys)^H * (variance^-1) * i_shunt
310-
rhs_block.eta() -=
311-
dot(hermitian_transpose(param.shunt_param[obj]),
312-
diagonal_inverse(static_cast<IndependentComplexRandVar<sym>>(m).variance),
313-
conj(m.value() / u[bus]));
310+
rhs_block.eta() -= dot(hermitian_transpose(param.shunt_param[obj]),
311+
diagonal_inverse(current.variance), current.value);
314312
}
315313
}
316314
// branch
317315
else {
318316
auto const add_branch_measurement = [&rhs_block, &param, obj,
319317
type](IntS measured_side,
320-
IndependentComplexRandVar<sym> const& current) {
318+
IndependentComplexRandVar<sym> const& branch_current) {
321319
// branch is either ff or tt
322320
IntS const b = static_cast<IntS>(type) / 2;
323321
// eta += Y{side, b}^H * (variance^-1) * i_branch_{f, t}
324322
rhs_block.eta() +=
325323
dot(hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b]),
326-
327-
diagonal_inverse(current.variance), current.value);
324+
diagonal_inverse(branch_current.variance), branch_current.value);
328325
};
329326
// measured at from-side: 0, to-side: 1
330327
for (IntS const measured_side : std::array<IntS, 2>{0, 1}) {
@@ -334,31 +331,22 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
334331

335332
// has measurement
336333
if (std::invoke(has_branch_power_[measured_side], measured_value, obj)) {
337-
PowerSensorCalcParam<sym> const& m =
338-
std::invoke(branch_power_[measured_side], measured_value, obj);
339-
auto const apparent_power = static_cast<IndependentComplexRandVar<sym>>(m);
334+
auto const& branch_power = std::invoke(branch_power_[measured_side], measured_value, obj);
340335
add_branch_measurement(measured_side,
341-
{.value = conj(apparent_power.value / u[measured_bus]),
342-
.variance = apparent_power.variance});
336+
power_to_global_current_measurement(branch_power, u[measured_bus]));
343337
}
344338
if (std::invoke(has_branch_current_[measured_side], measured_value, obj)) {
345-
CurrentSensorCalcParam<sym> const m =
339+
auto const& branch_current =
346340
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 = conj(measured_current.value) *
352-
phase_shift(u[measured_bus]); // offset with the phase shift
353-
}
354-
add_branch_measurement(measured_side, measured_current);
341+
add_branch_measurement(
342+
measured_side, current_to_global_current_measurement(branch_current, u[measured_bus]));
355343
}
356344
}
357345
}
358346
}
359347
// fill block with injection measurement, need to convert to current
360348
if (measured_value.has_bus_injection(bus)) {
361-
rhs_block.tau() = conj(measured_value.bus_injection(bus).value() / u[bus]);
349+
rhs_block.tau() = power_to_global_current_measurement(measured_value.bus_injection(bus), u[bus]).value;
362350
}
363351
}
364352
}
@@ -394,9 +382,62 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
394382
return max_dev;
395383
}
396384

397-
auto linearize_measurements(ComplexValueVector<sym> const& current_u, MeasuredValues<sym> const& measured_values) {
385+
auto linearize_measurements(ComplexValueVector<sym> const& current_u,
386+
MeasuredValues<sym> const& measured_values) const {
398387
return measured_values.combine_voltage_iteration_with_measurements(current_u);
399388
}
389+
390+
// The variance is not scaled as an approximation under the assumptions of:
391+
// - linearization of obtaining the current from power measurements
392+
// - voltages are ~1pu
393+
// - power sensor variances are often an approximation dominated by heuristics in the first place
394+
// See also https://github.com/PowerGridModel/power-grid-model/pull/951#issuecomment-2805154436
395+
IndependentComplexRandVar<sym>
396+
power_to_global_current_measurement(PowerSensorCalcParam<sym> const& power_measurement,
397+
ComplexValue<sym> const& voltage) const {
398+
auto measurement = static_cast<IndependentComplexRandVar<sym>>(power_measurement);
399+
measurement.value = conj(measurement.value / voltage);
400+
return measurement;
401+
}
402+
403+
// Overload when the voltage is not present: the value can't be determined, but the variance assumption still holds.
404+
// The variance is not scaled as an approximation under the assumptions of:
405+
// - linearization of obtaining the current from power measurements
406+
// - voltages are ~1pu
407+
// - power sensor variances are often an approximation dominated by heuristics in the first place
408+
// See also https://github.com/PowerGridModel/power-grid-model/pull/951#issuecomment-2805154436
409+
IndependentComplexRandVar<sym>
410+
power_to_global_current_measurement(PowerSensorCalcParam<sym> const& power_measurement) const {
411+
auto measurement = static_cast<IndependentComplexRandVar<sym>>(power_measurement);
412+
measurement.value = ComplexValue<sym>{nan};
413+
return measurement;
414+
}
415+
416+
IndependentComplexRandVar<sym>
417+
current_to_global_current_measurement(CurrentSensorCalcParam<sym> const& current_measurement,
418+
ComplexValue<sym> const& voltage) const {
419+
using statistics::scale;
420+
421+
auto const measurement = static_cast<IndependentComplexRandVar<sym>>(current_measurement.measurement);
422+
423+
switch (current_measurement.angle_measurement_type) {
424+
case AngleMeasurementType::global_angle:
425+
return measurement; // no offset
426+
case AngleMeasurementType::local_angle:
427+
return scale(conj(measurement),
428+
phase_shift(voltage)); // offset with the phase shift
429+
default:
430+
throw MissingCaseForEnumError{"AngleMeasurementType", current_measurement.angle_measurement_type};
431+
}
432+
}
433+
434+
// Overload when the voltage is not present: the value can't be determined, but the variance assumption still holds.
435+
IndependentComplexRandVar<sym>
436+
current_to_global_current_measurement(CurrentSensorCalcParam<sym> const& current_measurement) const {
437+
auto measurement = static_cast<IndependentComplexRandVar<sym>>(current_measurement.measurement);
438+
measurement.value = ComplexValue<sym>{nan};
439+
return measurement;
440+
}
400441
};
401442

402443
} // namespace iterative_linear_se

0 commit comments

Comments
 (0)