current sensor: math solver implementation

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp

@@ -66,6 +66,8 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
     using sym = sym_type;
 
     static constexpr auto is_iterative = true;
+    static constexpr auto has_current_sensor_implemented =
+        true; // TODO(mgovers): for testing purposes; remove after NRSE has current sensor implemented
 
   private:
     // block size 2 for symmetric, 6 for asym
@@ -146,6 +148,10 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
                                                   &MeasuredValues<sym>::has_branch_to_power};
     static constexpr std::array branch_power_{&MeasuredValues<sym>::branch_from_power,
                                               &MeasuredValues<sym>::branch_to_power};
+    static constexpr std::array has_branch_current_{&MeasuredValues<sym>::has_branch_from_current,
+                                                    &MeasuredValues<sym>::has_branch_to_current};
+    static constexpr std::array branch_current_{&MeasuredValues<sym>::branch_from_current,
+                                                &MeasuredValues<sym>::branch_to_current};
 
     Idx n_bus_;
     // shared topo data
@@ -205,19 +211,32 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
                     }
                     // branch
                     else {
-                        // branch from- and to-side index at 0, and 1 position
-                        IntS const b0 = static_cast<IntS>(type) / 2;
-                        IntS const b1 = static_cast<IntS>(type) % 2;
+                        auto const add_branch_measurement = [&block, &param, obj,
+                                                             type](IntS measured_side,
+                                                                   RealValue<sym> const& branch_current_variance) {
+                            // branch from- and to-side index at 0, and 1 position
+                            IntS const b0 = static_cast<IntS>(type) / 2;
+                            IntS const b1 = static_cast<IntS>(type) % 2;
+                            block.g() += dot(hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b0]),
+                                             diagonal_inverse(branch_current_variance),
+                                             param.branch_param[obj].value[measured_side * 2 + b1]);
+                        };
                         // measured at from-side: 0, to-side: 1
                         for (IntS const measured_side : std::array<IntS, 2>{0, 1}) {
                             // has measurement
                             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() +=
-                                    dot(hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b0]),
-                                        diagonal_inverse(static_cast<IndependentComplexRandVar<sym>>(power).variance),
-                                        param.branch_param[obj].value[measured_side * 2 + b1]);
+                                // assume voltage ~ 1 p.u. => current variance = power variance * 1² = power variance
+                                add_branch_measurement(measured_side,
+                                                       static_cast<IndependentComplexRandVar<sym>>(power).variance);
+                            }
+                            if (std::invoke(has_branch_current_[measured_side], measured_value, obj)) {
+                                // G += (variance^-1)
+                                auto const& current =
+                                    std::invoke(branch_current_[measured_side], measured_value, obj).measurement;
+                                add_branch_measurement(measured_side,
+                                                       static_cast<IndependentComplexRandVar<sym>>(current).variance);
                             }
                         }
                     }
@@ -296,23 +315,42 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
                 }
                 // branch
                 else {
-                    // branch is either ff or tt
-                    IntS const b = static_cast<IntS>(type) / 2;
-                    assert(b == static_cast<IntS>(type) % 2);
+                    auto const add_branch_measurement = [&rhs_block, &param, obj,
+                                                         type](IntS measured_side,
+                                                               IndependentComplexRandVar<sym> const& current) {
+                        // branch is either ff or tt
+                        IntS const b = static_cast<IntS>(type) / 2;
+                        // eta += Y{side, b}^H * (variance^-1) * i_branch_{f, t}
+                        rhs_block.eta() +=
+                            dot(hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b]),
+
+                                diagonal_inverse(current.variance), current.value);
+                    };
                     // measured at from-side: 0, to-side: 1
                     for (IntS const measured_side : std::array<IntS, 2>{0, 1}) {
+                        // the current needs to be calculated with the voltage of the measured bus side
+                        // NOTE: not the bus that is currently being processed!
+                        Idx const measured_bus = branch_bus_idx[obj][measured_side];
+
                         // has measurement
                         if (std::invoke(has_branch_power_[measured_side], measured_value, obj)) {
                             PowerSensorCalcParam<sym> 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
-                            // NOTE: not the current bus!
-                            Idx const measured_bus = branch_bus_idx[obj][measured_side];
-                            // eta += Y{side, b}^H * (variance^-1) * i_branch_{f, t}
-                            rhs_block.eta() +=
-                                dot(hermitian_transpose(param.branch_param[obj].value[measured_side * 2 + b]),
-                                    diagonal_inverse(static_cast<IndependentComplexRandVar<sym>>(m).variance),
-                                    conj(m.value() / u[measured_bus]));
+                            auto const apparent_power = static_cast<IndependentComplexRandVar<sym>>(m);
+                            add_branch_measurement(measured_side,
+                                                   {.value = conj(apparent_power.value / u[measured_bus]),
+                                                    .variance = apparent_power.variance});
+                        }
+                        if (std::invoke(has_branch_current_[measured_side], measured_value, obj)) {
+                            CurrentSensorCalcParam<sym> const m =
+                                std::invoke(branch_current_[measured_side], measured_value, obj);
+                            // for local angle current sensors, the current needs to be offset with the phase offset
+                            // of the measured bus side. NOTE: not the bus that is currently being processed!
+                            auto measured_current = static_cast<IndependentComplexRandVar<sym>>(m.measurement);
+                            if (m.angle_measurement_type == AngleMeasurementType::local_angle) {
+                                measured_current.value *= phase_shift(u[measured_bus]); // offset with the phase shift
+                            }
+                            add_branch_measurement(measured_side, measured_current);
                         }
                     }
                 }

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp

