diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp index 772b9edf93..904b64401a 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp @@ -182,6 +182,10 @@ struct MathModelTopology { Idx n_bus_power_sensor() const { return power_sensors_per_bus.element_size(); } + Idx n_branch_from_current_sensor() const { return current_sensors_per_branch_from.element_size(); } + + Idx n_branch_to_current_sensor() const { return current_sensors_per_branch_to.element_size(); } + Idx n_transformer_tap_regulator() const { return tap_regulators_per_branch.element_size(); } }; @@ -235,6 +239,8 @@ template struct StateEstimationInput { std::vector> measured_branch_from_power; std::vector> measured_branch_to_power; std::vector> measured_bus_injection; + std::vector> measured_branch_from_current; + std::vector> measured_branch_to_current; }; struct ShortCircuitInput { diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/common/exception.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/common/exception.hpp index 084108a4be..a5ed7eec80 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/common/exception.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/common/exception.hpp @@ -39,13 +39,13 @@ class InvalidArguments : public PowerGridError { }; template ... Options> - InvalidArguments(std::string const& method, std::string const& arguments) { + InvalidArguments(std::string_view method, std::string_view arguments) { append_msg(std::format("{} is not implemented for {}!\n", method, arguments)); } template requires(std::same_as, TypeValuePair> && ...) - InvalidArguments(std::string const& method, Options const&... options) + InvalidArguments(std::string_view method, Options const&... options) : InvalidArguments{method, "the following combination of options"} { (append_msg(std::format(" {}: {}\n", options.name, options.value)), ...); } @@ -54,7 +54,7 @@ class InvalidArguments : public PowerGridError { class MissingCaseForEnumError : public InvalidArguments { public: template - MissingCaseForEnumError(std::string const& method, T const& value) + MissingCaseForEnumError(std::string_view method, T const& value) : InvalidArguments{method, std::format("{} #{}", typeid(T).name(), detail::to_string(static_cast(value)))} {} }; @@ -97,7 +97,7 @@ class InvalidTransformerClock : public PowerGridError { class SparseMatrixError : public PowerGridError { public: - SparseMatrixError(Idx err, std::string const& msg = "") { + SparseMatrixError(Idx err, std::string_view msg = "") { append_msg( std::format("Sparse matrix error with error code #{} (possibly singular)\n", detail::to_string(err))); if (!msg.empty()) { @@ -117,7 +117,7 @@ class SparseMatrixError : public PowerGridError { class NotObservableError : public PowerGridError { public: - NotObservableError(std::string const& msg = "") { + NotObservableError(std::string_view msg = "") { append_msg("Not enough measurements available for state estimation.\n"); if (!msg.empty()) { append_msg(std::format("{}\n", msg)); @@ -137,7 +137,7 @@ class IterationDiverge : public PowerGridError { class MaxIterationReached : public IterationDiverge { public: - MaxIterationReached(std::string const& msg = "") { + MaxIterationReached(std::string_view msg = "") { append_msg(std::format("Maximum number of iterations reached{}\n", msg)); } }; @@ -161,14 +161,14 @@ class Idx2DNotFound : public PowerGridError { class InvalidMeasuredObject : public PowerGridError { public: - InvalidMeasuredObject(std::string const& object, std::string const& sensor) { + InvalidMeasuredObject(std::string_view object, std::string_view sensor) { append_msg(std::format("{} measurement is not supported for object of type {}", sensor, object)); } }; class InvalidMeasuredTerminalType : public PowerGridError { public: - InvalidMeasuredTerminalType(MeasuredTerminalType const terminal_type, std::string const& sensor) { + InvalidMeasuredTerminalType(MeasuredTerminalType const terminal_type, std::string_view sensor) { append_msg(std::format("{} measurement is not supported for object of type {}", sensor, detail::to_string(static_cast(terminal_type)))); } @@ -176,10 +176,10 @@ class InvalidMeasuredTerminalType : public PowerGridError { class InvalidRegulatedObject : public PowerGridError { public: - InvalidRegulatedObject(std::string const& object, std::string const& regulator) { + InvalidRegulatedObject(std::string_view object, std::string_view regulator) { append_msg(std::format("{} regulator is not supported for object of type {}", regulator, object)); } - InvalidRegulatedObject(ID id, std::string const& regulator) { + InvalidRegulatedObject(ID id, std::string_view regulator) { append_msg( std::format("{} regulator is not supported for object with ID {}", regulator, detail::to_string(id))); } @@ -203,7 +203,7 @@ class AutomaticTapCalculationError : public PowerGridError { class AutomaticTapInputError : public PowerGridError { public: - AutomaticTapInputError(std::string const& msg) { + AutomaticTapInputError(std::string_view msg) { append_msg(std::format("Automatic tap changer has invalid configuration. {}", msg)); } }; @@ -215,14 +215,21 @@ class IDWrongType : public PowerGridError { } }; +class ConflictingAngleMeasurementType : public PowerGridError { + public: + ConflictingAngleMeasurementType(std::string_view msg) { + append_msg(std::format("Conflicting angle measurement type. {}", msg)); + } +}; + class CalculationError : public PowerGridError { public: - explicit CalculationError(std::string const& msg) { append_msg(msg); } + explicit CalculationError(std::string_view msg) { append_msg(msg); } }; class BatchCalculationError : public CalculationError { public: - BatchCalculationError(std::string const& msg, IdxVector failed_scenarios, std::vector err_msgs) + BatchCalculationError(std::string_view msg, IdxVector failed_scenarios, std::vector err_msgs) : CalculationError(msg), failed_scenarios_{std::move(failed_scenarios)}, err_msgs_(std::move(err_msgs)) {} IdxVector const& failed_scenarios() const { return failed_scenarios_; } @@ -270,12 +277,12 @@ class InvalidShortCircuitPhaseOrType : public PowerGridError { class SerializationError : public PowerGridError { public: - explicit SerializationError(std::string const& msg) { append_msg(msg); } + explicit SerializationError(std::string_view msg) { append_msg(msg); } }; class DatasetError : public PowerGridError { public: - explicit DatasetError(std::string const& msg) { append_msg(std::format("Dataset error: {}", msg)); } + explicit DatasetError(std::string_view msg) { append_msg(std::format("Dataset error: {}", msg)); } }; class ExperimentalFeature : public InvalidArguments { @@ -289,7 +296,7 @@ class NotImplementedError : public PowerGridError { class UnreachableHit : public PowerGridError { public: - UnreachableHit(std::string const& method, std::string const& reason_for_assumption) { + UnreachableHit(std::string_view method, std::string_view reason_for_assumption) { append_msg(std::format("Unreachable code hit when executing {}.\n The following assumption for unreachability " "was not met: {}.\n This may be a bug in the library\n", method, reason_for_assumption)); @@ -299,7 +306,7 @@ class UnreachableHit : public PowerGridError { class TapSearchStrategyIncompatibleError : public InvalidArguments { public: template - TapSearchStrategyIncompatibleError(std::string const& method, T1 const& value1, T2 const& value2) + TapSearchStrategyIncompatibleError(std::string_view method, T1 const& value1, T2 const& value2) : InvalidArguments{method, std::format("{} #{} and {} #{}", typeid(T1).name(), detail::to_string(static_cast(value1)), typeid(T2).name(), detail::to_string(static_cast(value2)))} {} diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp index a4fdab9695..247a9b6764 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/common/statistics.hpp @@ -7,6 +7,8 @@ #include "common.hpp" #include "three_phase_tensor.hpp" +#include + /** * @file statistics.hpp * @brief This file contains various structures and functions for handling statistical representations of @@ -204,4 +206,87 @@ template struct PolarComplexRandVar { 9.0}}; } }; + +namespace statistics { +// combine multiple random variables of one quantity using Kalman filter +template + requires std::same_as, + UniformRealRandVar::sym>> || + std::same_as, + IndependentRealRandVar::sym>> || + std::same_as, + UniformComplexRandVar::sym>> || + std::same_as, + IndependentComplexRandVar::sym>> +constexpr auto combine(RandVarsView rand_vars) { + using RandVarType = std::ranges::range_value_t; + using ValueType = decltype(RandVarType::value); + using VarianceType = decltype(RandVarType::variance); + + VarianceType accumulated_inverse_variance{}; + ValueType weighted_accumulated_value{}; + + std::ranges::for_each(rand_vars, + [&accumulated_inverse_variance, &weighted_accumulated_value](auto const& measurement) { + accumulated_inverse_variance += VarianceType{1.0} / measurement.variance; + weighted_accumulated_value += measurement.value / measurement.variance; + }); + + if (!is_normal(accumulated_inverse_variance)) { + return RandVarType{.value = weighted_accumulated_value, + .variance = VarianceType{std::numeric_limits::infinity()}}; + } + return RandVarType{.value = weighted_accumulated_value / accumulated_inverse_variance, + .variance = VarianceType{1.0} / accumulated_inverse_variance}; +} + +template + requires std::same_as, + DecomposedComplexRandVar::sym>> +constexpr auto combine(RandVarsView rand_vars) { + using sym = std::ranges::range_value_t::sym; + + DecomposedComplexRandVar result{ + .real_component = + combine(rand_vars | std::views::transform([](auto const& x) -> auto& { return x.real_component; })), + .imag_component = + combine(rand_vars | std::views::transform([](auto const& x) -> auto& { return x.imag_component; }))}; + + if (!(is_normal(result.real_component.variance) && is_normal(result.imag_component.variance))) { + result.real_component.variance = RealValue{std::numeric_limits::infinity()}; + result.imag_component.variance = RealValue{std::numeric_limits::infinity()}; + } + return result; +} + +namespace detail { +template inline RealValue cabs_or_real(ComplexValue const& value) { + if (is_nan(imag(value))) { + return real(value); // only keep real part + } + return cabs(value); // get abs of the value +} +} // namespace detail + +template + requires std::same_as, + UniformComplexRandVar::sym>> +constexpr auto combine_magnitude(RandVarsView rand_vars) { + using sym = std::ranges::range_value_t::sym; + + auto const weighted_average_magnitude_measurement = + statistics::combine(rand_vars | std::views::transform([](auto const& measurement) { + return UniformRealRandVar{.value = detail::cabs_or_real(measurement.value), + .variance = measurement.variance}; + })); + return UniformComplexRandVar{.value = + [&weighted_average_magnitude_measurement]() { + ComplexValue abs_value = + piecewise_complex_value(DoubleComplex{0.0, nan}); + abs_value += weighted_average_magnitude_measurement.value; + return abs_value; + }(), + .variance = weighted_average_magnitude_measurement.variance}; +} +} // namespace statistics } // namespace power_grid_model diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/main_model_impl.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/main_model_impl.hpp index cd54e8e7b8..4607e7f6ff 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/main_model_impl.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/main_model_impl.hpp @@ -1171,6 +1171,8 @@ class MainModelImpl, ComponentLis se_input[i].measured_branch_from_power.resize(state.math_topology[i]->n_branch_from_power_sensor()); se_input[i].measured_branch_to_power.resize(state.math_topology[i]->n_branch_to_power_sensor()); se_input[i].measured_bus_injection.resize(state.math_topology[i]->n_bus_power_sensor()); + se_input[i].measured_branch_from_current.resize(state.math_topology[i]->n_branch_from_current_sensor()); + se_input[i].measured_branch_to_current.resize(state.math_topology[i]->n_branch_to_current_sensor()); } prepare_input_status::shunt_status, Shunt>(state, state.topo_comp_coup->shunt, @@ -1217,6 +1219,22 @@ class MainModelImpl, ComponentLis state, state.topo_comp_coup->power_sensor, se_input, [&state](Idx i) { return state.comp_topo->power_sensor_terminal_type[i] == MeasuredTerminalType::node; }); + prepare_input, CurrentSensorCalcParam, + &StateEstimationInput::measured_branch_from_current, GenericCurrentSensor>( + state, state.topo_comp_coup->current_sensor, se_input, [&state](Idx i) { + using enum MeasuredTerminalType; + return state.comp_topo->current_sensor_terminal_type[i] == branch_from || + // all branch3 sensors are at from side in the mathematical model + state.comp_topo->current_sensor_terminal_type[i] == branch3_1 || + state.comp_topo->current_sensor_terminal_type[i] == branch3_2 || + state.comp_topo->current_sensor_terminal_type[i] == branch3_3; + }); + prepare_input, CurrentSensorCalcParam, + &StateEstimationInput::measured_branch_to_current, GenericCurrentSensor>( + state, state.topo_comp_coup->current_sensor, se_input, [&state](Idx i) { + return state.comp_topo->current_sensor_terminal_type[i] == MeasuredTerminalType::branch_to; + }); + return se_input; } diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp index f99003a5d4..c10d31c1e5 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp @@ -142,7 +142,8 @@ template class IterativeLinearSESolver { private: // array selection function pointer - static constexpr std::array has_branch_{&MeasuredValues::has_branch_from, &MeasuredValues::has_branch_to}; + static constexpr std::array has_branch_power_{&MeasuredValues::has_branch_from_power, + &MeasuredValues::has_branch_to_power}; static constexpr std::array branch_power_{&MeasuredValues::branch_from_power, &MeasuredValues::branch_to_power}; @@ -210,7 +211,7 @@ template class IterativeLinearSESolver { // measured at from-side: 0, to-side: 1 for (IntS const measured_side : std::array{0, 1}) { // has measurement - if (std::invoke(has_branch_[measured_side], measured_value, obj)) { + 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); block.g() += @@ -301,7 +302,7 @@ template class IterativeLinearSESolver { // measured at from-side: 0, to-side: 1 for (IntS const measured_side : std::array{0, 1}) { // has measurement - if (std::invoke(has_branch_[measured_side], measured_value, obj)) { + if (std::invoke(has_branch_power_[measured_side], measured_value, obj)) { PowerSensorCalcParam const& m = std::invoke(branch_power_[measured_side], measured_value, obj); // the current needs to be calculated with the voltage of the measured bus side diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp index 2c69821f06..01704ec7ee 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp @@ -9,21 +9,12 @@ Collect all measured Values */ #include "../calculation_parameters.hpp" +#include "../common/exception.hpp" #include "../common/three_phase_tensor.hpp" #include namespace power_grid_model::math_solver { - -namespace detail { -template inline RealValue cabs_or_real(ComplexValue const& value) { - if (is_nan(imag(value))) { - return real(value); // only keep real part - } - return cabs(value); // get abs of the value -} -} // namespace detail - // processed measurement struct // combined all measurement of the same quantity // accumulate for bus injection measurement @@ -54,6 +45,8 @@ template class MeasuredValues { idx_shunt_power_(math_topology().n_shunt()), idx_load_gen_power_(math_topology().n_load_gen()), idx_source_power_(math_topology().n_source()), + idx_branch_from_current_(math_topology().n_branch()), + idx_branch_to_current_(math_topology().n_branch()), // default angle shift // sym: 0 // asym: 0, -120deg, -240deg @@ -68,12 +61,15 @@ template class MeasuredValues { constexpr bool has_angle() const { return n_voltage_angle_measurements_ > 0; } constexpr bool has_voltage_measurements() const { return n_voltage_measurements_ > 0; } + constexpr bool has_global_angle_current() const { return n_global_angle_current_measurements_ > 0; } constexpr bool has_voltage(Idx bus) const { return idx_voltage_[bus] >= 0; } constexpr bool has_angle_measurement(Idx bus) const { return !is_nan(imag(voltage(bus))); } constexpr bool has_bus_injection(Idx bus) const { return bus_injection_[bus].idx_bus_injection >= 0; } - constexpr bool has_branch_from(Idx branch) const { return idx_branch_from_power_[branch] >= 0; } - constexpr bool has_branch_to(Idx branch) const { return idx_branch_to_power_[branch] >= 0; } + constexpr bool has_branch_from_power(Idx branch) const { return idx_branch_from_power_[branch] >= 0; } + constexpr bool has_branch_to_power(Idx branch) const { return idx_branch_to_power_[branch] >= 0; } + constexpr bool has_branch_from_current(Idx branch) const { return idx_branch_from_current_[branch] >= 0; } + constexpr bool has_branch_to_current(Idx branch) const { return idx_branch_to_current_[branch] >= 0; } constexpr bool has_shunt(Idx shunt) const { return idx_shunt_power_[shunt] >= 0; } constexpr bool has_load_gen(Idx load_gen) const { return idx_load_gen_power_[load_gen] >= 0; } constexpr bool has_source(Idx source) const { return idx_source_power_[source] >= 0; } @@ -91,6 +87,12 @@ template class MeasuredValues { return power_main_value_[idx_branch_from_power_[branch]]; } constexpr auto const& branch_to_power(Idx branch) const { return power_main_value_[idx_branch_to_power_[branch]]; } + constexpr auto const& branch_from_current(Idx branch) const { + return current_main_value_[idx_branch_from_current_[branch]]; + } + constexpr auto const& branch_to_current(Idx branch) const { + return current_main_value_[idx_branch_to_current_[branch]]; + } constexpr auto const& shunt_power(Idx shunt) const { return power_main_value_[idx_shunt_power_[shunt]]; } constexpr auto const& load_gen_power(Idx load_gen) const { return extra_value_[idx_load_gen_power_[load_gen]]; } constexpr auto const& source_power(Idx source) const { return extra_value_[idx_source_power_[source]]; } @@ -177,6 +179,7 @@ template class MeasuredValues { // branch/shunt flow, bus voltage, injection flow std::vector> voltage_main_value_; std::vector> power_main_value_; + std::vector> current_main_value_; // flat array of all the load_gen/source measurement // not relevant for the main calculation, as extra data for load_gen/source calculation @@ -197,9 +200,13 @@ template class MeasuredValues { // relevant for extra value IdxVector idx_load_gen_power_; IdxVector idx_source_power_; + // current measurement + IdxVector idx_branch_from_current_; + IdxVector idx_branch_to_current_; Idx n_voltage_measurements_{}; Idx n_voltage_angle_measurements_{}; + Idx n_global_angle_current_measurements_{}; // average angle shift of voltages with angle measurement // default is zero is no voltage has angle measurement @@ -211,7 +218,7 @@ template class MeasuredValues { void process_bus_related_measurements(StateEstimationInput const& input) { /* - The main purpose of this function is to aggregate all voltage and power sensor values to + The main purpose of this function is to aggregate all voltage and power/current sensor values to one voltage sensor value per bus. one injection power sensor value per bus. one power sensor value per shunt (in injection reference direction, note shunt itself is not considered as @@ -383,36 +390,49 @@ template class MeasuredValues { void process_branch_measurements(StateEstimationInput const& input) { /* - The main purpose of this function is to aggregate all power sensor values to one power sensor value per branch - side. + The main purpose of this function is to aggregate all power/current sensor values to one power/current sensor + value per branch side. This function loops through all branches. The branch_bus_idx contains the from and to bus indexes of the branch, or disconnected if the branch is not connected at that side. For each branch the checker checks if the from and to side are connected by checking if branch_bus_idx = disconnected. - If the branch_bus_idx = disconnected, idx_branch_to_power_/idx_branch_from_power_ is set to disconnected. + If the branch_bus_idx = disconnected, idx_branch_(to/from)_(power/current)_ is set to disconnected. If the side is connected, but there are no measurements in this branch side - idx_branch_to_power_/idx_branch_from_power_ is set to disconnected. - Else, idx_branch_to_power_/idx_branch_from_power_ is set to the index of the aggregated data in - power_main_value_. + idx_branch_(to/from)_(power/current)_ is set to disconnected. + Else, idx_branch_(to/from)_(power/current)_ is set to the index of the aggregated data in + power/current_main_value_. All measurement values for a single side of a branch are combined in a weighted average, which is appended to - power_main_value_. The values in power_main_value_ can be found using - idx_branch_to_power_/idx_branch_from_power_. + power/current_main_value_. The values in power/current_main_value_ can be found using + idx_branch_(to/from)_(power/current)_. */ MathModelTopology const& topo = math_topology(); static constexpr auto branch_from_checker = [](BranchIdx x) { return x[0] != -1; }; static constexpr auto branch_to_checker = [](BranchIdx x) { return x[1] != -1; }; for (Idx const branch : IdxRange{topo.n_branch()}) { - // from side + // from side power idx_branch_from_power_[branch] = process_one_object(branch, topo.power_sensors_per_branch_from, topo.branch_bus_idx, input.measured_branch_from_power, power_main_value_, branch_from_checker); - // to side + // to side power idx_branch_to_power_[branch] = process_one_object(branch, topo.power_sensors_per_branch_to, topo.branch_bus_idx, input.measured_branch_to_power, power_main_value_, branch_to_checker); + // from side current + idx_branch_from_current_[branch] = + process_one_object(branch, topo.current_sensors_per_branch_from, topo.branch_bus_idx, + input.measured_branch_from_current, current_main_value_, branch_from_checker); + // to side current + idx_branch_to_current_[branch] = + process_one_object(branch, topo.current_sensors_per_branch_to, topo.branch_bus_idx, + input.measured_branch_to_current, current_main_value_, branch_to_checker); + + n_global_angle_current_measurements_ = + std::ranges::count_if(current_main_value_, [](auto const& measurement) { + return measurement.angle_measurement_type == AngleMeasurementType::global_angle; + }); } } @@ -423,70 +443,36 @@ template class MeasuredValues { template static VoltageSensorCalcParam combine_measurements(std::vector> const& data, IdxRange const& sensors) { - double accumulated_inverse_variance{}; - ComplexValue accumulated_value{}; - for (auto pos : sensors) { - auto const& measurement = data[pos]; - auto const inv_variance = 1.0 / measurement.variance; - - accumulated_inverse_variance += inv_variance; - if constexpr (only_magnitude) { - ComplexValue abs_value = piecewise_complex_value(DoubleComplex{0.0, nan}); - abs_value += detail::cabs_or_real(measurement.value); - accumulated_value += abs_value * inv_variance; - } else { - // accumulate value - accumulated_value += measurement.value * inv_variance; - } - } - - if (!is_normal(accumulated_inverse_variance)) { - return VoltageSensorCalcParam{accumulated_value, std::numeric_limits::infinity()}; + auto complex_measurements = sensors | std::views::transform([&data](Idx pos) -> auto& { return data[pos]; }); + if constexpr (only_magnitude) { + return statistics::combine_magnitude(complex_measurements); + } else { + return statistics::combine(complex_measurements); } - - return VoltageSensorCalcParam{accumulated_value / accumulated_inverse_variance, - 1.0 / accumulated_inverse_variance}; } - - // combine multiple measurements of one quantity - // using Kalman filter - // if only_magnitude = true, combine the abs value of the individual data - // set imag part to nan, to signal this is a magnitude only measurement - // - // Since p and q are entirely decoupled, the real and imaginary components accumulate separately template requires(!only_magnitude) static PowerSensorCalcParam combine_measurements(std::vector> const& data, IdxRange const& sensors) { - struct AccumulatedValues { - RealValue inverse_p_variance{}; - RealValue inverse_q_variance{}; - RealValue p_value{}; - RealValue q_value{}; - }; - - AccumulatedValues accumulated; - std::ranges::for_each(sensors, [&data, &accumulated](auto pos) { - DecomposedComplexRandVar const& measurement = data[pos]; - accumulated.inverse_p_variance += RealValue{1.0} / measurement.real_component.variance; - accumulated.inverse_q_variance += RealValue{1.0} / measurement.imag_component.variance; - - accumulated.p_value += measurement.real_component.value / measurement.real_component.variance; - accumulated.q_value += measurement.imag_component.value / measurement.imag_component.variance; - }); - - if (is_normal(accumulated.inverse_p_variance) && is_normal(accumulated.inverse_q_variance)) { - return PowerSensorCalcParam{ - .real_component = {.value = accumulated.p_value / accumulated.inverse_p_variance, - .variance = RealValue{1.0} / accumulated.inverse_p_variance}, - .imag_component = {.value = accumulated.q_value / accumulated.inverse_q_variance, - .variance = RealValue{1.0} / accumulated.inverse_q_variance}}; + return statistics::combine(sensors | std::views::transform([&data](Idx pos) -> auto& { return data[pos]; })); + } + template + requires(!only_magnitude) + static CurrentSensorCalcParam combine_measurements(std::vector> const& data, + IdxRange const& sensors) { + auto const params = sensors | std::views::transform([&data](Idx pos) -> auto& { return data[pos]; }); + auto const angle_measurement_type = sensors.empty() ? AngleMeasurementType::local_angle // fallback + : params.front().angle_measurement_type; + if (std::ranges::any_of(params, [angle_measurement_type](auto const& param) { + return param.angle_measurement_type != angle_measurement_type; + })) { + throw ConflictingAngleMeasurementType{ + "Cannot mix local and global angle current measurements on the same terminal."}; } - return PowerSensorCalcParam{ - .real_component = {.value = accumulated.p_value, - .variance = RealValue{std::numeric_limits::infinity()}}, - .imag_component = {.value = accumulated.q_value, - .variance = RealValue{std::numeric_limits::infinity()}}}; + + return {.angle_measurement_type = angle_measurement_type, + .measurement = statistics::combine( + params | std::views::transform([](auto const& param) -> auto& { return param.measurement; }))}; } template @@ -512,11 +498,10 @@ template class MeasuredValues { static constexpr DefaultStatusChecker default_status_checker{}; - template + template static Idx process_one_object(Idx const object, grouped_idx_vector_type auto const& sensors_per_object, - std::vector const& object_status, - std::vector> const& input_data, - std::vector>& result_data, + std::vector const& object_status, std::vector const& input_data, + std::vector& result_data, StatusChecker status_checker = default_status_checker) { if (!status_checker(object_status[object])) { return disconnected; @@ -555,6 +540,16 @@ template class MeasuredValues { } } } + for (auto const& x : current_main_value_) { + auto const variance = x.measurement.real_component.variance + x.measurement.imag_component.variance; + if constexpr (is_symmetric_v) { + unconstrained_min(variance); + } else { + for (Idx const phase : {0, 1, 2}) { + unconstrained_min(variance[phase]); + } + } + } // scale auto const inv_norm_var = 1.0 / min_var; @@ -563,6 +558,10 @@ template class MeasuredValues { x.real_component.variance *= inv_norm_var; x.imag_component.variance *= inv_norm_var; }); + std::ranges::for_each(current_main_value_, [inv_norm_var](auto& x) { + x.measurement.real_component.variance *= inv_norm_var; + x.measurement.imag_component.variance *= inv_norm_var; + }); } void calculate_non_over_determined_injection(Idx n_unmeasured, IdxRange const& load_gens, IdxRange const& sources, diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp index 809734124a..5878bb3b17 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp @@ -212,6 +212,8 @@ template class NewtonRaphsonSESolver { typename SparseLUSolver, NRSERhs, NRSEUnknown>::BlockPermArray perm_; void initialize_unknown(ComplexValueVector& initial_u, MeasuredValues const& measured_values) { + using statistics::detail::cabs_or_real; + reset_unknown(); RealValue const mean_angle_shift = measured_values.mean_angle_shift(); for (Idx bus = 0; bus != n_bus_; ++bus) { @@ -222,7 +224,7 @@ template class NewtonRaphsonSESolver { if (measured_values.has_angle_measurement(bus)) { estimated_result.theta() = arg(measured_values.voltage(bus)); } - estimated_result.v() = detail::cabs_or_real(measured_values.voltage(bus)); + estimated_result.v() = cabs_or_real(measured_values.voltage(bus)); } initial_u[bus] = estimated_result.v() * exp(1.0i * estimated_result.theta()); } @@ -304,14 +306,14 @@ template class NewtonRaphsonSESolver { [[fallthrough]]; case YBusElementType::btf: { auto const& y_branch = param.branch_param[obj]; - if (measured_values.has_branch_from(obj)) { + if (measured_values.has_branch_from_power(obj)) { auto const ij_voltage_order = (type == YBusElementType::bft) ? Order::row_major : Order::column_major; process_branch_measurement(block, diag_block, rhs_block, y_branch.yff(), y_branch.yft(), u_state, ij_voltage_order, measured_values.branch_from_power(obj)); } - if (measured_values.has_branch_to(obj)) { + if (measured_values.has_branch_to_power(obj)) { auto const ij_voltage_order = (type == YBusElementType::btf) ? Order::row_major : Order::column_major; process_branch_measurement(block, diag_block, rhs_block, y_branch.ytt(), y_branch.ytf(), @@ -530,6 +532,8 @@ template class NewtonRaphsonSESolver { /// @param bus bus with voltage measurement void process_voltage_measurements(NRSEGainBlock& block, NRSERhs& rhs_block, MeasuredValues const& measured_values, Idx const& bus) { + using statistics::detail::cabs_or_real; + if (!measured_values.has_voltage(bus)) { return; } @@ -537,7 +541,7 @@ template class NewtonRaphsonSESolver { // G += 1.0 / variance // for 3x3 tensor, fill diagonal auto const w_v = RealTensor{1.0 / measured_values.voltage_var(bus)}; - auto const abs_measured_v = detail::cabs_or_real(measured_values.voltage(bus)); + auto const abs_measured_v = cabs_or_real(measured_values.voltage(bus)); auto const delta_v = abs_measured_v - x_[bus].v(); auto const virtual_angle_measurement_bus = measured_values.has_voltage(math_topo_->slack_bus) diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp index 7e855f5d95..1c9e1493ef 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp @@ -63,7 +63,9 @@ std::pair, bool> count_flow_sensors(MeasuredValues cons // we only need one flow sensor, so the loop will break if (element.element_type != YBusElementType::shunt) { Idx const branch = element.idx; - if ((measured_values.has_branch_from(branch) || measured_values.has_branch_to(branch)) && + if ((measured_values.has_branch_from_power(branch) || measured_values.has_branch_to_power(branch) || + measured_values.has_branch_from_current(branch) || + measured_values.has_branch_to_current(branch)) && topo.branch_bus_idx[branch][0] != -1 && topo.branch_bus_idx[branch][1] != -1) { flow_sensors[ybus_index] = 1; has_at_least_one_sensor = true; @@ -153,6 +155,11 @@ inline ObservabilityResult necessary_observability_check(MeasuredValues con throw NotObservableError{}; } + if (measured_values.has_global_angle_current() && n_voltage_phasor_sensor == 0) { + throw NotObservableError{ + "Global angle current sensors require at least one voltage angle measurement as a reference point.\n"}; + } + // for radial grid without phasor measurement, try to assign injection sensor to branch sensor // we can then check sufficient condition for observability if (topo.is_radial && n_voltage_phasor_sensor == 0) { @@ -161,8 +168,9 @@ inline ObservabilityResult necessary_observability_check(MeasuredValues con if (Idx const n_flow_sensor_new = std::reduce(flow_sensors.cbegin(), flow_sensors.cend(), Idx{}, std::plus{}); n_flow_sensor_new < n_bus - 1) { - throw NotObservableError{"The number of power sensors appears sufficient, but they are not independent " - "enough. The system is still not observable.\n"}; + throw NotObservableError{ + "The number of power/current sensors appears sufficient, but they are not independent " + "enough. The system is still not observable.\n"}; } result.is_sufficiently_observable = true; } diff --git a/src/power_grid_model/_core/error_handling.py b/src/power_grid_model/_core/error_handling.py index 17d544fd01..17248926b3 100644 --- a/src/power_grid_model/_core/error_handling.py +++ b/src/power_grid_model/_core/error_handling.py @@ -16,6 +16,7 @@ AutomaticTapCalculationError, AutomaticTapInputError, ConflictID, + ConflictingAngleMeasurementType, ConflictVoltage, IDNotFound, IDWrongType, @@ -75,6 +76,7 @@ _AUTOMATIC_TAP_INPUT_ERROR_RE = re.compile(r"Automatic tap changer has invalid configuration") _ID_WRONG_TYPE_RE = re.compile(r"Wrong type for object with id (-?\d+)\n") +_CONFLICTING_ANGLE_MEASUREMENT_TYPE_RE = re.compile(r"Conflicting angle measurement type") _INVALID_CALCULATION_METHOD_RE = re.compile(r"The calculation method is invalid for this calculation!") _INVALID_SHORT_CIRCUIT_PHASE_OR_TYPE_RE = re.compile(r"short circuit type") # multiple different flavors _POWER_GRID_DATASET_ERROR_RE = re.compile(r"Dataset error: ") # multiple different flavors @@ -100,6 +102,7 @@ _AUTOMATIC_TAP_CALCULATION_ERROR_RE: AutomaticTapCalculationError, _AUTOMATIC_TAP_INPUT_ERROR_RE: AutomaticTapInputError, _ID_WRONG_TYPE_RE: IDWrongType, + _CONFLICTING_ANGLE_MEASUREMENT_TYPE_RE: ConflictingAngleMeasurementType, _INVALID_CALCULATION_METHOD_RE: InvalidCalculationMethod, _INVALID_SHORT_CIRCUIT_PHASE_OR_TYPE_RE: InvalidShortCircuitPhaseOrType, _POWER_GRID_DATASET_ERROR_RE: PowerGridDatasetError, diff --git a/src/power_grid_model/enum.py b/src/power_grid_model/enum.py index d6ae356007..7eefd773bc 100644 --- a/src/power_grid_model/enum.py +++ b/src/power_grid_model/enum.py @@ -137,6 +137,19 @@ class MeasuredTerminalType(IntEnum): """ +class AngleMeasurementType(IntEnum): + """The type of the angle measurement for current sensors.""" + + local_angle = 0 + """ + The angle is relative to the local voltage angle + """ + global_angle = 1 + """ + The angle is relative to the global voltage angle. + """ + + class FaultType(IntEnum): """The type of fault represented by a fault component""" diff --git a/src/power_grid_model/errors.py b/src/power_grid_model/errors.py index 4576253906..cbca6bf97c 100644 --- a/src/power_grid_model/errors.py +++ b/src/power_grid_model/errors.py @@ -101,6 +101,10 @@ class AutomaticTapInputError(PowerGridError): """Automatic tap changer has invalid configuration.""" +class ConflictingAngleMeasurementType(PowerGridError): + """Conflicting angle measurement types found.""" + + class InvalidShortCircuitPhaseOrType(PowerGridError): """Invalid (combination of) short circuit types and phase(s) provided.""" diff --git a/tests/cpp_unit_tests/test_measured_values.cpp b/tests/cpp_unit_tests/test_measured_values.cpp index d59e7e87e2..f63592acb6 100644 --- a/tests/cpp_unit_tests/test_measured_values.cpp +++ b/tests/cpp_unit_tests/test_measured_values.cpp @@ -39,7 +39,7 @@ TEST_CASE("Measured Values") { MeasuredValues const values{std::make_shared(std::move(topo)), input}; - CHECK(values.has_bus_injection(0)); + REQUIRE(values.has_bus_injection(0)); auto const& injection = values.bus_injection(0); check_close(injection.value(), 1.0 + 0.1i); check_close(injection.real_component.variance, 0.75); @@ -61,7 +61,7 @@ TEST_CASE("Measured Values") { MeasuredValues const values{std::make_shared(std::move(topo)), input}; - CHECK(values.has_bus_injection(0)); + REQUIRE(values.has_bus_injection(0)); auto const& injection = values.bus_injection(0); check_close(injection.value(), ComplexValue{1.0 + 0.0i, 1.1 + 0.1i, 1.2 - 0.2i}); check_close(injection.real_component.variance, RealValue{0.75, 1.5, 1.625}); @@ -84,7 +84,7 @@ TEST_CASE("Measured Values") { MeasuredValues const values{std::make_shared(std::move(topo)), input}; - CHECK(values.has_bus_injection(0)); + REQUIRE(values.has_bus_injection(0)); auto const& injection = values.bus_injection(0); check_close(injection.value(), 2.0 + 1.0i); check_close(injection.real_component.variance, 16.0 / 61.0); @@ -107,7 +107,7 @@ TEST_CASE("Measured Values") { MeasuredValues const values{std::make_shared(std::move(topo)), input}; - CHECK(values.has_bus_injection(0)); + REQUIRE(values.has_bus_injection(0)); auto const& injection = values.bus_injection(0); check_close(injection.value(), 5.0 + 2.2i); check_close(injection.real_component.variance, 0.3); @@ -131,7 +131,7 @@ TEST_CASE("Measured Values") { MeasuredValues const values{std::make_shared(std::move(topo)), input}; - CHECK(values.has_bus_injection(0)); + REQUIRE(values.has_bus_injection(0)); auto const& injection = values.bus_injection(0); check_close(injection.value(), ComplexValue{5.0 + 2.2i, 0.1 + 0.2i, 0.8 + 3.0i}); check_close(injection.real_component.variance, RealValue{0.3, 0.2, 0.85}); @@ -139,4 +139,159 @@ TEST_CASE("Measured Values") { } } -} // namespace power_grid_model::math_solver \ No newline at end of file +TEST_CASE_TEMPLATE("Measured Values - Accumulate branch flow sensors", sym, symmetric_t, asymmetric_t) { + using enum AngleMeasurementType; + + auto const get_single_branch_topology = [] { + auto topo = MathModelTopology{}; + topo.phase_shift = {0.0, 0.0}; + topo.branch_bus_idx = {{0, 1}}; + topo.shunts_per_bus = {from_dense, {}, 1}; + topo.load_gens_per_bus = {from_dense, {}, 1}; + topo.sources_per_bus = {from_dense, {}, 1}; + return topo; + }; + auto const measurement_a = [] { + if constexpr (is_symmetric_v) { + return DecomposedComplexRandVar{.real_component = {.value = 1.0, .variance = 1.0}, + .imag_component = {.value = 1.5, .variance = 4.0}}; + } else { + return DecomposedComplexRandVar{ + .real_component = {.value = RealValue{1.0, 1.0, 1.0}, + .variance = RealValue{1.0, 1.0, 1.0}}, + .imag_component = {.value = RealValue{1.5, 1.5, 1.5}, + .variance = RealValue{4.0, 4.0, 4.0}}}; + } + }(); + auto const measurement_b = [] { + if constexpr (is_symmetric_v) { + return DecomposedComplexRandVar{.real_component = {.value = 4.0, .variance = 2.0}, + .imag_component = {.value = 0.7, .variance = 3.0}}; + } else { + return DecomposedComplexRandVar{ + .real_component = {.value = RealValue{4.0, 4.0, 4.0}, + .variance = RealValue{2.0, 2.0, 2.0}}, + .imag_component = {.value = RealValue{0.7, 0.7, 0.7}, + .variance = RealValue{3.0, 3.0, 3.0}}}; + } + }(); + + auto const check_accumulated = [](DecomposedComplexRandVar const& value) { + auto const expected_value = 2.0 + (73.0i / 70.0); + auto const expected_real_variance = 7.0 / 25.0; + auto const expected_imag_variance = 18.0 / 25.0; + if constexpr (is_symmetric_v) { + check_close(value.value(), expected_value); + check_close(value.real_component.variance, expected_real_variance); + check_close(value.imag_component.variance, expected_imag_variance); + } else { + for (Idx phase : IdxRange(3)) { + // eigen index-based element access + check_close(value.value()(phase), expected_value); + check_close(value.real_component.variance(phase), expected_real_variance); + check_close(value.imag_component.variance(phase), expected_imag_variance); + } + } + }; + + SUBCASE("Accumulate power sensors on branch") { + auto topo = get_single_branch_topology(); + topo.power_sensors_per_branch_from = {from_dense, {0, 0}, 1}; + + StateEstimationInput input{}; + input.measured_branch_from_power = {measurement_a, measurement_b}; + + MeasuredValues const values{std::make_shared(std::move(topo)), input}; + + REQUIRE(values.has_branch_from_power(0)); + CHECK_FALSE(values.has_branch_from_current(0)); + CHECK_FALSE(values.has_branch_to_power(0)); + CHECK_FALSE(values.has_branch_to_current(0)); + + check_accumulated(values.branch_from_power(0)); + } + + SUBCASE("Accumulate local angle current sensors on branch") { + auto topo = get_single_branch_topology(); + topo.current_sensors_per_branch_from = {from_dense, {0, 0}, 1}; + + StateEstimationInput input{}; + input.measured_branch_from_current = {{.angle_measurement_type = local_angle, .measurement = measurement_a}, + {.angle_measurement_type = local_angle, .measurement = measurement_b}}; + + MeasuredValues const values{std::make_shared(std::move(topo)), input}; + + REQUIRE(values.has_branch_from_current(0)); + CHECK_FALSE(values.has_branch_from_power(0)); + CHECK_FALSE(values.has_branch_to_power(0)); + CHECK_FALSE(values.has_branch_to_current(0)); + + CHECK(values.branch_from_current(0).angle_measurement_type == local_angle); + check_accumulated(values.branch_from_current(0).measurement); + } + + SUBCASE("Accumulate global angle current sensors on branch") { + auto topo = get_single_branch_topology(); + topo.current_sensors_per_branch_from = {from_dense, {0, 0}, 1}; + topo.voltage_sensors_per_bus = {from_dense, {0}, 2}; + + StateEstimationInput input{}; + ComplexValue const measured_voltage = [] { + if constexpr (is_symmetric_v) { + return ComplexValue{1.0, 1.0}; + } else { + return ComplexValue{RealValue{1.0, 1.0, 1.0}, + RealValue{1.0, 1.0, 1.0}}; + } + }(); + input.measured_voltage = {{.value = measured_voltage, .variance = 10.0}}; + input.measured_branch_from_current = {{.angle_measurement_type = global_angle, .measurement = measurement_a}, + {.angle_measurement_type = global_angle, .measurement = measurement_b}}; + + MeasuredValues const values{std::make_shared(std::move(topo)), input}; + + REQUIRE(values.has_branch_from_current(0)); + CHECK_FALSE(values.has_branch_from_power(0)); + CHECK_FALSE(values.has_branch_to_power(0)); + CHECK_FALSE(values.has_branch_to_current(0)); + + CHECK(values.branch_from_current(0).angle_measurement_type == global_angle); + check_accumulated(values.branch_from_current(0).measurement); + } + + SUBCASE("Accumulate local angle current sensors on branch") { + auto topo = get_single_branch_topology(); + topo.current_sensors_per_branch_from = {from_dense, {0, 0}, 1}; + + StateEstimationInput input{}; + input.measured_branch_from_current = {{.angle_measurement_type = local_angle, .measurement = measurement_a}, + {.angle_measurement_type = local_angle, .measurement = measurement_b}}; + + MeasuredValues const values{std::make_shared(std::move(topo)), input}; + + REQUIRE(values.has_branch_from_current(0)); + CHECK_FALSE(values.has_branch_from_power(0)); + CHECK_FALSE(values.has_branch_to_power(0)); + CHECK_FALSE(values.has_branch_to_current(0)); + + CHECK(values.branch_from_current(0).angle_measurement_type == local_angle); + check_accumulated(values.branch_from_current(0).measurement); + } + + SUBCASE("Cannot accumulate different angle measurement types on same terminal") { + auto topo = get_single_branch_topology(); + topo.current_sensors_per_branch_from = {from_dense, {0, 0}, 1}; + + StateEstimationInput input{}; + input.measured_branch_from_current = { + {.angle_measurement_type = AngleMeasurementType::local_angle, .measurement = measurement_a}, + {.angle_measurement_type = AngleMeasurementType::global_angle, .measurement = measurement_b}}; + + auto const create_measured_values = [&topo, &input] { + return MeasuredValues{std::make_shared(topo), input}; + }; + CHECK_THROWS_AS(create_measured_values(), ConflictingAngleMeasurementType); + } +} + +} // namespace power_grid_model::math_solver diff --git a/tests/cpp_unit_tests/test_observability.cpp b/tests/cpp_unit_tests/test_observability.cpp index fec8450788..9bc377a64e 100644 --- a/tests/cpp_unit_tests/test_observability.cpp +++ b/tests/cpp_unit_tests/test_observability.cpp @@ -11,16 +11,32 @@ namespace power_grid_model { namespace { -void check_not_observable(MathModelTopology const& topo, MathModelParam const& param, - StateEstimationInput const& se_input) { +void check_whether_observable(bool is_observable, MathModelTopology const& topo, + MathModelParam const& param, + StateEstimationInput const& se_input) { auto topo_ptr = std::make_shared(topo); auto param_ptr = std::make_shared const>(param); YBus const y_bus{topo_ptr, param_ptr}; math_solver::MeasuredValues const measured_values{y_bus.shared_topology(), se_input}; - CHECK_THROWS_AS( - math_solver::necessary_observability_check(measured_values, y_bus.math_topology(), y_bus.y_bus_structure()), - NotObservableError); + if (is_observable) { + CHECK_NOTHROW(math_solver::necessary_observability_check(measured_values, y_bus.math_topology(), + y_bus.y_bus_structure())); + } else { + CHECK_THROWS_AS( + math_solver::necessary_observability_check(measured_values, y_bus.math_topology(), y_bus.y_bus_structure()), + NotObservableError); + } +} + +void check_observable(MathModelTopology const& topo, MathModelParam const& param, + StateEstimationInput const& se_input) { + check_whether_observable(true, topo, param, se_input); +} + +void check_not_observable(MathModelTopology const& topo, MathModelParam const& param, + StateEstimationInput const& se_input) { + check_whether_observable(false, topo, param, se_input); } } // namespace @@ -45,12 +61,13 @@ TEST_CASE("Necessary observability check") { topo.power_sensors_per_shunt = {from_sparse, {0}}; topo.power_sensors_per_branch_from = {from_sparse, {0, 1, 1, 1}}; topo.power_sensors_per_branch_to = {from_sparse, {0, 0, 0, 0}}; + topo.current_sensors_per_branch_from = {from_sparse, {0, 0, 0, 0}}; + topo.current_sensors_per_branch_to = {from_sparse, {0, 0, 0, 0}}; topo.voltage_sensors_per_bus = {from_sparse, {0, 1, 1, 1}}; MathModelParam param; param.source_param = {SourceCalcParam{.y1 = 10.0 - 50.0i, .y0 = 10.0 - 50.0i}}; param.branch_param = {{1.0, -1.0, -1.0, 1.0}, {1.0, -1.0, -1.0, 1.0}, {1.0, -1.0, -1.0, 1.0}}; - auto param_ptr = std::make_shared const>(param); StateEstimationInput se_input; se_input.source_status = {1}; @@ -61,13 +78,7 @@ TEST_CASE("Necessary observability check") { se_input.measured_branch_from_power = { {.real_component = {.value = 1.0, .variance = 1.0}, .imag_component = {.value = 0.0, .variance = 1.0}}}; - SUBCASE("Observable grid") { - auto topo_ptr = std::make_shared(topo); - YBus const y_bus{topo_ptr, param_ptr}; - math_solver::MeasuredValues const measured_values{y_bus.shared_topology(), se_input}; - CHECK_NOTHROW(math_solver::necessary_observability_check(measured_values, y_bus.math_topology(), - y_bus.y_bus_structure())); - } + SUBCASE("Observable grid") { check_observable(topo, param, se_input); } SUBCASE("No voltage sensor") { topo.voltage_sensors_per_bus = {from_sparse, {0, 0, 0, 0}}; @@ -107,6 +118,59 @@ TEST_CASE("Necessary observability check") { // this will throw NotObservableError check_not_observable(topo, param, se_input); } -} + SUBCASE("Current sensors also measure branch flow") { + using enum AngleMeasurementType; + + topo.power_sensors_per_branch_from = {from_dense, {}, 3}; + se_input.measured_branch_from_power = {}; + topo.current_sensors_per_branch_from = {from_dense, {2}, 3}; + + DecomposedComplexRandVar const current_measurement{ + .real_component = {.value = 1.0, .variance = 1.0}, .imag_component = {.value = 0.0, .variance = 1.0}}; + + SUBCASE("With voltage phasor measurement") { + SUBCASE("Local current sensor") { + se_input.measured_branch_from_current = { + {.angle_measurement_type = local_angle, .measurement = current_measurement}}; + check_observable(topo, param, se_input); + } + SUBCASE("Global angle current sensor") { + se_input.measured_branch_from_current = { + {.angle_measurement_type = global_angle, .measurement = current_measurement}}; + check_observable(topo, param, se_input); + } + } + SUBCASE("No voltage phasor measurement and single current sensor") { + se_input.measured_voltage = {{.value = {1.0, nan}, .variance = 1.0}}; + SUBCASE("Local current sensor") { + se_input.measured_branch_from_current = { + {.angle_measurement_type = local_angle, .measurement = current_measurement}}; + check_not_observable(topo, param, se_input); + } + SUBCASE("Global angle current sensor") { + se_input.measured_branch_from_current = { + {.angle_measurement_type = global_angle, .measurement = current_measurement}}; + check_not_observable(topo, param, se_input); + } + } + SUBCASE("No voltage phasor measurement and two current sensors") { + se_input.measured_voltage = {{.value = {1.0, nan}, .variance = 1.0}}; + topo.current_sensors_per_branch_from = {from_dense, {0, 2}, 3}; + + SUBCASE("Local current sensor") { + se_input.measured_branch_from_current = { + {.angle_measurement_type = local_angle, .measurement = current_measurement}, + {.angle_measurement_type = local_angle, .measurement = current_measurement}}; + check_observable(topo, param, se_input); + } + SUBCASE("Global angle current sensor") { + se_input.measured_branch_from_current = { + {.angle_measurement_type = global_angle, .measurement = current_measurement}, + {.angle_measurement_type = global_angle, .measurement = current_measurement}}; + check_not_observable(topo, param, se_input); + } + } + } +} } // namespace power_grid_model diff --git a/tests/cpp_unit_tests/test_statistics.cpp b/tests/cpp_unit_tests/test_statistics.cpp index 6a8462f73a..4df368497b 100644 --- a/tests/cpp_unit_tests/test_statistics.cpp +++ b/tests/cpp_unit_tests/test_statistics.cpp @@ -815,6 +815,400 @@ TEST_CASE("Test statistics") { } } +TEST_CASE("Test statistics - combine") { + using statistics::combine; + using std::views::take; + + SUBCASE("UniformRealRandVar | IndependentRealRandVar") { + // using a template lambda to avoid code duplication and to avoid having to create a separate test case + auto const check = []() { + std::vector const measurements{ + {.value = 1.0, .variance = 0.2}, {.value = 2.0, .variance = 0.3}, {.value = 3.0, .variance = 0.6}}; + + CHECK(combine(measurements | take(0)).value == 0.0); + CHECK(is_inf(combine(measurements | take(0)).variance)); + + CHECK(combine(measurements | take(1)).value == measurements.front().value); + CHECK(combine(measurements | take(1)).variance == measurements.front().variance); + + CHECK(combine(measurements | take(2)).value == doctest::Approx(7.0 / 5.0)); + CHECK(combine(measurements | take(2)).variance == doctest::Approx(3.0 / 25.0)); + + CHECK(combine(measurements | take(3)).value == doctest::Approx(5.0 / 3.0)); + CHECK(combine(measurements | take(3)).variance == doctest::Approx(1.0 / 10.0)); + }; + SUBCASE("UniformRealRandVar") { check.template operator()>(); } + SUBCASE("IndependentRealRandVar") { + check.template operator()>(); + } + } + + SUBCASE("UniformRealRandVar") { + std::vector> const measurements{{.value = {1.0, 2.0, -1.0}, .variance = 0.2}, + {.value = {2.0, 4.0, 3.0}, .variance = 0.3}, + {.value = {4.0, 5.0, 6.0}, .variance = 0.6}}; + + CHECK(combine(measurements | take(0)).value(0) == 0.0); + CHECK(combine(measurements | take(0)).value(1) == 0.0); + CHECK(combine(measurements | take(0)).value(2) == 0.0); + CHECK(is_inf(combine(measurements | take(0)).variance)); + + CHECK(combine(measurements | take(1)).value(0) == measurements.front().value(0)); + CHECK(combine(measurements | take(1)).value(1) == measurements.front().value(1)); + CHECK(combine(measurements | take(1)).value(2) == measurements.front().value(2)); + CHECK(combine(measurements | take(1)).variance == measurements.front().variance); + + CHECK(combine(measurements | take(2)).value(0) == doctest::Approx(7.0 / 5.0)); + CHECK(combine(measurements | take(2)).value(1) == doctest::Approx(14.0 / 5.0)); + CHECK(combine(measurements | take(2)).value(2) == doctest::Approx(3.0 / 5.0)); + CHECK(combine(measurements | take(2)).variance == doctest::Approx(3.0 / 25.0)); + + CHECK(combine(measurements | take(3)).value(0) == doctest::Approx(11.0 / 6.0)); + CHECK(combine(measurements | take(3)).value(1) == doctest::Approx(19.0 / 6.0)); + CHECK(combine(measurements | take(3)).value(2) == doctest::Approx(3.0 / 2.0)); + CHECK(combine(measurements | take(3)).variance == doctest::Approx(1.0 / 10.0)); + } + + SUBCASE("IndependentRealRandVar") { + std::vector> const measurements{ + {.value = {1.0, 2.0, -1.0}, .variance = {0.2, 0.3, 0.4}}, + {.value = {2.0, 4.0, 3.0}, .variance = {0.3, 0.4, 0.5}}, + {.value = {4.0, 5.0, 6.0}, .variance = {0.6, 0.7, 0.8}}}; + + CHECK(combine(measurements | take(0)).value(0) == 0.0); + CHECK(combine(measurements | take(0)).value(1) == 0.0); + CHECK(combine(measurements | take(0)).value(2) == 0.0); + CHECK(is_inf(combine(measurements | take(0)).variance(0))); + CHECK(is_inf(combine(measurements | take(0)).variance(1))); + CHECK(is_inf(combine(measurements | take(0)).variance(2))); + + CHECK(combine(measurements | take(1)).value(0) == measurements.front().value(0)); + CHECK(combine(measurements | take(1)).value(1) == measurements.front().value(1)); + CHECK(combine(measurements | take(1)).value(2) == measurements.front().value(2)); + CHECK(combine(measurements | take(1)).variance(0) == measurements.front().variance(0)); + CHECK(combine(measurements | take(1)).variance(1) == measurements.front().variance(1)); + CHECK(combine(measurements | take(1)).variance(2) == measurements.front().variance(2)); + + CHECK(combine(measurements | take(2)).value(0) == doctest::Approx(7.0 / 5.0)); + CHECK(combine(measurements | take(2)).value(1) == doctest::Approx(20.0 / 7.0)); + CHECK(combine(measurements | take(2)).value(2) == doctest::Approx(7.0 / 9.0)); + CHECK(combine(measurements | take(2)).variance(0) == doctest::Approx(3.0 / 25.0)); + CHECK(combine(measurements | take(2)).variance(1) == doctest::Approx(6.0 / 35.0)); + CHECK(combine(measurements | take(2)).variance(2) == doctest::Approx(2.0 / 9.0)); + + CHECK(combine(measurements | take(3)).value(0) == doctest::Approx(11.0 / 6.0)); + CHECK(combine(measurements | take(3)).value(1) == doctest::Approx(200.0 / 61.0)); + CHECK(combine(measurements | take(3)).value(2) == doctest::Approx(44.0 / 23.0)); + CHECK(combine(measurements | take(3)).variance(0) == doctest::Approx(1.0 / 10.0)); + CHECK(combine(measurements | take(3)).variance(1) == doctest::Approx(42.0 / 305.0)); + CHECK(combine(measurements | take(3)).variance(2) == doctest::Approx(4.0 / 23.0)); + } + + SUBCASE("UniformComplexRandVar | IndependentComplexRandVar") { + // using a template lambda to avoid code duplication and to avoid having to create a separate test case + auto const check = []() { + std::vector const measurements{T{.value = 1.0 + 5.0i, .variance = 0.2}, + T{.value = 2.0 + 6.0i, .variance = 0.3}, + T{.value = 4.0 + 3.0i, .variance = 0.6}}; + + CHECK(combine(measurements | take(0)).value.real() == 0.0); + CHECK(combine(measurements | take(0)).value.imag() == 0.0); + CHECK(is_inf(combine(measurements | take(0)).variance)); + + CHECK(combine(measurements | take(1)).value.real() == measurements.front().value.real()); + CHECK(combine(measurements | take(1)).value.imag() == measurements.front().value.imag()); + CHECK(combine(measurements | take(1)).variance == measurements.front().variance); + + CHECK(combine(measurements | take(2)).value.real() == doctest::Approx(7.0 / 5.0)); + CHECK(combine(measurements | take(2)).value.imag() == doctest::Approx(27.0 / 5.0)); + CHECK(combine(measurements | take(2)).variance == doctest::Approx(3.0 / 25.0)); + + CHECK(combine(measurements | take(3)).value.real() == doctest::Approx(11.0 / 6.0)); + CHECK(combine(measurements | take(3)).value.imag() == doctest::Approx(30.0 / 6.0)); + CHECK(combine(measurements | take(3)).variance == doctest::Approx(1.0 / 10.0)); + }; + SUBCASE("UniformComplexRandVar") { + check.template operator()>(); + } + SUBCASE("IndependentComplexRandVar") { + check.template operator()>(); + } + } + + SUBCASE("UniformComplexRandVar") { + std::vector> const measurements{ + {.value = {RealValue{1.0, 2.0, -1.0}, RealValue{5.0, 6.0, 7.0}}, + .variance = 0.2}, + {.value = {RealValue{2.0, 4.0, 3.0}, RealValue{6.0, -7.0, 2.0}}, + .variance = 0.3}, + {.value = {RealValue{4.0, 5.0, 6.0}, RealValue{3.0, 1.0, 2.0}}, + .variance = 0.6}}; + + CHECK(combine(measurements | take(0)).value(0).real() == 0.0); + CHECK(combine(measurements | take(0)).value(1).real() == 0.0); + CHECK(combine(measurements | take(0)).value(2).real() == 0.0); + CHECK(combine(measurements | take(0)).value(0).imag() == 0.0); + CHECK(combine(measurements | take(0)).value(1).imag() == 0.0); + CHECK(combine(measurements | take(0)).value(2).imag() == 0.0); + CHECK(is_inf(combine(measurements | take(0)).variance)); + + CHECK(combine(measurements | take(1)).value(0).real() == measurements.front().value(0).real()); + CHECK(combine(measurements | take(1)).value(1).real() == measurements.front().value(1).real()); + CHECK(combine(measurements | take(1)).value(2).real() == measurements.front().value(2).real()); + CHECK(combine(measurements | take(1)).value(0).imag() == measurements.front().value(0).imag()); + CHECK(combine(measurements | take(1)).value(1).imag() == measurements.front().value(1).imag()); + CHECK(combine(measurements | take(1)).value(2).imag() == measurements.front().value(2).imag()); + CHECK(combine(measurements | take(1)).variance == measurements.front().variance); + + CHECK(combine(measurements | take(2)).value(0).real() == doctest::Approx(7.0 / 5.0)); + CHECK(combine(measurements | take(2)).value(1).real() == doctest::Approx(14.0 / 5.0)); + CHECK(combine(measurements | take(2)).value(2).real() == doctest::Approx(3.0 / 5.0)); + CHECK(combine(measurements | take(2)).value(0).imag() == doctest::Approx(27.0 / 5.0)); + CHECK(combine(measurements | take(2)).value(1).imag() == doctest::Approx(4.0 / 5.0)); + CHECK(combine(measurements | take(2)).value(2).imag() == doctest::Approx(25.0 / 5.0)); + CHECK(combine(measurements | take(2)).variance == doctest::Approx(3.0 / 25.0)); + + CHECK(combine(measurements | take(3)).value(0).real() == doctest::Approx(11.0 / 6.0)); + CHECK(combine(measurements | take(3)).value(1).real() == doctest::Approx(19.0 / 6.0)); + CHECK(combine(measurements | take(3)).value(2).real() == doctest::Approx(9.0 / 6.0)); + CHECK(combine(measurements | take(3)).value(0).imag() == doctest::Approx(30.0 / 6.0)); + CHECK(combine(measurements | take(3)).value(1).imag() == doctest::Approx(5.0 / 6.0)); + CHECK(combine(measurements | take(3)).value(2).imag() == doctest::Approx(27.0 / 6.0)); + CHECK(combine(measurements | take(3)).variance == doctest::Approx(1.0 / 10.0)); + } + + SUBCASE("IndependentComplexRandVar") { + std::vector> const measurements{ + {.value = {RealValue{1.0, 2.0, -1.0}, RealValue{5.0, 6.0, 7.0}}, + .variance = {0.2, 0.3, 0.4}}, + {.value = {RealValue{2.0, 4.0, 3.0}, RealValue{6.0, -7.0, 2.0}}, + .variance = {0.3, 0.4, 0.5}}, + {.value = {RealValue{4.0, 5.0, 6.0}, RealValue{3.0, 1.0, 2.0}}, + .variance = {0.6, 0.7, 0.8}}}; + + CHECK(combine(measurements | take(0)).value(0).real() == 0.0); + CHECK(combine(measurements | take(0)).value(1).real() == 0.0); + CHECK(combine(measurements | take(0)).value(2).real() == 0.0); + CHECK(combine(measurements | take(0)).value(0).imag() == 0.0); + CHECK(combine(measurements | take(0)).value(1).imag() == 0.0); + CHECK(combine(measurements | take(0)).value(2).imag() == 0.0); + CHECK(is_inf(combine(measurements | take(0)).variance(0))); + CHECK(is_inf(combine(measurements | take(0)).variance(1))); + CHECK(is_inf(combine(measurements | take(0)).variance(2))); + + CHECK(combine(measurements | take(1)).value(0).real() == measurements.front().value(0).real()); + CHECK(combine(measurements | take(1)).value(1).real() == measurements.front().value(1).real()); + CHECK(combine(measurements | take(1)).value(2).real() == measurements.front().value(2).real()); + CHECK(combine(measurements | take(1)).value(0).imag() == measurements.front().value(0).imag()); + CHECK(combine(measurements | take(1)).value(1).imag() == measurements.front().value(1).imag()); + CHECK(combine(measurements | take(1)).value(2).imag() == measurements.front().value(2).imag()); + CHECK(combine(measurements | take(1)).variance(0) == measurements.front().variance(0)); + CHECK(combine(measurements | take(1)).variance(1) == measurements.front().variance(1)); + CHECK(combine(measurements | take(1)).variance(2) == measurements.front().variance(2)); + + CHECK(combine(measurements | take(2)).value(0).real() == doctest::Approx(7.0 / 5.0)); + CHECK(combine(measurements | take(2)).value(1).real() == doctest::Approx(20.0 / 7.0)); + CHECK(combine(measurements | take(2)).value(2).real() == doctest::Approx(7.0 / 9.0)); + CHECK(combine(measurements | take(2)).value(0).imag() == doctest::Approx(27.0 / 5.0)); + CHECK(combine(measurements | take(2)).value(1).imag() == doctest::Approx(3.0 / 7.0)); + CHECK(combine(measurements | take(2)).value(2).imag() == doctest::Approx(43.0 / 9.0)); + CHECK(combine(measurements | take(2)).variance(0) == doctest::Approx(3.0 / 25.0)); + CHECK(combine(measurements | take(2)).variance(1) == doctest::Approx(6.0 / 35.0)); + CHECK(combine(measurements | take(2)).variance(2) == doctest::Approx(2.0 / 9.0)); + + CHECK(combine(measurements | take(3)).value(0).real() == doctest::Approx(11.0 / 6.0)); + CHECK(combine(measurements | take(3)).value(1).real() == doctest::Approx(200.0 / 61.0)); + CHECK(combine(measurements | take(3)).value(2).real() == doctest::Approx(44.0 / 23.0)); + CHECK(combine(measurements | take(3)).value(0).imag() == doctest::Approx(30.0 / 6.0)); + CHECK(combine(measurements | take(3)).value(1).imag() == doctest::Approx(33.0 / 61.0)); + CHECK(combine(measurements | take(3)).value(2).imag() == doctest::Approx(96.0 / 23.0)); + CHECK(combine(measurements | take(3)).variance(0) == doctest::Approx(1.0 / 10.0)); + CHECK(combine(measurements | take(3)).variance(1) == doctest::Approx(42.0 / 305.0)); + CHECK(combine(measurements | take(3)).variance(2) == doctest::Approx(4.0 / 23.0)); + } + + SUBCASE("DecomposedComplexRandVar") { + std::vector> const measurements{ + DecomposedComplexRandVar{.real_component = {.value = 1.0, .variance = 0.2}, + .imag_component = {.value = 5.0, .variance = 0.1}}, + DecomposedComplexRandVar{.real_component = {.value = 2.0, .variance = 0.3}, + .imag_component = {.value = 6.0, .variance = 0.2}}, + DecomposedComplexRandVar{.real_component = {.value = 4.0, .variance = 0.6}, + .imag_component = {.value = 3.0, .variance = 0.3}}}; + + CHECK(combine(measurements | take(0)).real_component.value == 0.0); + CHECK(combine(measurements | take(0)).imag_component.value == 0.0); + CHECK(is_inf(combine(measurements | take(0)).real_component.variance)); + CHECK(is_inf(combine(measurements | take(0)).imag_component.variance)); + + CHECK(combine(measurements | take(1)).real_component.value == measurements.front().real_component.value); + CHECK(combine(measurements | take(1)).imag_component.value == measurements.front().imag_component.value); + CHECK(combine(measurements | take(1)).real_component.variance == measurements.front().real_component.variance); + CHECK(combine(measurements | take(1)).imag_component.variance == measurements.front().imag_component.variance); + + CHECK(combine(measurements | take(2)).real_component.value == doctest::Approx(7.0 / 5.0)); + CHECK(combine(measurements | take(2)).imag_component.value == doctest::Approx(80.0 / 15.0)); + CHECK(combine(measurements | take(2)).real_component.variance == doctest::Approx(3.0 / 25.0)); + CHECK(combine(measurements | take(2)).imag_component.variance == doctest::Approx(1.0 / 15.0)); + + CHECK(combine(measurements | take(3)).real_component.value == doctest::Approx(11.0 / 6.0)); + CHECK(combine(measurements | take(3)).imag_component.value == doctest::Approx(270.0 / 55.0)); + CHECK(combine(measurements | take(3)).real_component.variance == doctest::Approx(1.0 / 10.0)); + CHECK(combine(measurements | take(3)).imag_component.variance == doctest::Approx(3.0 / 55.0)); + } + + SUBCASE("DecomposedComplexRandVar") { + std::vector> const measurements{ + DecomposedComplexRandVar{ + .real_component = {.value = RealValue{1.0, 2.0, -1.0}, .variance = {0.2, 0.3, 0.4}}, + .imag_component = {.value = RealValue{5.0, 6.0, 7.0}, .variance = {0.1, 0.2, 0.3}}}, + DecomposedComplexRandVar{ + .real_component = {.value = RealValue{2.0, 4.0, 3.0}, .variance = {0.3, 0.4, 0.5}}, + .imag_component = {.value = RealValue{6.0, -7.0, 2.0}, .variance = {0.2, 0.3, 0.4}}}, + DecomposedComplexRandVar{ + .real_component = {.value = RealValue{4.0, 5.0, 6.0}, .variance = {0.6, 0.7, 0.8}}, + .imag_component = {.value = RealValue{3.0, 1.0, 2.0}, .variance = {0.3, 0.4, 0.5}}}}; + + CHECK(combine(measurements | take(0)).real_component.value(0) == 0.0); + CHECK(combine(measurements | take(0)).real_component.value(2) == 0.0); + CHECK(combine(measurements | take(0)).real_component.value(1) == 0.0); + CHECK(combine(measurements | take(0)).imag_component.value(0) == 0.0); + CHECK(combine(measurements | take(0)).imag_component.value(2) == 0.0); + CHECK(combine(measurements | take(0)).imag_component.value(1) == 0.0); + CHECK(is_inf(combine(measurements | take(0)).real_component.variance(0))); + CHECK(is_inf(combine(measurements | take(0)).real_component.variance(1))); + CHECK(is_inf(combine(measurements | take(0)).real_component.variance(2))); + CHECK(is_inf(combine(measurements | take(0)).imag_component.variance(0))); + CHECK(is_inf(combine(measurements | take(0)).imag_component.variance(1))); + CHECK(is_inf(combine(measurements | take(0)).imag_component.variance(2))); + + CHECK(combine(measurements | take(1)).real_component.value(0) == measurements.front().real_component.value(0)); + CHECK(combine(measurements | take(1)).real_component.value(1) == measurements.front().real_component.value(1)); + CHECK(combine(measurements | take(1)).real_component.value(2) == measurements.front().real_component.value(2)); + CHECK(combine(measurements | take(1)).imag_component.value(0) == measurements.front().imag_component.value(0)); + CHECK(combine(measurements | take(1)).imag_component.value(1) == measurements.front().imag_component.value(1)); + CHECK(combine(measurements | take(1)).imag_component.value(2) == measurements.front().imag_component.value(2)); + CHECK(combine(measurements | take(1)).real_component.variance(0) == + measurements.front().real_component.variance(0)); + CHECK(combine(measurements | take(1)).real_component.variance(1) == + measurements.front().real_component.variance(1)); + CHECK(combine(measurements | take(1)).real_component.variance(2) == + measurements.front().real_component.variance(2)); + CHECK(combine(measurements | take(1)).imag_component.variance(0) == + measurements.front().imag_component.variance(0)); + CHECK(combine(measurements | take(1)).imag_component.variance(1) == + measurements.front().imag_component.variance(1)); + CHECK(combine(measurements | take(1)).imag_component.variance(2) == + measurements.front().imag_component.variance(2)); + + CHECK(combine(measurements | take(2)).real_component.value(0) == doctest::Approx(7.0 / 5.0)); + CHECK(combine(measurements | take(2)).real_component.value(1) == doctest::Approx(20.0 / 7.0)); + CHECK(combine(measurements | take(2)).real_component.value(2) == doctest::Approx(7.0 / 9.0)); + CHECK(combine(measurements | take(2)).imag_component.value(0) == doctest::Approx(80.0 / 15.0)); + CHECK(combine(measurements | take(2)).imag_component.value(1) == doctest::Approx(20 / 25.0)); + CHECK(combine(measurements | take(2)).imag_component.value(2) == doctest::Approx(170.0 / 35.0)); + CHECK(combine(measurements | take(2)).real_component.variance(0) == doctest::Approx(3.0 / 25.0)); + CHECK(combine(measurements | take(2)).real_component.variance(1) == doctest::Approx(6.0 / 35.0)); + CHECK(combine(measurements | take(2)).real_component.variance(2) == doctest::Approx(2.0 / 9.0)); + CHECK(combine(measurements | take(2)).imag_component.variance(0) == doctest::Approx(1.0 / 15.0)); + CHECK(combine(measurements | take(2)).imag_component.variance(1) == doctest::Approx(3.0 / 25.0)); + CHECK(combine(measurements | take(2)).imag_component.variance(2) == doctest::Approx(6.0 / 35.0)); + + CHECK(combine(measurements | take(3)).real_component.value(0) == doctest::Approx(11.0 / 6.0)); + CHECK(combine(measurements | take(3)).real_component.value(1) == doctest::Approx(200.0 / 61.0)); + CHECK(combine(measurements | take(3)).real_component.value(2) == doctest::Approx(44.0 / 23.0)); + CHECK(combine(measurements | take(3)).imag_component.value(0) == doctest::Approx(270.0 / 55.0)); + CHECK(combine(measurements | take(3)).imag_component.value(1) == doctest::Approx(55.0 / 65.0)); + CHECK(combine(measurements | take(3)).imag_component.value(2) == doctest::Approx(194.0 / 47.0)); + CHECK(combine(measurements | take(3)).real_component.variance(0) == doctest::Approx(1.0 / 10.0)); + CHECK(combine(measurements | take(3)).real_component.variance(1) == doctest::Approx(42.0 / 305.0)); + CHECK(combine(measurements | take(3)).real_component.variance(2) == doctest::Approx(4.0 / 23.0)); + CHECK(combine(measurements | take(3)).imag_component.variance(0) == doctest::Approx(3.0 / 55.0)); + CHECK(combine(measurements | take(3)).imag_component.variance(1) == doctest::Approx(6.0 / 65)); + CHECK(combine(measurements | take(3)).imag_component.variance(2) == doctest::Approx(6.0 / 47.0)); + } +} + +TEST_CASE("Test statistics - combine_magnitude") { + using statistics::combine_magnitude; + using std::views::take; + + SUBCASE("UniformComplexRandVar") { + // using a template lambda to avoid code duplication and to avoid having to create a separate test case + std::vector> const measurements{ + {.value = ComplexValue{1.0, 5.0}, .variance = 0.2}, + {.value = ComplexValue{2.0, nan}, .variance = 0.3}, + {.value = ComplexValue{4.0, nan}, .variance = 0.6}}; + + CHECK(combine_magnitude(measurements | take(0)).value.real() == 0.0); + CHECK(is_nan(combine_magnitude(measurements | take(0)).value.imag())); + CHECK(is_inf(combine_magnitude(measurements | take(0)).variance)); + + CHECK(combine_magnitude(measurements | take(1)).value.real() == cabs(measurements.front().value)); + CHECK(is_nan(combine_magnitude(measurements | take(1)).value.imag())); + CHECK(combine_magnitude(measurements | take(1)).variance == measurements.front().variance); + + CHECK(combine_magnitude(measurements | take(2)).value.real() == + doctest::Approx((3.0 * std::sqrt(26.0) + 4.0) / 5.0)); + CHECK(is_nan(combine_magnitude(measurements | take(2)).value.imag())); + CHECK(combine_magnitude(measurements | take(2)).variance == doctest::Approx(3.0 / 25.0)); + + CHECK(combine_magnitude(measurements | take(3)).value.real() == + doctest::Approx((8.0 + 3.0 * std::sqrt(26.0)) / 6.0)); + CHECK(is_nan(combine_magnitude(measurements | take(3)).value.imag())); + CHECK(combine_magnitude(measurements | take(3)).variance == doctest::Approx(1.0 / 10.0)); + } + + SUBCASE("UniformComplexRandVar") { + std::vector> const measurements{ + {.value = {RealValue{1.0, 2.0, -1.0}, RealValue{5.0, 6.0, 7.0}}, + .variance = 0.2}, + {.value = {RealValue{2.0, 4.0, 3.0}, RealValue{nan}}, .variance = 0.3}, + {.value = {RealValue{4.0, 5.0, 6.0}, RealValue{nan}}, .variance = 0.6}}; + + CHECK(combine_magnitude(measurements | take(0)).value(0).real() == 0.0); + CHECK(combine_magnitude(measurements | take(0)).value(1).real() == 0.0); + CHECK(combine_magnitude(measurements | take(0)).value(2).real() == 0.0); + CHECK(is_nan(combine_magnitude(measurements | take(0)).value(0).imag())); + CHECK(is_nan(combine_magnitude(measurements | take(0)).value(0).imag())); + CHECK(is_nan(combine_magnitude(measurements | take(0)).value(0).imag())); + CHECK(is_inf(combine_magnitude(measurements | take(0)).variance)); + + CHECK(combine_magnitude(measurements | take(1)).value(0).real() == + doctest::Approx(cabs(measurements.front().value(0)))); + CHECK(combine_magnitude(measurements | take(1)).value(1).real() == + doctest::Approx(cabs(measurements.front().value(1)))); + CHECK(combine_magnitude(measurements | take(1)).value(2).real() == + doctest::Approx(cabs(measurements.front().value(2)))); + CHECK(is_nan(combine_magnitude(measurements | take(1)).value(0).imag())); + CHECK(is_nan(combine_magnitude(measurements | take(1)).value(1).imag())); + CHECK(is_nan(combine_magnitude(measurements | take(1)).value(2).imag())); + CHECK(combine_magnitude(measurements | take(1)).variance == measurements.front().variance); + + CHECK(combine_magnitude(measurements | take(2)).value(0).real() == + doctest::Approx((3.0 * std::sqrt(26.0) + 4.0) / 5.0)); + CHECK(combine_magnitude(measurements | take(2)).value(1).real() == + doctest::Approx((8.0 + 6.0 * std::sqrt(10.0)) / 5.0)); + CHECK(combine_magnitude(measurements | take(2)).value(2).real() == + doctest::Approx((6.0 + 15.0 * std::sqrt(2.0)) / 5.0)); + CHECK(is_nan(combine_magnitude(measurements | take(2)).value(0).imag())); + CHECK(is_nan(combine_magnitude(measurements | take(2)).value(1).imag())); + CHECK(is_nan(combine_magnitude(measurements | take(2)).value(2).imag())); + CHECK(combine_magnitude(measurements | take(2)).variance == doctest::Approx(3.0 / 25.0)); + + CHECK(combine_magnitude(measurements | take(3)).value(0).real() == + doctest::Approx((8.0 + 3.0 * std::sqrt(26.0)) / 6.0)); + CHECK(combine_magnitude(measurements | take(3)).value(1).real() == + doctest::Approx((13.0 + 6.0 * std::sqrt(10.0)) / 6.0)); + CHECK(combine_magnitude(measurements | take(3)).value(2).real() == + doctest::Approx((4.0 + 5.0 * std::sqrt(2.0)) / 2.0)); + CHECK(is_nan(combine_magnitude(measurements | take(3)).value(0).imag())); + CHECK(is_nan(combine_magnitude(measurements | take(3)).value(1).imag())); + CHECK(is_nan(combine_magnitude(measurements | take(3)).value(2).imag())); + CHECK(combine_magnitude(measurements | take(3)).variance == doctest::Approx(1.0 / 10.0)); + } +} + TEST_SUITE_END(); } // namespace power_grid_model diff --git a/tests/data/state_estimation/basic-current-sensor/params.json b/tests/data/state_estimation/basic-current-sensor/params.json index b8769a93ef..1646fba785 100644 --- a/tests/data/state_estimation/basic-current-sensor/params.json +++ b/tests/data/state_estimation/basic-current-sensor/params.json @@ -12,14 +12,14 @@ "iterative_linear": { "experimental_features": "enabled", "fail": { - "raises": "NotObservableError", + "raises": "SparseMatrixError", "reason": "Current sensors are not yet implemented for this calculation method" } }, "newton_raphson": { "experimental_features": "enabled", "fail": { - "raises": "NotObservableError", + "raises": "SparseMatrixError", "reason": "Current sensors are not yet implemented for this calculation method" } } diff --git a/tests/unit/test_error_handling.py b/tests/unit/test_error_handling.py index b322757e5a..7e53009020 100644 --- a/tests/unit/test_error_handling.py +++ b/tests/unit/test_error_handling.py @@ -8,10 +8,18 @@ from power_grid_model import PowerGridModel from power_grid_model._core.power_grid_meta import initialize_array -from power_grid_model.enum import CalculationMethod, LoadGenType, MeasuredTerminalType, TapChangingStrategy +from power_grid_model.enum import ( + AngleMeasurementType, + CalculationMethod, + LoadGenType, + MeasuredTerminalType, + TapChangingStrategy, + _ExperimentalFeatures, +) from power_grid_model.errors import ( AutomaticTapInputError, ConflictID, + ConflictingAngleMeasurementType, ConflictVoltage, IDNotFound, IDWrongType, @@ -349,6 +357,107 @@ def test_automatic_tap_changing(): model.calculate_power_flow(tap_changing_strategy=TapChangingStrategy.min_voltage_tap) +def test_conflicting_angle_measurement_type() -> None: + node_input = initialize_array("input", "node", 2) + node_input["id"] = [0, 1] + node_input["u_rated"] = [10e3, 10e3] + + line_input = initialize_array("input", "line", 1) + line_input["id"] = [2] + line_input["from_node"] = [0] + line_input["to_node"] = [1] + + source_input = initialize_array("input", "source", 1) + source_input["id"] = [3] + source_input["node"] = [0] + source_input["status"] = [1] + source_input["u_ref"] = [1.0] + + sym_voltage_sensor_input = initialize_array("input", "sym_voltage_sensor", 1) + sym_voltage_sensor_input["id"] = [4] + sym_voltage_sensor_input["measured_object"] = [0] + sym_voltage_sensor_input["u_sigma"] = [1.0] + sym_voltage_sensor_input["u_measured"] = [1.0] + sym_voltage_sensor_input["u_angle_measured"] = [0.0] + + sym_current_sensor_input = initialize_array("input", "sym_current_sensor", 2) + sym_current_sensor_input["id"] = [5, 6] + sym_current_sensor_input["measured_object"] = [2, 2] + sym_current_sensor_input["measured_terminal_type"] = [ + MeasuredTerminalType.branch_from, + MeasuredTerminalType.branch_from, + ] + sym_current_sensor_input["angle_measurement_type"] = [ + AngleMeasurementType.local_angle, + AngleMeasurementType.global_angle, + ] + sym_current_sensor_input["i_sigma"] = [1.0, 1.0] + sym_current_sensor_input["i_angle_sigma"] = [1.0, 1.0] + sym_current_sensor_input["i_measured"] = [0.0, 0.0] + sym_current_sensor_input["i_angle_measured"] = [0.0, 0.0] + + model = PowerGridModel( + input_data={ + "node": node_input, + "line": line_input, + "source": source_input, + "sym_voltage_sensor": sym_voltage_sensor_input, + "sym_current_sensor": sym_current_sensor_input, + } + ) + + with pytest.raises(ConflictingAngleMeasurementType): + model._calculate_state_estimation(decode_error=True, experimental_features=_ExperimentalFeatures.enabled) + + +def test_global_current_measurement_without_voltage_angle() -> None: + node_input = initialize_array("input", "node", 2) + node_input["id"] = [0, 1] + node_input["u_rated"] = [10e3, 10e3] + + line_input = initialize_array("input", "line", 1) + line_input["id"] = [2] + line_input["from_node"] = [0] + line_input["to_node"] = [1] + + source_input = initialize_array("input", "source", 1) + source_input["id"] = [3] + source_input["node"] = [0] + source_input["status"] = [1] + source_input["u_ref"] = [1.0] + + sym_voltage_sensor_input = initialize_array("input", "sym_voltage_sensor", 1) + sym_voltage_sensor_input["id"] = [4] + sym_voltage_sensor_input["measured_object"] = [0] + sym_voltage_sensor_input["u_sigma"] = [1.0] + sym_voltage_sensor_input["u_measured"] = [0.0] + + sym_current_sensor_input = initialize_array("input", "sym_current_sensor", 1) + sym_current_sensor_input["id"] = [5] + sym_current_sensor_input["measured_object"] = [2] + sym_current_sensor_input["measured_terminal_type"] = [ + MeasuredTerminalType.branch_from, + ] + sym_current_sensor_input["angle_measurement_type"] = [ + AngleMeasurementType.global_angle, + ] + sym_current_sensor_input["i_sigma"] = [1.0] + sym_current_sensor_input["i_measured"] = [0.0] + + model = PowerGridModel( + input_data={ + "node": node_input, + "line": line_input, + "source": source_input, + "sym_voltage_sensor": sym_voltage_sensor_input, + "sym_current_sensor": sym_current_sensor_input, + } + ) + + with pytest.raises(NotObservableError): + model._calculate_state_estimation(decode_error=True, experimental_features=_ExperimentalFeatures.enabled) + + @pytest.mark.skip(reason="TODO") def test_handle_power_grid_dataset_error(): pass