Merge pull request #961 from PowerGridModel/fix-current-sensor-output-residuals

Fix output residuals for current sensor

Author: mgovers

power_grid_model_c/power_grid_model/include/power_grid_model/all_components.hpp

@@ -31,16 +31,4 @@ using AllComponents =
     ComponentList<Node, Line, Link, GenericBranch, Transformer, ThreeWindingTransformer, Shunt, Source, SymGenerator,
                   AsymGenerator, SymLoad, AsymLoad, SymPowerSensor, AsymPowerSensor, SymVoltageSensor,
                   AsymVoltageSensor, SymCurrentSensor, AsymCurrentSensor, Fault, TransformerTapRegulator>;
-
-template <typename T>
-concept power_or_current_sensor_c =
-    std::derived_from<T, GenericPowerSensor> || std::derived_from<T, GenericCurrentSensor>;
-
-static_assert(power_or_current_sensor_c<SymPowerSensor>);
-static_assert(power_or_current_sensor_c<AsymPowerSensor>);
-static_assert(power_or_current_sensor_c<SymCurrentSensor>);
-static_assert(power_or_current_sensor_c<AsymCurrentSensor>);
-static_assert(!power_or_current_sensor_c<SymVoltageSensor>);
-static_assert(!power_or_current_sensor_c<AsymVoltageSensor>);
-
 } // namespace power_grid_model

power_grid_model_c/power_grid_model/include/power_grid_model/component/current_sensor.hpp