@@ -91,6 +91,8 @@ template <symmetry_tag sym_type> class NewtonRaphsonSESolver {
     using sym = sym_type;
 
     static constexpr auto is_iterative = true;
+    static constexpr auto has_current_sensor_implemented =
+        false; // TODO(mgovers): for testing purposes; remove after NRSE has current sensor implemented
 
   private:
     enum class Order : IntS { row_major = 0, column_major = 1 };

tests/cpp_unit_tests/test_math_solver_se.hpp

@@ -344,6 +344,80 @@ TEST_CASE_TEMPLATE_DEFINE("Test math solver - SE, measurements", SolverType, tes
         CHECK(real(output.branch[0].s_f) == doctest::Approx(1.95));
     }
 
+    SUBCASE("Source and branch current") {
+        /*
+        network, v means voltage measured, c means current measured
+
+         bus_0(v) -(c)-branch_0-- bus_1
+            |                       |
+        source_0(p)               load_0
+
+        */
+        // to make the distinction between global and local angle current measurement
+        auto const global_shift = phase_shift(1.0 + 1.0i);
+
+        topo.power_sensors_per_source = {from_sparse, {0, 1}};
+        topo.current_sensors_per_branch_from = {from_sparse, {0, 1}};
+        se_input.measured_voltage = {{.value = 1.0 * global_shift, .variance = 0.1}};
+
+        auto param_ptr = std::make_shared<MathModelParam<symmetric_t> const>(param);
+        auto topo_ptr = std::make_shared<MathModelTopology const>(topo);
+        YBus<symmetric_t> const y_bus_sym{topo_ptr, param_ptr};
+
+        SolverType solver{y_bus_sym, topo_ptr};
+
+        SUBCASE("Local angle current sensor") {
+            se_input.measured_source_power = {{.real_component = {.value = 1.93, .variance = 0.05},
+                                               .imag_component = {.value = 0.0, .variance = 0.05}}};
+            se_input.measured_branch_from_current = {
+                {.angle_measurement_type = AngleMeasurementType::local_angle,
+                 .measurement = {.real_component = {.value = 1.97, .variance = 0.05},
+                                 .imag_component = {.value = 0.0, .variance = 0.05}}}};
+
+            output = run_state_estimation(solver, y_bus_sym, se_input, error_tolerance, num_iter, info);
+
+            if (SolverType::has_current_sensor_implemented) { // TODO(mgovers): for testing purposes; remove if
+                                                              // statement after NRSE has current sensor implemented
+                CHECK(real(output.bus_injection[0]) == doctest::Approx(1.95));
+                CHECK(real(output.source[0].s) == doctest::Approx(1.95));
+                CHECK(real(output.branch[0].s_f) == doctest::Approx(1.95));
+                CHECK(real(output.branch[0].i_f) == doctest::Approx(real(1.95 * global_shift)));
+                CHECK(imag(output.branch[0].i_f) == doctest::Approx(imag(1.95 * global_shift)));
+            } else {
+                CHECK_FALSE(real(output.bus_injection[0]) == doctest::Approx(1.95));
+                CHECK_FALSE(real(output.source[0].s) == doctest::Approx(1.95));
+                CHECK_FALSE(real(output.branch[0].s_f) == doctest::Approx(1.95));
+                CHECK_FALSE(real(output.branch[0].i_f) == doctest::Approx(real(1.95 * global_shift)));
+                CHECK_FALSE(imag(output.branch[0].i_f) == doctest::Approx(imag(1.95 * global_shift)));
+            }
+        }
+        SUBCASE("Global angle current sensor") {
+            se_input.measured_source_power = {{.real_component = {.value = 1.93, .variance = 0.05},
+                                               .imag_component = {.value = 0.0, .variance = 0.05}}};
+            se_input.measured_branch_from_current = {
+                {.angle_measurement_type = AngleMeasurementType::global_angle,
+                 .measurement = {.real_component = {.value = real(1.97 * global_shift), .variance = 0.05},
+                                 .imag_component = {.value = imag(1.97 * global_shift), .variance = 0.05}}}};
+
+            output = run_state_estimation(solver, y_bus_sym, se_input, error_tolerance, num_iter, info);
+
+            if (SolverType::has_current_sensor_implemented) { // TODO(mgovers): for testing purposes; remove if
+                                                              // statement after NRSE has current sensor implemented
+                CHECK(real(output.bus_injection[0]) == doctest::Approx(1.95));
+                CHECK(real(output.source[0].s) == doctest::Approx(1.95));
+                CHECK(real(output.branch[0].s_f) == doctest::Approx(1.95));
+                CHECK(real(output.branch[0].i_f) == doctest::Approx(real(1.95 * global_shift)));
+                CHECK(imag(output.branch[0].i_f) == doctest::Approx(imag(1.95 * global_shift)));
+            } else {
+                CHECK_FALSE(real(output.bus_injection[0]) == doctest::Approx(1.95));
+                CHECK_FALSE(real(output.source[0].s) == doctest::Approx(1.95));
+                CHECK_FALSE(real(output.branch[0].s_f) == doctest::Approx(1.95));
+                CHECK_FALSE(real(output.branch[0].i_f) == doctest::Approx(real(1.95 * global_shift)));
+                CHECK_FALSE(imag(output.branch[0].i_f) == doctest::Approx(imag(1.95 * global_shift)));
+            }
+        }
+    }
+
     SUBCASE("Load and branch") {
         /*
         network, v means voltage measured, p means power measured

tests/cpp_validation_tests/test_validation.cpp

@@ -295,6 +295,8 @@ std::map<std::string, PGM_TapChangingStrategy, std::less<>> const optimizer_stra
     {"min_voltage_tap", PGM_tap_changing_strategy_min_voltage_tap},
     {"max_voltage_tap", PGM_tap_changing_strategy_max_voltage_tap},
     {"fast_any_tap", PGM_tap_changing_strategy_fast_any_tap}};
+std::map<std::string, PGM_ExperimentalFeatures, std::less<>> const experimental_features_mapping = {
+    {"disabled", PGM_experimental_features_disabled}, {"enabled", PGM_experimental_features_enabled}};
 
 // case parameters
 struct CaseParam {
@@ -304,6 +306,7 @@ struct CaseParam {
     std::string calculation_method;
     std::string short_circuit_voltage_scaling;
     std::string tap_changing_strategy;
+    std::string experimental_features;
     double err_tol = 1e-8;
     Idx max_iter = 20;
     bool sym{};
@@ -329,6 +332,7 @@ Options get_options(CaseParam const& param, Idx threading = -1) {
     options.set_threading(threading);
     options.set_short_circuit_voltage_scaling(sc_voltage_scaling_mapping.at(param.short_circuit_voltage_scaling));
     options.set_tap_changing_strategy(optimizer_strategy_mapping.at(param.tap_changing_strategy));
+    options.set_experimental_features(experimental_features_mapping.at(param.experimental_features));
     return options;
 }
 
@@ -390,6 +394,7 @@ std::optional<CaseParam> construct_case(std::filesystem::path const& case_dir, j
     }
 
     param.tap_changing_strategy = calculation_method_params.value("tap_changing_strategy", "disabled");
+    param.experimental_features = calculation_method_params.value("experimental_features", "disabled");
     param.case_name += sym ? "-sym"s : "-asym"s;
     param.case_name += "-"s + param.calculation_method;
     param.case_name += is_batch ? "_batch"s : ""s;

tests/data/state_estimation/basic-current-sensor/params.json

@@ -10,11 +10,7 @@
   },
   "extra_params": {
     "iterative_linear": {
-      "experimental_features": "enabled",
-      "fail": {
-        "raises": "SparseMatrixError",
-        "reason": "Current sensors are not yet implemented for this calculation method"
-      }
+      "experimental_features": "enabled"
     },
     "newton_raphson": {
       "experimental_features": "enabled",