Basic PGM support for creating current sensors (experimental in state estimation)

code_generation/data/dataset_class_maps/dataset_definitions.json

@@ -52,6 +52,10 @@
                     "names": ["sym_power_sensor", "asym_power_sensor"],
                     "class_name": "PowerSensorInput"
                 },
+                {
+                    "names": ["sym_current_sensor", "asym_current_sensor"],
+                    "class_name": "CurrentSensorInput"
+                },
                 {
                     "names": ["fault"],
                     "class_name": "FaultInput"
@@ -90,6 +94,10 @@
                     "names": ["sym_power_sensor", "asym_power_sensor"],
                     "class_name": "PowerSensorOutput"
                 },
+                {
+                    "names": ["sym_current_sensor", "asym_current_sensor"],
+                    "class_name": "CurrentSensorOutput"
+                },
                 {
                     "names": ["fault"],
                     "class_name": "FaultOutput"
@@ -144,6 +152,10 @@
                     "names": ["sym_power_sensor", "asym_power_sensor"],
                     "class_name": "PowerSensorUpdate"
                 },
+                {
+                    "names": ["sym_current_sensor", "asym_current_sensor"],
+                    "class_name": "CurrentSensorUpdate"
+                },
                 {
                     "names": ["fault"],
                     "class_name": "FaultUpdate"
@@ -178,6 +190,10 @@
                     "names": ["sym_power_sensor", "asym_power_sensor"],
                     "class_name": "SensorShortCircuitOutput"
                 },
