Skip to content

Iterative Linear State Estimation: power sensor accuracy #951

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

Closed
wants to merge 9 commits into from
Closed

Iterative Linear State Estimation: power sensor accuracy #951

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
c153797
first implementation of ILSE with rescaled variances
mgovers Apr 9, 2025
f57ccb3
fix eigen operations
mgovers Apr 10, 2025
076c6f8
update tests to reflect migration to scaled power sensor variance
mgovers Apr 10, 2025
e976694
minor additional requirement
mgovers Apr 10, 2025
f621ab2
make clang happy
mgovers Apr 10, 2025
83b4c48
improve naming
mgovers Apr 11, 2025
dfd8d17
Merge branch 'main' into feature/ilse-power-sensor-accuracy
mgovers Apr 14, 2025
a71a4ee
add scaling functionality to statistics.hpp
mgovers Apr 14, 2025
a9109dd
Merge branch 'feature/add-scaling-to-statistics' into feature/ilse-po…
mgovers Apr 15, 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
45 changes: 45 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 @@ -208,6 +208,51 @@ template <symmetry_tag sym_type> struct PolarComplexRandVar {
};

namespace statistics {
template <symmetry_tag sym, template <symmetry_tag> typename RandVarType>
requires is_in_list_c<RandVarType<sym>, UniformRealRandVar<sym>, IndependentRealRandVar<sym>,
UniformComplexRandVar<sym>, IndependentComplexRandVar<sym>>
inline auto scale(RandVarType<sym> const& var, double scale_factor) {
return RandVarType<sym>{.value = var.value * scale_factor, .variance = var.variance * abs2(scale_factor)};
}

template <typename RandVarType>
requires is_in_list_c<RandVarType, IndependentRealRandVar<asymmetric_t>, IndependentComplexRandVar<asymmetric_t>>
inline auto scale(RandVarType const& var, RealValue<asymmetric_t> const& scale_factor) {
return RandVarType{.value = var.value * scale_factor, .variance = var.variance * abs2(scale_factor)};
}

template <symmetry_tag sym, template <symmetry_tag> typename RandVarType>
requires is_in_list_c<RandVarType<sym>, UniformComplexRandVar<sym>, IndependentComplexRandVar<sym>>
inline auto scale(RandVarType<sym> const& var, DoubleComplex const& scale_factor) {
return RandVarType<sym>{.value = var.value * scale_factor, .variance = var.variance * abs2(scale_factor)};
}

inline auto scale(IndependentComplexRandVar<asymmetric_t> const& var, ComplexValue<asymmetric_t> const& scale_factor) {
return IndependentComplexRandVar<asymmetric_t>{.value = var.value * scale_factor,
.variance = var.variance * abs2(scale_factor)};
}

template <symmetry_tag sym, typename ScaleType>
requires is_in_list_c<ScaleType, RealValue<symmetric_t>, RealValue<asymmetric_t>>
inline auto scale(DecomposedComplexRandVar<sym> const& var, ScaleType const& scale_factor) {
return DecomposedComplexRandVar<sym>{.real_component = scale(var.real_component, scale_factor),
.imag_component = scale(var.imag_component, scale_factor)};
}

template <symmetry_tag sym, typename ScaleType>
requires(std::same_as<ScaleType, ComplexValue<symmetric_t>> ||
(is_asymmetric_v<sym> && std::same_as<ScaleType, ComplexValue<asymmetric_t>>))
inline auto scale(DecomposedComplexRandVar<sym> const& var, ScaleType const& scale_factor) {
ComplexValue<sym> const scaled_value = var.value() * scale_factor;
return DecomposedComplexRandVar<sym>{
.real_component = {.value = real(scaled_value),
.variance = var.real_component.variance * abs2(real(scale_factor)) +
var.imag_component.variance * abs2(imag(scale_factor))},
.imag_component = {.value = imag(scaled_value),
.variance = var.real_component.variance * abs2(imag(scale_factor)) +
var.imag_component.variance * abs2(real(scale_factor))}};
}

// combine multiple random variables of one quantity using Kalman filter
template <std::ranges::view RandVarsView>
requires std::same_as<std::ranges::range_value_t<RandVarsView>,
Expand Down
74 changes: 49 additions & 25 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 @@ -98,12 +98,6 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
auto const observability_result =
necessary_observability_check(measured_values, y_bus.math_topology(), y_bus.y_bus_structure());

// prepare matrix
sub_timer = Timer(calculation_info, 2222, "Prepare matrix, including pre-factorization");
prepare_matrix(y_bus, measured_values);
// prefactorize
sparse_solver_.prefactorize(data_gain_, perm_, observability_result.use_perturbation());

// initialize voltage with initial angle
sub_timer = Timer(calculation_info, 2223, "Initialize voltages");
RealValue<sym> const mean_angle_shift = measured_values.mean_angle_shift();
Expand All @@ -117,8 +111,21 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
if (num_iter++ == max_iter) {
throw IterationDiverge{max_iter, max_dev, err_tol};
}

// prepare matrix
sub_timer = Timer(calculation_info, 2222, "Linearize measurements");
// get generated (measured/estimated) voltage phasor
// with current result voltage angle
ComplexValueVector<sym> u = linearize_measurements(output.u, measured_values);

// prepare matrix
sub_timer = Timer(calculation_info, 2222, "Prepare matrix, including pre-factorization");
prepare_matrix(y_bus, measured_values, u);
// prefactorize
sparse_solver_.prefactorize(data_gain_, perm_, observability_result.use_perturbation());

sub_timer = Timer(calculation_info, 2224, "Calculate rhs");
prepare_rhs(y_bus, measured_values, output.u);
prepare_rhs(y_bus, measured_values, u);
// solve with prefactorization
sub_timer = Timer(calculation_info, 2225, "Solve sparse linear equation (pre-factorized)");
sparse_solver_.solve_with_prefactorized_matrix(data_gain_, perm_, x_rhs_, x_rhs_);
Expand Down Expand Up @@ -163,10 +170,12 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
return ComplexDiagonalTensor<sym>{static_cast<ComplexValue<sym>>(RealValue<sym>{1.0} / value)};
}

