Skip to content

Current sensor: minor refactor ilse math solver #958

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 4 commits into from
Apr 16, 2025
Merged

Current sensor: minor refactor ilse math solver #958

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
bd448eb
add conjugation to statistics
mgovers Apr 16, 2025
d66e476
add conjugation to statistics
mgovers Apr 16, 2025
7d04eaf
more clearly reflect maths behind ILSE
mgovers Apr 16, 2025
841a8ef
update readme to get rid of duplicate CI badges
mgovers Apr 16, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 1 addition & 4 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
18 changes: 18 additions & 0 deletions power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -207,6 +207,24 @@ template <symmetry_tag sym_type> struct PolarComplexRandVar {
}
};

template <symmetry_tag sym, template <symmetry_tag> typename RandVarType>
requires std::same_as<RandVarType<sym>, UniformComplexRandVar<sym>> ||
std::same_as<RandVarType<sym>, IndependentComplexRandVar<sym>>
inline auto conj(RandVarType<sym> var) {
var.value = conj(var.value);
return var;
}

template <symmetry_tag sym> inline auto conj(DecomposedComplexRandVar<sym> var) {
var.imag_component.value = -var.imag_component.value;
return var;
}

template <symmetry_tag sym> inline auto conj(PolarComplexRandVar<sym> var) {
var.angle.value = -var.angle.value;
return var;
}

namespace statistics {
// Var(s x) ≈ Var(x) * ||s||²
template <symmetry_tag sym, template <symmetry_tag> typename RandVarType>
Expand Down
125 changes: 83 additions & 42 deletions ...el_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -203,40 +203,39 @@ template <symmetry_tag sym_type> 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<IndependentComplexRandVar<sym>>(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, &param, obj,
type](IntS measured_side,
RealValue<sym> const& branch_current_variance) {
auto const add_branch_measurement = [&block, &param, obj, type](
IntS measured_side,
IndependentComplexRandVar<sym> const& branch_current) {
// branch from- and to-side index at 0, and 1 position
IntS const b0 = static_cast<IntS>(type) / 2;
IntS const b1 = static_cast<IntS>(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<IntS, 2>{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<IndependentComplexRandVar<sym>>(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<IndependentComplexRandVar<sym>>(current).variance);
current_to_global_current_measurement(branch_current));
}
}
}
Expand All @@ -250,8 +249,8 @@ template <symmetry_tag sym_type> 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<sym>{static_cast<ComplexValue<sym>>(
-(static_cast<IndependentComplexRandVar<sym>>(injection).variance))};
block.r() = ComplexTensor<sym>{
static_cast<ComplexValue<sym>>(-power_to_global_current_measurement(injection).variance)};
}
}
// injection measurement not exist
Expand Down Expand Up @@ -305,26 +304,24 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
// shunt
if (type == YBusElementType::shunt) {
if (measured_value.has_shunt(obj)) {
PowerSensorCalcParam<sym> const& m = measured_value.shunt_power(obj);
PowerSensorCalcParam<sym> 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<IndependentComplexRandVar<sym>>(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, &param, obj,
type](IntS measured_side,
IndependentComplexRandVar<sym> const& current) {
IndependentComplexRandVar<sym> const& branch_current) {
// branch is either ff or tt
IntS const b = static_cast<IntS>(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<IntS, 2>{0, 1}) {
Expand All @@ -334,31 +331,22 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {

// has measurement
if (std::invoke(has_branch_power_[measured_side], measured_value, obj)) {
PowerSensorCalcParam<sym> const& m =
std::invoke(branch_power_[measured_side], measured_value, obj);
auto const apparent_power = static_cast<IndependentComplexRandVar<sym>>(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<sym> 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<IndependentComplexRandVar<sym>>(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;
}
}
}
Expand Down Expand Up @@ -394,9 +382,62 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
return max_dev;
}

auto linearize_measurements(ComplexValueVector<sym> const& current_u, MeasuredValues<sym> const& measured_values) {
auto linearize_measurements(ComplexValueVector<sym> const& current_u,
MeasuredValues<sym> 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<sym>
power_to_global_current_measurement(PowerSensorCalcParam<sym> const& power_measurement,
ComplexValue<sym> const& voltage) const {
auto measurement = static_cast<IndependentComplexRandVar<sym>>(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<sym>
power_to_global_current_measurement(PowerSensorCalcParam<sym> const& power_measurement) const {
auto measurement = static_cast<IndependentComplexRandVar<sym>>(power_measurement);
measurement.value = ComplexValue<sym>{nan};
return measurement;
}

IndependentComplexRandVar<sym>
current_to_global_current_measurement(CurrentSensorCalcParam<sym> const& current_measurement,
ComplexValue<sym> const& voltage) const {
using statistics::scale;

auto const measurement = static_cast<IndependentComplexRandVar<sym>>(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<sym>
current_to_global_current_measurement(CurrentSensorCalcParam<sym> const& current_measurement) const {
auto measurement = static_cast<IndependentComplexRandVar<sym>>(current_measurement.measurement);
measurement.value = ComplexValue<sym>{nan};
return measurement;
}
};

} // namespace iterative_linear_se
Expand Down
Loading
Loading