@@ -26,11 +26,12 @@ class GenericCurrentSensor : public Sensor {
     MeasuredTerminalType get_terminal_type() const { return terminal_type_; }
     AngleMeasurementType get_angle_measurement_type() const { return angle_measurement_type_; }
 
-    template <symmetry_tag sym> CurrentSensorOutput<sym> get_output(ComplexValue<sym> const& i) const {
+    template <symmetry_tag sym>
+    CurrentSensorOutput<sym> get_output(ComplexValue<sym> const& i, ComplexValue<sym> const& u) const {
         if constexpr (is_symmetric_v<sym>) {
-            return get_sym_output(i);
+            return get_sym_output(i, u);
         } else {
-            return get_asym_output(i);
+            return get_asym_output(i, u);
         }
     }
 
@@ -62,8 +63,10 @@ class GenericCurrentSensor : public Sensor {
     virtual CurrentSensorCalcParam<symmetric_t> sym_calc_param() const = 0;
     virtual CurrentSensorCalcParam<asymmetric_t> asym_calc_param() const = 0;
 
-    virtual CurrentSensorOutput<symmetric_t> get_sym_output(ComplexValue<symmetric_t> const& i) const = 0;
-    virtual CurrentSensorOutput<asymmetric_t> get_asym_output(ComplexValue<asymmetric_t> const& i) const = 0;
+    virtual CurrentSensorOutput<symmetric_t> get_sym_output(ComplexValue<symmetric_t> const& i,
+                                                            ComplexValue<symmetric_t> const& u) const = 0;
+    virtual CurrentSensorOutput<asymmetric_t> get_asym_output(ComplexValue<asymmetric_t> const& i,
+                                                              ComplexValue<asymmetric_t> const& u) const = 0;
 };
 
 template <symmetry_tag current_sensor_symmetry_> class CurrentSensor : public GenericCurrentSensor {
@@ -137,20 +140,43 @@ template <symmetry_tag current_sensor_symmetry_> class CurrentSensor : public Ge
         return CurrentSensorCalcParam<sym_calc>{.angle_measurement_type = angle_measurement_type(),
                                                 .measurement = DecomposedComplexRandVar<sym_calc>(i_polar)};
     }
-    CurrentSensorOutput<symmetric_t> get_sym_output(ComplexValue<symmetric_t> const& i) const final {
-        return get_generic_output<symmetric_t>(i);
+    CurrentSensorOutput<symmetric_t> get_sym_output(ComplexValue<symmetric_t> const& i,
+                                                    ComplexValue<symmetric_t> const& u) const final {
+        return get_generic_output<symmetric_t>(i, u);
     }
-    CurrentSensorOutput<asymmetric_t> get_asym_output(ComplexValue<asymmetric_t> const& i) const final {
-        return get_generic_output<asymmetric_t>(i);
+    CurrentSensorOutput<asymmetric_t> get_asym_output(ComplexValue<asymmetric_t> const& i,
+                                                      ComplexValue<asymmetric_t> const& u) const final {
+        return get_generic_output<asymmetric_t>(i, u);
     }
     template <symmetry_tag sym_calc>
-    CurrentSensorOutput<sym_calc> get_generic_output(ComplexValue<sym_calc> const& i) const {
+    CurrentSensorOutput<sym_calc> get_generic_output(ComplexValue<sym_calc> const& i,
+                                                     ComplexValue<sym_calc> const& u) const {
         CurrentSensorOutput<sym_calc> output{};
         output.id = id();
-        ComplexValue<sym_calc> const i_residual{process_mean_val<sym_calc>(i_measured_ - i) * base_current_};
         output.energized = 1; // current sensor is always energized
-        output.i_residual = cabs(i_residual);
-        output.i_angle_residual = arg(i_residual);
+        auto const i_calc_param = calc_param<sym_calc>();
+        auto const angle_measurement_type = i_calc_param.angle_measurement_type;
+        auto const i_measured_complex = i_calc_param.measurement.value();
+        ComplexValue<sym_calc> i_output = [&i, &u, &angle_measurement_type] {
+            switch (angle_measurement_type) {
+            case AngleMeasurementType::global_angle: {
+                return i;
+            }
+            case AngleMeasurementType::local_angle: {
+                // I_l = conj(I_g) * exp(i * u_angle)
+                // Tranform back the output angle to the local angle frame of reference
+                auto const u_angle = arg(u);
+                return ComplexValue<sym_calc>{conj(i) * (cos(u_angle) + 1i * sin(u_angle))};
+            }
+            default: {
+                throw MissingCaseForEnumError{"generic output angle measurement type", angle_measurement_type};
+            }
+            }
+        }();
+        output.i_residual = (cabs(i_measured_complex) - cabs(i_output)) * base_current_;
+        // arg(e^i * u_angle) = u_angle in (-pi, pi]
+        auto const unbounded_i_angle_residual = arg(i_measured_complex) - arg(i_output);
+        output.i_angle_residual = arg(ComplexValue<sym_calc>{exp(1.0i * unbounded_i_angle_residual)});
         return output;
     }
 };

power_grid_model_c/power_grid_model/include/power_grid_model/main_core/output.hpp

@@ -242,14 +242,14 @@ inline auto output_result(Component const& voltage_sensor, MainModelState<Compon
 }
 
 // output power sensor
-template <power_or_current_sensor_c Component, class ComponentContainer,
+template <std::derived_from<GenericPowerSensor> Component, class ComponentContainer,
           steady_state_solver_output_type SolverOutputType>
     requires model_component_state_c<MainModelState, ComponentContainer, Component>
-constexpr auto output_result(Component const& power_or_current_sensor, MainModelState<ComponentContainer> const& state,
+constexpr auto output_result(Component const& power_sensor, MainModelState<ComponentContainer> const& state,
                              std::vector<SolverOutputType> const& solver_output, Idx const obj_seq) {
     using sym = typename SolverOutputType::sym;
 
-    auto const terminal_type = power_or_current_sensor.get_terminal_type();
+    auto const terminal_type = power_sensor.get_terminal_type();
     Idx2D const obj_math_id = [&]() {
         switch (terminal_type) {
             using enum MeasuredTerminalType;
@@ -281,7 +281,7 @@ constexpr auto output_result(Component const& power_or_current_sensor, MainModel
     }();
 
     if (obj_math_id.group == -1) {
-        return power_or_current_sensor.template get_null_output<sym>();
+        return power_sensor.template get_null_output<sym>();
     }
 
     switch (terminal_type) {
@@ -295,30 +295,91 @@ constexpr auto output_result(Component const& power_or_current_sensor, MainModel
     case branch3_2:
         [[fallthrough]];
     case branch3_3:
-        return power_or_current_sensor.template get_output<sym>(
-            solver_output[obj_math_id.group].branch[obj_math_id.pos].s_f);
+        return power_sensor.template get_output<sym>(solver_output[obj_math_id.group].branch[obj_math_id.pos].s_f);
     case branch_to:
-        return power_or_current_sensor.template get_output<sym>(
-            solver_output[obj_math_id.group].branch[obj_math_id.pos].s_t);
+        return power_sensor.template get_output<sym>(solver_output[obj_math_id.group].branch[obj_math_id.pos].s_t);
     case source:
-        return power_or_current_sensor.template get_output<sym>(
-            solver_output[obj_math_id.group].source[obj_math_id.pos].s);
+        return power_sensor.template get_output<sym>(solver_output[obj_math_id.group].source[obj_math_id.pos].s);
     case shunt:
-        return power_or_current_sensor.template get_output<sym>(
-            solver_output[obj_math_id.group].shunt[obj_math_id.pos].s);
+        return power_sensor.template get_output<sym>(solver_output[obj_math_id.group].shunt[obj_math_id.pos].s);
     case load:
         [[fallthrough]];
     case generator:
-        return power_or_current_sensor.template get_output<sym>(
-            solver_output[obj_math_id.group].load_gen[obj_math_id.pos].s);
+        return power_sensor.template get_output<sym>(solver_output[obj_math_id.group].load_gen[obj_math_id.pos].s);
     case node:
-        return power_or_current_sensor.template get_output<sym>(
-            solver_output[obj_math_id.group].bus_injection[obj_math_id.pos]);
+        return power_sensor.template get_output<sym>(solver_output[obj_math_id.group].bus_injection[obj_math_id.pos]);
     default:
         throw MissingCaseForEnumError{std::format("{} output_result()", Component::name), terminal_type};
     }
 }
-template <power_or_current_sensor_c Component, class ComponentContainer,
+template <std::derived_from<GenericPowerSensor> Component, class ComponentContainer,
+          short_circuit_solver_output_type SolverOutputType>
+    requires model_component_state_c<MainModelState, ComponentContainer, Component>
+constexpr auto output_result(Component const& power_or_current_sensor,
+                             MainModelState<ComponentContainer> const& /* state */,
+                             std::vector<SolverOutputType> const& /* solver_output */, Idx const /* obj_seq */) {
+    return power_or_current_sensor.get_null_sc_output();
+}
+
+// output current sensor
+template <std::derived_from<GenericCurrentSensor> Component, class ComponentContainer,
+          steady_state_solver_output_type SolverOutputType>
+    requires model_component_state_c<MainModelState, ComponentContainer, Component>
+constexpr auto output_result(Component const& current_sensor, MainModelState<ComponentContainer> const& state,
+                             std::vector<SolverOutputType> const& solver_output, Idx const obj_seq) {
+    using sym = typename SolverOutputType::sym;
+
+    auto const terminal_type = current_sensor.get_terminal_type();
+    Idx2D const obj_math_id = [&]() {
+        switch (terminal_type) {
+            using enum MeasuredTerminalType;
+
+        case branch_from:
+            [[fallthrough]];
+        case branch_to:
+            return state.topo_comp_coup->branch[obj_seq];
+        // from branch3, get relevant math object branch based on the measured side
+        case branch3_1:
+            return Idx2D{state.topo_comp_coup->branch3[obj_seq].group, state.topo_comp_coup->branch3[obj_seq].pos[0]};
+        case branch3_2:
+            return Idx2D{state.topo_comp_coup->branch3[obj_seq].group, state.topo_comp_coup->branch3[obj_seq].pos[1]};
+        case branch3_3:
+            return Idx2D{state.topo_comp_coup->branch3[obj_seq].group, state.topo_comp_coup->branch3[obj_seq].pos[2]};
+        default:
+            throw MissingCaseForEnumError{std::format("{} output_result()", Component::name), terminal_type};
+        }
+    }();
+
+    if (obj_math_id.group == -1) {
+        return current_sensor.template get_null_output<sym>();
+    }
+
+    auto const topological_index = get_topology_index<Branch>(state, obj_math_id);
+    auto const branch_nodes = get_branch_nodes<Branch>(state, topological_index);
+    auto const node_from_math_id = get_math_id<Node>(state, branch_nodes[0]);
+    auto const node_to_math_id = get_math_id<Node>(state, branch_nodes[1]);
+
+    switch (terminal_type) {
+        using enum MeasuredTerminalType;
+
+    case branch_from:
+        // all power sensors in branch3 are at from side in the mathematical model
+        [[fallthrough]];
+    case branch3_1:
+        [[fallthrough]];
+    case branch3_2:
+        [[fallthrough]];
+    case branch3_3:
+        return current_sensor.template get_output<sym>(solver_output[obj_math_id.group].branch[obj_math_id.pos].i_f,
+                                                       solver_output[node_from_math_id.group].u[node_from_math_id.pos]);
+    case branch_to:
+        return current_sensor.template get_output<sym>(solver_output[obj_math_id.group].branch[obj_math_id.pos].i_t,
+                                                       solver_output[node_to_math_id.group].u[node_to_math_id.pos]);
+    default:
+        throw MissingCaseForEnumError{std::format("{} output_result()", Component::name), terminal_type};
+    }
+}
+template <std::derived_from<GenericCurrentSensor> Component, class ComponentContainer,
           short_circuit_solver_output_type SolverOutputType>
     requires model_component_state_c<MainModelState, ComponentContainer, Component>
 constexpr auto output_result(Component const& power_or_current_sensor,