void prepare_matrix(YBus<sym> const& y_bus, MeasuredValues<sym> const& measured_value) {
void prepare_matrix(YBus<sym> const& y_bus, MeasuredValues<sym> const& measured_value,
ComplexValueVector<sym> const& u) {
MathModelParam<sym> const& param = y_bus.math_model_param();
IdxVector const& row_indptr = y_bus.row_indptr_lu();
IdxVector const& col_indices = y_bus.col_indices_lu();
std::vector<BranchIdx> const& branch_bus_idx = y_bus.math_topology().branch_bus_idx;

// loop data index, all rows and columns
for (Idx row = 0; row != n_bus_; ++row) {
Expand Down Expand Up @@ -197,9 +206,12 @@ 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);
auto const shunt_current_conj =
statistics::scale(shunt_power, ComplexValue<sym>{1.0 / u[row]});
block.g() +=
dot(hermitian_transpose(param.shunt_param[obj]),
diagonal_inverse(static_cast<IndependentComplexRandVar<sym>>(shunt_power).variance),
diagonal_inverse(
static_cast<IndependentComplexRandVar<sym>>(shunt_current_conj).variance),
param.shunt_param[obj]);
}
}
Expand All @@ -214,9 +226,15 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
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);
// the current needs to be calculated with the voltage of the measured bus side
// NOTE: not the current bus!
Idx const measured_bus = branch_bus_idx[obj][measured_side];
auto const current_conj =
statistics::scale(power, ComplexValue<sym>{1.0 / u[measured_bus]});
block.g() +=
dot(hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b0]),
diagonal_inverse(static_cast<IndependentComplexRandVar<sym>>(power).variance),
diagonal_inverse(
static_cast<IndependentComplexRandVar<sym>>(current_conj).variance),
param.branch_param[obj].value[measured_side * 2 + b1]);
}
}
Expand All @@ -232,7 +250,9 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
// 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))};
-(static_cast<IndependentComplexRandVar<sym>>(
statistics::scale(injection, ComplexValue<sym>{1.0 / u[row]}))
.variance))};
}
}
// injection measurement not exist
Expand Down Expand Up @@ -260,12 +280,9 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
}

void prepare_rhs(YBus<sym> const& y_bus, MeasuredValues<sym> const& measured_value,
ComplexValueVector<sym> const& current_u) {
ComplexValueVector<sym> const& u) {
MathModelParam<sym> const& param = y_bus.math_model_param();
std::vector<BranchIdx> const branch_bus_idx = y_bus.math_topology().branch_bus_idx;
// get generated (measured/estimated) voltage phasor
// with current result voltage angle
ComplexValueVector<sym> u = linearize_measurements(current_u, measured_value);
std::vector<BranchIdx> const& branch_bus_idx = y_bus.math_topology().branch_bus_idx;

// loop all bus to fill rhs
for (Idx bus = 0; bus != n_bus_; ++bus) {
Expand All @@ -286,12 +303,14 @@ 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& power = measured_value.shunt_power(obj);
auto const current_conj = statistics::scale(power, ComplexValue<sym>{1.0 / u[bus]});
auto const rot_invariant_current_conj =
static_cast<IndependentComplexRandVar<sym>>(current_conj);
// 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(rot_invariant_current_conj.variance),
conj(rot_invariant_current_conj.value));
}
}
// branch
Expand All @@ -303,23 +322,28 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
for (IntS const measured_side : std::array<IntS, 2>{0, 1}) {
// has measurement
if (std::invoke(has_branch_power_[measured_side], measured_value, obj)) {
PowerSensorCalcParam<sym> const& m =
PowerSensorCalcParam<sym> const& power =
std::invoke(branch_power_[measured_side], measured_value, obj);
// the current needs to be calculated with the voltage of the measured bus side
// NOTE: not the current bus!
Idx const measured_bus = branch_bus_idx[obj][measured_side];
auto const current_conj =
statistics::scale(power, ComplexValue<sym>{1.0 / u[measured_bus]});
auto const rot_invariant_current_conj =
static_cast<IndependentComplexRandVar<sym>>(current_conj);
// 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(static_cast<IndependentComplexRandVar<sym>>(m).variance),
conj(m.value() / u[measured_bus]));
diagonal_inverse(rot_invariant_current_conj.variance),
conj(rot_invariant_current_conj.value));
}
}
}
}
// 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() =
conj(statistics::scale(measured_value.bus_injection(bus), ComplexValue<sym>{1.0 / u[bus]}).value());
}
}
}
Expand Down
Loading
Loading