+                {
+                    "names": ["sym_current_sensor", "asym_current_sensor"],
+                    "class_name": "SensorShortCircuitOutput"
+                },
                 {
                     "names": ["fault"],
                     "class_name": "FaultShortCircuitOutput"

power_grid_model_c/power_grid_model/include/power_grid_model/all_components.hpp

@@ -9,6 +9,7 @@
 #include "common/component_list.hpp"
 // component include
 #include "component/appliance.hpp"
+#include "component/current_sensor.hpp"
 #include "component/fault.hpp"
 #include "component/generic_branch.hpp"
 #include "component/line.hpp"
@@ -29,6 +30,17 @@ namespace power_grid_model {
 using AllComponents =
     ComponentList<Node, Line, Link, GenericBranch, Transformer, ThreeWindingTransformer, Shunt, Source, SymGenerator,
                   AsymGenerator, SymLoad, AsymLoad, SymPowerSensor, AsymPowerSensor, SymVoltageSensor,
-                  AsymVoltageSensor, Fault, TransformerTapRegulator>;
+                  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/calculation_parameters.hpp

@@ -103,15 +103,19 @@ template <symmetry_tag sym_type> struct CurrentSensorCalcParam {
 };
 
 template <typename T>
-concept sensor_calc_param_type =
-    std::derived_from<T, VoltageSensorCalcParam<symmetric_t>> ||
-    std::derived_from<T, VoltageSensorCalcParam<asymmetric_t>> ||
-    std::derived_from<T, PowerSensorCalcParam<symmetric_t>> || std::derived_from<T, PowerSensorCalcParam<asymmetric_t>>;
+concept sensor_calc_param_type = std::derived_from<T, VoltageSensorCalcParam<symmetric_t>> ||
+                                 std::derived_from<T, VoltageSensorCalcParam<asymmetric_t>> ||
+                                 std::derived_from<T, PowerSensorCalcParam<symmetric_t>> ||
+                                 std::derived_from<T, PowerSensorCalcParam<asymmetric_t>> ||
+                                 std::derived_from<T, CurrentSensorCalcParam<symmetric_t>> ||
+                                 std::derived_from<T, CurrentSensorCalcParam<asymmetric_t>>;
 
 static_assert(sensor_calc_param_type<VoltageSensorCalcParam<symmetric_t>>);
 static_assert(sensor_calc_param_type<VoltageSensorCalcParam<asymmetric_t>>);
 static_assert(sensor_calc_param_type<PowerSensorCalcParam<symmetric_t>>);
 static_assert(sensor_calc_param_type<PowerSensorCalcParam<asymmetric_t>>);
+static_assert(sensor_calc_param_type<CurrentSensorCalcParam<symmetric_t>>);
+static_assert(sensor_calc_param_type<CurrentSensorCalcParam<asymmetric_t>>);
 
 struct TransformerTapRegulatorCalcParam {
     double u_set{};
@@ -150,6 +154,8 @@ struct MathModelTopology {
     DenseGroupedIdxVector power_sensors_per_branch_from;
     DenseGroupedIdxVector power_sensors_per_branch_to;
     DenseGroupedIdxVector power_sensors_per_bus;
+    DenseGroupedIdxVector current_sensors_per_branch_from;
+    DenseGroupedIdxVector current_sensors_per_branch_to;
     DenseGroupedIdxVector tap_regulators_per_branch;
 
     Idx n_bus() const { return static_cast<Idx>(phase_shift.size()); }
@@ -359,6 +365,8 @@ struct ComponentTopology {
     IdxVector voltage_sensor_node_idx;
     IdxVector power_sensor_object_idx; // the index is relative to branch, source, shunt or load_gen
     std::vector<MeasuredTerminalType> power_sensor_terminal_type;
+    IdxVector current_sensor_object_idx; // the index is relative to branch
+    std::vector<MeasuredTerminalType> current_sensor_terminal_type;
     std::vector<ComponentType> regulator_type;
     IdxVector regulated_object_idx; // the index is relative to branch or branch3
     std::vector<ComponentType> regulated_object_type;
@@ -418,7 +426,8 @@ struct TopologicalComponentToMathCoupling {
     std::vector<Idx2D> load_gen;
     std::vector<Idx2D> source;
     std::vector<Idx2D> voltage_sensor;
-    std::vector<Idx2D> power_sensor; // can be coupled to branch-from/to, source, load_gen, or shunt sensor
+    std::vector<Idx2D> power_sensor;   // can be coupled to branch-from/to, source, load_gen, or shunt sensor
+    std::vector<Idx2D> current_sensor; // can be coupled to branch-from/to
     std::vector<Idx2D> regulator;
 };
 

power_grid_model_c/power_grid_model/include/power_grid_model/main_core/input.hpp

@@ -74,7 +74,9 @@ inline void add_component(MainModelState<ComponentContainer>& state, ForwardIter
                 get_component<Branch>(state, measured_object);
                 break;
             case branch3_1:
+                [[fallthrough]];
             case branch3_2:
+                [[fallthrough]];
             case branch3_3:
                 get_component<Branch3>(state, measured_object);
                 break;
@@ -94,11 +96,42 @@ inline void add_component(MainModelState<ComponentContainer>& state, ForwardIter
                 get_component<Node>(state, measured_object);
                 break;
             default:
-                throw MissingCaseForEnumError(std::string(GenericPowerSensor::name) + " item retrieval",
-                                              input.measured_terminal_type);
+                throw MissingCaseForEnumError{std::format("{} item retrieval", GenericPowerSensor::name),
+                                              input.measured_terminal_type};
             }
 
             emplace_component<Component>(state, id, input);
+        } else if constexpr (std::derived_from<Component, GenericCurrentSensor>) {
+            // it is not allowed to place a sensor at a link
+            if (get_component_idx_by_id(state, input.measured_object).group == get_component_type_index<Link>(state)) {
+                throw InvalidMeasuredObject("Link", "CurrentSensor");
+            }
+            // check correctness and get node based on measured terminal type
+            ID const node = [&state, measured_object = input.measured_object,
+                             measured_terminal_type = input.measured_terminal_type] {
+                switch (measured_terminal_type) {
+                    using enum MeasuredTerminalType;
+                    using enum Branch3Side;
+
+                case branch_from:
+                    return get_component<Branch>(state, measured_object).node(BranchSide::from);
+                case branch_to:
+                    return get_component<Branch>(state, measured_object).node(BranchSide::to);
+                case branch3_1:
+                    return get_component<Branch3>(state, measured_object).node(side_1);
+                case branch3_2:
+                    return get_component<Branch3>(state, measured_object).node(side_2);
+                case branch3_3:
+                    return get_component<Branch3>(state, measured_object).node(side_2);
+                default:
+                    throw MissingCaseForEnumError{std::format("{} item retrieval", GenericCurrentSensor::name),
+                                                  measured_terminal_type};
+                }
+            }();
+
+            double const u_rated = get_component<Node>(state, node).u_rated();
+
+            emplace_component<Component>(state, id, input, u_rated);
         } else if constexpr (std::derived_from<Component, Fault>) {
             // check that fault object exists (currently, only faults at nodes are supported)
             get_component<Node>(state, input.fault_object);
@@ -118,18 +151,20 @@ inline void add_component(MainModelState<ComponentContainer>& state, ForwardIter
                     case from:
                         return regulated_object.node(static_cast<BranchSide>(input.control_side));
                     default:
-                        throw MissingCaseForEnumError{std::string{Component::name} + " item retrieval",
+                        throw MissingCaseForEnumError{std::format("{} item retrieval", Component::name),
                                                       input.control_side};
                     }
                 } else if (regulated_object_idx.group == get_component_type_index<ThreeWindingTransformer>(state)) {
                     auto const& regulated_object = get_component<ThreeWindingTransformer>(state, regulated_object_idx);
                     switch (input.control_side) {
                     case side_1:
+                        [[fallthrough]];
                     case side_2:
+                        [[fallthrough]];
                     case side_3:
                         return regulated_object.node(static_cast<Branch3Side>(input.control_side));
                     default:
-                        throw MissingCaseForEnumError{std::string{Component::name} + " item retrieval",
+                        throw MissingCaseForEnumError{std::format("{} item retrieval", Component::name),
                                                       input.control_side};
                     }
                 } else {

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

@@ -66,6 +66,13 @@ constexpr auto comp_base_sequence_cbegin(MainModelState<ComponentContainer> cons
            get_component_sequence_offset<GenericPowerSensor, Component>(state);
 }
 
+template <std::derived_from<GenericCurrentSensor> Component, class ComponentContainer>
+    requires model_component_state_c<MainModelState, ComponentContainer, Component>
+constexpr auto comp_base_sequence_cbegin(MainModelState<ComponentContainer> const& state) {
+    return state.comp_topo->current_sensor_object_idx.cbegin() +
+           get_component_sequence_offset<GenericCurrentSensor, Component>(state);
+}
+
 template <std::same_as<Fault> Component, class ComponentContainer>
     requires model_component_state_c<MainModelState, ComponentContainer, Component>
 constexpr auto comp_base_sequence_cbegin(MainModelState<ComponentContainer> const& state) {
@@ -235,14 +242,14 @@ inline auto output_result(Component const& voltage_sensor, MainModelState<Compon
 }
 
 // output power sensor
-template <std::derived_from<GenericPowerSensor> Component, class ComponentContainer,
+template <power_or_current_sensor_c Component, class ComponentContainer,
           steady_state_solver_output_type SolverOutputType>
     requires model_component_state_c<MainModelState, ComponentContainer, Component>
-constexpr auto output_result(Component const& power_sensor, MainModelState<ComponentContainer> const& state,
+constexpr auto output_result(Component const& power_or_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 = power_sensor.get_terminal_type();
+    auto const terminal_type = power_or_current_sensor.get_terminal_type();
     Idx2D const obj_math_id = [&]() {
         switch (terminal_type) {
             using enum MeasuredTerminalType;
@@ -269,45 +276,55 @@ constexpr auto output_result(Component const& power_sensor, MainModelState<Compo
         case node:
             return state.topo_comp_coup->node[obj_seq];
         default:
-            throw MissingCaseForEnumError(std::string(GenericPowerSensor::name) + " output_result()", terminal_type);
+            throw MissingCaseForEnumError{std::format("{} output_result()", Component::name), terminal_type};
         }
     }();
 
     if (obj_math_id.group == -1) {
-        return power_sensor.template get_null_output<sym>();
+        return power_or_current_sensor.template get_null_output<sym>();
     }
 
     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 power_sensor.template get_output<sym>(solver_output[obj_math_id.group].branch[obj_math_id.pos].s_f);
+        return power_or_current_sensor.template get_output<sym>(
+            solver_output[obj_math_id.group].branch[obj_math_id.pos].s_f);
     case branch_to:
-        return power_sensor.template get_output<sym>(solver_output[obj_math_id.group].branch[obj_math_id.pos].s_t);
+        return power_or_current_sensor.template get_output<sym>(
+            solver_output[obj_math_id.group].branch[obj_math_id.pos].s_t);
     case source:
-        return power_sensor.template get_output<sym>(solver_output[obj_math_id.group].source[obj_math_id.pos].s);
+        return power_or_current_sensor.template get_output<sym>(
+            solver_output[obj_math_id.group].source[obj_math_id.pos].s);
     case shunt:
-        return power_sensor.template get_output<sym>(solver_output[obj_math_id.group].shunt[obj_math_id.pos].s);
+        return power_or_current_sensor.template get_output<sym>(
+            solver_output[obj_math_id.group].shunt[obj_math_id.pos].s);
     case load:
         [[fallthrough]];
     case generator:
-        return power_sensor.template get_output<sym>(solver_output[obj_math_id.group].load_gen[obj_math_id.pos].s);
+        return power_or_current_sensor.template get_output<sym>(
+            solver_output[obj_math_id.group].load_gen[obj_math_id.pos].s);
     case node:
-        return power_sensor.template get_output<sym>(solver_output[obj_math_id.group].bus_injection[obj_math_id.pos]);
+        return power_or_current_sensor.template get_output<sym>(
+            solver_output[obj_math_id.group].bus_injection[obj_math_id.pos]);
     default:
-        throw MissingCaseForEnumError(std::string(GenericPowerSensor::name) + " output_result()", terminal_type);
+        throw MissingCaseForEnumError{std::format("{} output_result()", Component::name), terminal_type};
     }
 }
-template <std::derived_from<GenericPowerSensor> Component, class ComponentContainer,
+template <power_or_current_sensor_c Component, class ComponentContainer,
           short_circuit_solver_output_type SolverOutputType>
     requires model_component_state_c<MainModelState, ComponentContainer, Component>
-constexpr auto output_result(Component const& power_sensor, MainModelState<ComponentContainer> const& /* state */,
+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_sensor.get_null_sc_output();
+    return power_or_current_sensor.get_null_sc_output();
 }
 
 // output fault
@@ -450,5 +467,4 @@ constexpr ResIt output_result(MainModelState<ComponentContainer> const& state,
     res_it = output_result<Shunt>(state, math_output, res_it);
     return res_it;
 }
-
 } // namespace power_grid_model::main_core

power_grid_model_c/power_grid_model/include/power_grid_model/main_core/topology.hpp

@@ -120,7 +120,9 @@ constexpr void register_topology_components(MainModelState<ComponentContainer> c
             case generator:
                 return get_component_sequence_idx<GenericLoadGen>(state, measured_object);
             case branch3_1:
+                [[fallthrough]];
             case branch3_2:
+                [[fallthrough]];
             case branch3_3:
                 return get_component_sequence_idx<Branch3>(state, measured_object);
             case node:
@@ -136,6 +138,38 @@ constexpr void register_topology_components(MainModelState<ComponentContainer> c
         [](GenericPowerSensor const& power_sensor) { return power_sensor.get_terminal_type(); });
 }
 
+template <std::same_as<GenericCurrentSensor> Component, class ComponentContainer>
+    requires model_component_state_c<MainModelState, ComponentContainer, Component>
+constexpr void register_topology_components(MainModelState<ComponentContainer> const& state,
+                                            ComponentTopology& comp_topo) {
+    detail::register_topo_components<Component>(
+        state, comp_topo.current_sensor_object_idx, [&state](GenericCurrentSensor const& current_sensor) {
+            using enum MeasuredTerminalType;
+
+            auto const measured_object = current_sensor.measured_object();
+
+            switch (current_sensor.get_terminal_type()) {
+            case branch_from:
+                [[fallthrough]];
+            case branch_to:
+                return get_component_sequence_idx<Branch>(state, measured_object);
+            case branch3_1:
+                [[fallthrough]];
+            case branch3_2:
+                [[fallthrough]];
+            case branch3_3:
+                return get_component_sequence_idx<Branch3>(state, measured_object);
+            default:
+                throw MissingCaseForEnumError("Current sensor idx to seq transformation",
+                                              current_sensor.get_terminal_type());
+            }
+        });
+
+    detail::register_topo_components<Component>(
+        state, comp_topo.current_sensor_terminal_type,
+        [](GenericCurrentSensor const& current_sensor) { return current_sensor.get_terminal_type(); });
+}
+
 template <std::derived_from<Regulator> Component, class ComponentContainer>
     requires model_component_state_c<MainModelState, ComponentContainer, Component>
 constexpr void register_topology_components(MainModelState<ComponentContainer> const& state,

power_grid_model_c/power_grid_model/include/power_grid_model/main_model.hpp

@@ -14,10 +14,10 @@ namespace power_grid_model {
 
 class MainModel {
   private:
-    using Impl =
-        MainModelImpl<ExtraRetrievableTypes<Base, Node, Branch, Branch3, Appliance, GenericLoadGen, GenericLoad,
-                                            GenericGenerator, GenericPowerSensor, GenericVoltageSensor, Regulator>,
-                      AllComponents>;
+    using Impl = MainModelImpl<
+        ExtraRetrievableTypes<Base, Node, Branch, Branch3, Appliance, GenericLoadGen, GenericLoad, GenericGenerator,
+                              GenericPowerSensor, GenericVoltageSensor, GenericCurrentSensor, Regulator>,
+        AllComponents>;
 
   public:
     using Options = MainModelOptions;
@@ -68,6 +68,10 @@ class MainModel {
 
     CalculationInfo calculation_info() const { return impl().calculation_info(); }
 
+    void check_no_experimental_features_used(Options const& options) const {
+        impl().check_no_experimental_features_used(options);
+    }
+
   private:
     Impl& impl() {
         assert(impl_ != nullptr);