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.

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
Show file tree
Hide file tree
Changes from all commits
Commits
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
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
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