initial implementation for local current sensor

Signed-off-by: Santiago Figueroa Manrique <figueroa1395@gmail.com>

Author: figueroa1395

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,25 +140,32 @@ 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();
         output.energized = 1; // current sensor is always energized
+        auto i_output = i;
         auto const i_calc_param = calc_param<sym_calc>();
-        auto const i_measured_complex = i_calc_param.measurement.value();
-        if (i_calc_param.angle_measurement_type == AngleMeasurementType::global_angle) {
-            output.i_residual = (cabs(i_measured_complex) - cabs(i)) * base_current_;
-            output.i_angle_residual = arg(i_measured_complex) - arg(i);
-        } else {
-            throw NotImplementedError();
+        auto const& angle_measurement_type = i_calc_param.angle_measurement_type;
+        auto const& i_measured_complex = i_calc_param.measurement.value();
+        if (angle_measurement_type == 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);
+            i_output = conj(i_output) * (cos(u_angle) + 1i * sin(u_angle));
         }
+        output.i_residual = (cabs(i_measured_complex) - cabs(i_output)) * base_current_;
+        output.i_angle_residual = arg(i_measured_complex) - arg(i_output);
         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,25 +281,9 @@ 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>();
     }
 
-    auto const get_current_or_power_branch_f = [&solver_output, &obj_math_id]<power_or_current_sensor_c Sensor>() {
-        if constexpr (std::derived_from<Sensor, GenericPowerSensor>) {
-            return solver_output[obj_math_id.group].branch[obj_math_id.pos].s_f;
-        } else if constexpr (std::derived_from<Sensor, GenericCurrentSensor>) {
-            return solver_output[obj_math_id.group].branch[obj_math_id.pos].i_f;
-        }
-    };
-
-    auto const get_current_or_power_branch_t = [&solver_output, &obj_math_id]<power_or_current_sensor_c Sensor>() {
-        if constexpr (std::derived_from<Sensor, GenericPowerSensor>) {
-            return solver_output[obj_math_id.group].branch[obj_math_id.pos].s_t;
-        } else if constexpr (std::derived_from<Sensor, GenericCurrentSensor>) {
-            return solver_output[obj_math_id.group].branch[obj_math_id.pos].i_t;
-        }
-    };
-
     switch (terminal_type) {
         using enum MeasuredTerminalType;
 
@@ -311,25 +295,78 @@ 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>(
-            get_current_or_power_branch_f.template operator()<Component>());
+        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>(
-            get_current_or_power_branch_t.template operator()<Component>());
+        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};
+    }
+}
+
+// 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};
     }
@@ -421,8 +458,7 @@ constexpr ResIt output_result(MainModelState<ComponentContainer> const& state,
                               MathOutput<std::vector<SolverOutputType>> const& math_output, ResIt res_it) {
     return detail::produce_output<Component, Idx2D>(
         state, res_it, [&state, &math_output](Component const& component, Idx2D const math_id) {
-            return output_result<Component>(component, state, math_output.solver_output,
-                                            math_id); // probably here is where all the math_id things are called
+            return output_result<Component>(component, state, math_output.solver_output, math_id);
         });
 }
 template <std::derived_from<Base> Component, class ComponentContainer, solver_output_type SolverOutputType,
@@ -431,8 +467,7 @@ template <std::derived_from<Base> Component, class ComponentContainer, solver_ou
              requires(Component const& component, MainModelState<ComponentContainer> const& state,
                       std::vector<SolverOutputType> const& solver_output, Idx obj_seq) {
                  {
-                     output_result<Component>(component, state, solver_output,
-                                              obj_seq) // probably here is where the sensors things are called
+                     output_result<Component>(component, state, solver_output, obj_seq)
                  } -> detail::assignable_to<std::add_lvalue_reference_t<std::iter_value_t<ResIt>>>;
              }
 constexpr ResIt output_result(MainModelState<ComponentContainer> const& state,

tests/cpp_unit_tests/test_current_sensor.cpp

@@ -70,9 +70,9 @@ TEST_CASE("Test current sensor") {
             CurrentSensorCalcParam<asymmetric_t> asym_sensor_param = sym_current_sensor.calc_param<asymmetric_t>();
 
             CurrentSensorOutput<symmetric_t> const sym_sensor_output =
-                sym_current_sensor.get_output<symmetric_t>(i_sym);
+                sym_current_sensor.get_output<symmetric_t>(i_sym, ComplexValue<symmetric_t>{0.0});
             CurrentSensorOutput<asymmetric_t> sym_sensor_output_asym_param =
-                sym_current_sensor.get_output<asymmetric_t>(i_asym);
+                sym_current_sensor.get_output<asymmetric_t>(i_asym, ComplexValue<asymmetric_t>{0.0});
 
             // Check symmetric sensor output for symmetric parameters
             if constexpr (angle_measurement_type == AngleMeasurementType::global_angle) {