Skip to content

Commit

Permalink
Merge pull request #478 from PowerGridModel/fix/nrse-line-measurements
Browse files Browse the repository at this point in the history 
Fix NRSE line measurement implementation
  • Loading branch information
@mgovers
mgovers authored Feb 6, 2024
2 parents 22feb91 + cdcd535 commit 5eddc40
Show file tree
Hide file tree
Showing 18 changed files with 591 additions and 305 deletions.
18 changes: 7 additions & 11 deletions ...el_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,7 @@
#ifndef POWER_GRID_MODEL_MATH_SOLVER_ITERATIVE_LINEAR_SE_SOLVER_HPP
#define POWER_GRID_MODEL_MATH_SOLVER_ITERATIVE_LINEAR_SE_SOLVER_HPP

/*
iterative linear state estimation solver
*/
// iterative linear state estimation solver

#include "block_matrix.hpp"
#include "common_solver_functions.hpp"
Expand Down Expand Up @@ -46,12 +44,10 @@ template <bool sym> struct ILSEUnknown : public Block<DoubleComplex, sym, false,
template <bool sym> using ILSERhs = ILSEUnknown<sym>;

// class of 2*2 (6*6) se gain block
/*
[
[G, QH]
[Q, R ]
]
*/
// [
// [G, QH]
// [Q, R ]
// ]
template <bool sym> class ILSEGainBlock : public Block<DoubleComplex, sym, true, 2> {
public:
template <int r, int c> using GetterType = typename Block<DoubleComplex, sym, true, 2>::template GetterType<r, c>;
Expand Down Expand Up @@ -189,7 +185,7 @@ template <bool sym> class IterativeLinearSESolver {
// shunt
if (type == YBusElementType::shunt) {
if (measured_value.has_shunt(obj)) {
// G += Ys^H * (variance^-1) * Ys
// 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(shunt_power.p_variance + shunt_power.q_variance),
Expand Down Expand Up @@ -282,7 +278,7 @@ template <bool sym> class IterativeLinearSESolver {
if (type == YBusElementType::shunt) {
if (measured_value.has_shunt(obj)) {
PowerSensorCalcParam<sym> const& m = measured_value.shunt_power(obj);
// eta -= Ys^H * (variance^-1) * i_shunt
// eta += (-Ys)^H * (variance^-1) * i_shunt
rhs_block.eta() -= dot(hermitian_transpose(param.shunt_param[obj]),
diagonal_inverse(m.p_variance + m.q_variance), conj(m.value / u[bus]));
}
Expand Down
6 changes: 3 additions & 3 deletions power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -242,7 +242,7 @@ template <bool sym> class MeasuredValues {

// assign a meaningful mean angle shift, if at least one voltage has angle measurement
if (has_angle()) {
mean_angle_shift_ = angle_cum / n_voltage_angle_measurements_;
mean_angle_shift_ = angle_cum / RealValue<sym>{static_cast<double>(n_voltage_angle_measurements_)};
}

static constexpr auto const is_measured = [](auto const& value) { return value >= 0; };
Expand Down Expand Up @@ -333,8 +333,8 @@ template <bool sym> class MeasuredValues {
// only direct injection
power_main_value_.push_back(direct_injection_measurement);
}
} else if (uncertain_direct_injection || any_zero(appliance_injection_measurement.p_variance) ||
any_zero(appliance_injection_measurement.q_variance)) {
} else if (uncertain_direct_injection || all_zero(appliance_injection_measurement.p_variance) ||
all_zero(appliance_injection_measurement.q_variance)) {
// only appliance injection if
// there is no direct injection measurement,
// or we have zero injection
Expand Down
597 changes: 379 additions & 218 deletions ...odel_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp
View file

Large diffs are not rendered by default.

6 changes: 5 additions & 1 deletion power_grid_model_c/power_grid_model/include/power_grid_model/three_phase_tensor.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,7 @@ template <scalar_value T> class Tensor : public Eigen3Tensor<T> {

template <scalar_value T> class DiagonalTensor : public Eigen3DiagonalTensor<T> {
public:
DiagonalTensor() { (*this) = Eigen3DiagonalTensor<T>::Zero(); }
DiagonalTensor() { (*this).setZero(); }
// additional constructors
explicit DiagonalTensor(T const& x) : Eigen3DiagonalTensor<T>{x, x, x} {}
explicit DiagonalTensor(Vector<T> const& v) : Eigen3DiagonalTensor<T>{v(0), v(1), v(2)} {}
Expand Down Expand Up @@ -304,6 +304,10 @@ inline auto is_inf(RealValue<false> const& value) { return is_inf(value(0)) || i
inline auto any_zero(std::floating_point auto value) { return value == 0.0; }
inline auto any_zero(RealValue<false> const& value) { return (value == RealValue<false>{0.0}).any(); }

// all_zero
inline auto all_zero(std::floating_point auto value) { return value == 0.0; }
inline auto all_zero(RealValue<false> const& value) { return (value == RealValue<false>{0.0}).all(); }

// update real values
//
// RealValue is only updated when the update value is not nan
Expand Down
45 changes: 45 additions & 0 deletions tests/cpp_unit_tests/test_main_model_se.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -211,6 +211,51 @@ TEST_CASE_TEMPLATE("Test main model - state estimation", CalculationMethod, Iter
}
}
}
SUBCASE("Line power sensor") {
main_model.add_component<Node>({{1, 10e3}, {2, 10e3}});
main_model.add_component<Line>({{3, 1, 2, 1, 1, 0.01, 0.01, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1e3}});
main_model.add_component<Source>({{4, 1, 1, 1.0, nan, nan, nan, nan}});
main_model.add_component<Shunt>({{6, 2, 1, 1800 / 10e3 / 10e3, -180 / 10e3 / 10e3, 0.0, 0.0}});
main_model.add_component<SymVoltageSensor>({{11, 1, 1e2, 10.0e3, 0.0}});
SUBCASE("Symmetric Power Sensor - Symmetric Calculation") {
main_model.add_component<SymPowerSensor>(
{{17, 3, MeasuredTerminalType::branch_from, 1e2, 1800.0, 180.0, nan, nan},
{18, 3, MeasuredTerminalType::branch_to, 1e2, -1800.0, -180.0, nan, nan},
{16, 6, MeasuredTerminalType::shunt, 1e2, 1800.0, 180.0, nan, nan}});
SUBCASE("Line flow") {
main_model.set_construction_complete();
std::vector<MathOutput<true>> const math_output =
main_model.calculate_state_estimation<true>(1e-8, 20, calculation_method);

std::vector<SymApplianceOutput> shunt_output(1);
std::vector<SymNodeOutput> node_output(2);
std::vector<SymPowerSensorOutput> power_sensor_output(3);
std::vector<BranchOutput<true>> line_output(1);
main_model.output_result<Shunt>(math_output, shunt_output.begin());
main_model.output_result<Node>(math_output, node_output.begin());
main_model.output_result<Line>(math_output, line_output.begin());
main_model.output_result<SymPowerSensor>(math_output, power_sensor_output.begin());

CHECK(shunt_output[0].p == doctest::Approx(1800.0).epsilon(0.01));
CHECK(shunt_output[0].q == doctest::Approx(180.0).epsilon(0.01));

CHECK(line_output[0].p_from == doctest::Approx(1800.0).epsilon(0.01));
CHECK(line_output[0].q_from == doctest::Approx(180.0).epsilon(0.01));
CHECK(line_output[0].p_to == doctest::Approx(-1800.0).epsilon(0.01));
CHECK(line_output[0].q_to == doctest::Approx(-180.0).epsilon(0.01));

// dealing with orders of magnitude kW / kVA and precision at W / VA level
auto const zero_at_order_of_magnitude = doctest::Approx(0.0).scale(1e3).epsilon(0.001);

CHECK(power_sensor_output[0].p_residual == zero_at_order_of_magnitude); // shunt
CHECK(power_sensor_output[0].q_residual == zero_at_order_of_magnitude); // shunt
CHECK(power_sensor_output[1].p_residual == zero_at_order_of_magnitude); // branch_from
CHECK(power_sensor_output[1].q_residual == zero_at_order_of_magnitude); // branch_from
CHECK(power_sensor_output[2].p_residual == zero_at_order_of_magnitude); // branch_to
CHECK(power_sensor_output[2].q_residual == zero_at_order_of_magnitude); // branch_to
}
}
}
SUBCASE("Forbid Link Power Measurements") {
main_model.add_component<Node>({{1, 10e3}, {2, 10e3}});
main_model.add_component<Link>({{3, 1, 2, 1, 1}});
Expand Down
Loading

0 comments on commit 5eddc40

Please sign in to comment.