Fix NRSE line measurement implementation

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/measured_values.hpp

@@ -242,7 +242,7 @@ template <bool sym> class MeasuredValues {
 
         // assign a meaningful mean angle shift, if at least one voltage has angle measurement
         if (has_angle()) {
-            mean_angle_shift_ = angle_cum / n_voltage_angle_measurements_;
+            mean_angle_shift_ = angle_cum / static_cast<RealValue<sym>>(n_voltage_angle_measurements_);
         }
 
         static constexpr auto const is_measured = [](auto const& value) { return value >= 0; };

power_grid_model_c/power_grid_model/include/power_grid_model/three_phase_tensor.hpp

@@ -76,7 +76,7 @@ template <scalar_value T> class Tensor : public Eigen3Tensor<T> {
 
 template <scalar_value T> class DiagonalTensor : public Eigen3DiagonalTensor<T> {
   public:
-    DiagonalTensor() { (*this) = Eigen3DiagonalTensor<T>::Zero(); }
+    DiagonalTensor() { (*this).setZero(); }
     // additional constructors
     explicit DiagonalTensor(T const& x) : Eigen3DiagonalTensor<T>{x, x, x} {}
     explicit DiagonalTensor(Vector<T> const& v) : Eigen3DiagonalTensor<T>{v(0), v(1), v(2)} {}

tests/cpp_unit_tests/test_main_model_se.cpp

@@ -211,6 +211,48 @@ TEST_CASE_TEMPLATE("Test main model - state estimation", CalculationMethod, Iter
                 }
             }
         }
+        SUBCASE("Line power sensor") {
+            main_model.add_component<Node>({{1, 10e3}, {2, 10e3}});
+            main_model.add_component<Line>({{3, 1, 2, 1, 1, 0.01, 0.01, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1e3}});
+            main_model.add_component<Source>({{4, 1, 1, 1.0, nan, nan, nan, nan}});
+            main_model.add_component<Shunt>({{6, 2, 1, 1800 / 10e3 / 10e3, -180 / 10e3 / 10e3, 0.0, 0.0}});
+            main_model.add_component<SymVoltageSensor>({{11, 1, 1e2, 10.0e3, 0.0}});
+            SUBCASE("Symmetric Power Sensor - Symmetric Calculation") {
+                main_model.add_component<SymPowerSensor>(
+                    {{17, 3, MeasuredTerminalType::branch_from, 1e2, 1800.0, 180.0, nan, nan},
+                     {18, 3, MeasuredTerminalType::branch_to, 1e2, -1800.0, -180.0, nan, nan},
+                     {16, 6, MeasuredTerminalType::shunt, 1e2, 1800.0, 180.0, nan, nan}});
+                SUBCASE("Line flow") {
+                    main_model.set_construction_complete();
+                    std::vector<MathOutput<true>> const math_output =
+                        main_model.calculate_state_estimation<true>(1e-8, 20, calculation_method);
+
+                    std::vector<SymApplianceOutput> shunt_output(1);
+                    std::vector<SymNodeOutput> node_output(2);
+                    std::vector<SymPowerSensorOutput> power_sensor_output(3);
+                    std::vector<BranchOutput<true>> line_output(1);
+                    main_model.output_result<Shunt>(math_output, shunt_output.begin());
+                    main_model.output_result<Node>(math_output, node_output.begin());
+                    main_model.output_result<Line>(math_output, line_output.begin());
+                    main_model.output_result<SymPowerSensor>(math_output, power_sensor_output.begin());
+
+                    CHECK(shunt_output[0].p == doctest::Approx(1800.0).scale(1e3));
+                    CHECK(shunt_output[0].q == doctest::Approx(180.0).scale(1e3));
+
+                    CHECK(line_output[0].p_from == doctest::Approx(1800.0).scale(1e3));
+                    CHECK(line_output[0].q_from == doctest::Approx(180.0).scale(1e3));
+                    CHECK(line_output[0].p_to == doctest::Approx(-1800.0).scale(1e3));
+                    CHECK(line_output[0].q_to == doctest::Approx(-180.0).scale(1e3));
+
+                    CHECK(power_sensor_output[0].p_residual == doctest::Approx(0.0).scale(1e3)); // shunt
+                    CHECK(power_sensor_output[0].q_residual == doctest::Approx(0.0).scale(1e3)); // shunt
+                    CHECK(power_sensor_output[1].p_residual == doctest::Approx(0.0).scale(1e3)); // branch_from
+                    CHECK(power_sensor_output[1].q_residual == doctest::Approx(0.0).scale(1e3)); // branch_from
+                    CHECK(power_sensor_output[2].p_residual == doctest::Approx(0.0).scale(1e3)); // branch_to
+                    CHECK(power_sensor_output[2].q_residual == doctest::Approx(0.0).scale(1e3)); // branch_to
+                }
+            }
+        }
         SUBCASE("Forbid Link Power Measurements") {
             main_model.add_component<Node>({{1, 10e3}, {2, 10e3}});
             main_model.add_component<Link>({{3, 1, 2, 1, 1}});

tests/cpp_unit_tests/test_math_solver.cpp

@@ -521,24 +521,45 @@ TEST_CASE("Test math solver") {
     SUBCASE("Test sym se with angle") {
         MathSolver<true> solver{topo_ptr};
         CalculationInfo info;
-        MathOutput<true> const output =
-            solver.run_state_estimation(se_input_angle, 1e-10, 20, info, iterative_linear, y_bus_sym);
+        MathOutput<true> output;
+
+        SUBCASE("iterative linear") {
+            output = solver.run_state_estimation(se_input_angle, 1e-10, 20, info, iterative_linear, y_bus_sym);
+        }
+        SUBCASE("Newton-Raphson") {
+            output = solver.run_state_estimation(se_input_angle, 1e-10, 20, info, newton_raphson, y_bus_sym);
+        }
+
         assert_output(output, output_ref);
     }
 
     SUBCASE("Test sym se without angle") {
         MathSolver<true> solver{topo_ptr};
         CalculationInfo info;
-        MathOutput<true> const output =
-            solver.run_state_estimation(se_input_no_angle, 1e-10, 20, info, iterative_linear, y_bus_sym);
+        MathOutput<true> output;
+
+        SUBCASE("iterative linear") {
+            output = solver.run_state_estimation(se_input_no_angle, 1e-10, 20, info, iterative_linear, y_bus_sym);
+        }
+        SUBCASE("Newton-Raphson") {
+            output = solver.run_state_estimation(se_input_no_angle, 1e-10, 20, info, newton_raphson, y_bus_sym);
+        }
+
         assert_output(output, output_ref, true);
     }
 
     SUBCASE("Test sym se with angle, const z") {
         MathSolver<true> solver{topo_ptr};
         CalculationInfo info;
-        MathOutput<true> const output =
-            solver.run_state_estimation(se_input_angle_const_z, 1e-10, 20, info, iterative_linear, y_bus_sym);
+        MathOutput<true> output;
+
+        SUBCASE("iterative linear") {
+            output = solver.run_state_estimation(se_input_angle_const_z, 1e-10, 20, info, iterative_linear, y_bus_sym);
+        }
+        SUBCASE("Newton-Raphson") {
+            output = solver.run_state_estimation(se_input_angle_const_z, 1e-10, 20, info, newton_raphson, y_bus_sym);
+        }
+
         assert_output(output, output_ref_z);
     }
 
@@ -548,32 +569,62 @@ TEST_CASE("Test math solver") {
         auto& branch_from_power = se_input_angle.measured_branch_from_power.front();
         branch_from_power.p_variance = 0.25;
         branch_from_power.q_variance = 0.75;
-        MathOutput<true> const output =
-            solver.run_state_estimation(se_input_angle, 1e-10, 20, info, iterative_linear, y_bus_sym);
+        MathOutput<true> output;
+
+        SUBCASE("iterative linear") {
+            output = solver.run_state_estimation(se_input_angle, 1e-10, 20, info, iterative_linear, y_bus_sym);
+        }
+        SUBCASE("Newton-Raphson") {
+            output = solver.run_state_estimation(se_input_angle, 1e-10, 20, info, newton_raphson, y_bus_sym);
+        }
+
         assert_output(output, output_ref);
     }
 
     SUBCASE("Test asym se with angle") {
         MathSolver<false> solver{topo_ptr};
         CalculationInfo info;
-        MathOutput<false> const output =
-            solver.run_state_estimation(se_input_asym_angle, 1e-10, 20, info, iterative_linear, y_bus_asym);
+        MathOutput<false> output;
+
+        SUBCASE("iterative linear") {
+            output = solver.run_state_estimation(se_input_asym_angle, 1e-10, 20, info, iterative_linear, y_bus_asym);
+        }
+        SUBCASE("Newton-Raphson") {
+            output = solver.run_state_estimation(se_input_asym_angle, 1e-10, 20, info, newton_raphson, y_bus_asym);
+        }
+
         assert_output(output, output_ref_asym);
     }
 
     SUBCASE("Test asym se without angle") {
         MathSolver<false> solver{topo_ptr};
         CalculationInfo info;
-        MathOutput<false> const output =
-            solver.run_state_estimation(se_input_asym_no_angle, 1e-10, 20, info, iterative_linear, y_bus_asym);
+        MathOutput<false> output;
+
+        SUBCASE("iterative linear") {
+            output = solver.run_state_estimation(se_input_asym_no_angle, 1e-10, 20, info, iterative_linear, y_bus_asym);
+        }
+        SUBCASE("Newton-Raphson") {
+            output = solver.run_state_estimation(se_input_asym_no_angle, 1e-10, 20, info, newton_raphson, y_bus_asym);
+        }
+
         assert_output(output, output_ref_asym, true);
     }
 
     SUBCASE("Test asym se with angle, const z") {
         MathSolver<false> solver{topo_ptr};
         CalculationInfo info;
-        MathOutput<false> const output =
-            solver.run_state_estimation(se_input_asym_angle_const_z, 1e-10, 20, info, iterative_linear, y_bus_asym);
+        MathOutput<false> output;
+
+        SUBCASE("iterative linear") {
+            output =
+                solver.run_state_estimation(se_input_asym_angle_const_z, 1e-10, 20, info, iterative_linear, y_bus_asym);
+        }
+        SUBCASE("Newton-Raphson") {
+            output =
+                solver.run_state_estimation(se_input_asym_angle_const_z, 1e-10, 20, info, newton_raphson, y_bus_asym);
+        }
+
         assert_output(output, output_ref_asym_z);
     }
 
@@ -583,8 +634,15 @@ TEST_CASE("Test math solver") {
         auto& branch_from_power = se_input_asym_angle.measured_branch_from_power.front();
         branch_from_power.p_variance = RealValue<false>{0.25};
         branch_from_power.q_variance = RealValue<false>{0.75};
-        MathOutput<false> const output =
-            solver.run_state_estimation(se_input_asym_angle, 1e-10, 20, info, iterative_linear, y_bus_asym);
+        MathOutput<false> output;
+
+        SUBCASE("iterative linear") {
+            output = solver.run_state_estimation(se_input_asym_angle, 1e-10, 20, info, iterative_linear, y_bus_asym);
+        }
+        SUBCASE("Newton-Raphson") {
+            output = solver.run_state_estimation(se_input_asym_angle, 1e-10, 20, info, newton_raphson, y_bus_asym);
+        }
+
         assert_output(output, output_ref_asym);
     }
 }
@@ -1065,7 +1123,15 @@ TEST_CASE("Math solver, zero variance test") {
 
     MathSolver<true> solver{topo_ptr};
     CalculationInfo info;
-    MathOutput<true> output = solver.run_state_estimation(se_input, 1e-10, 20, info, iterative_linear, y_bus_sym);
+    MathOutput<true> output;
+
+    SUBCASE("iterative linear") {
+        output = solver.run_state_estimation(se_input, 1e-10, 20, info, iterative_linear, y_bus_sym);
+    }
+    SUBCASE("Newton-Raphson") {
+        // TODO(mgovers): to be added in #478
+        output = solver.run_state_estimation(se_input, 1e-10, 20, info, newton_raphson, y_bus_sym);
+    }
 
     // check both voltage
     check_close(output.u[0], 1.0);
@@ -1126,8 +1192,15 @@ TEST_CASE("Math solver, measurements") {
         auto param_ptr = std::make_shared<MathModelParam<true> const>(param);
         auto topo_ptr = std::make_shared<MathModelTopology const>(topo);
         YBus<true> const y_bus_sym{topo_ptr, param_ptr};
+
         MathSolver<true> solver{topo_ptr};
-        output = solver.run_state_estimation(se_input, 1e-10, 20, info, iterative_linear, y_bus_sym);
+
+        SUBCASE("iterative linear") {
+            output = solver.run_state_estimation(se_input, 1e-10, 20, info, iterative_linear, y_bus_sym);
+        }
+        SUBCASE("Newton-Raphson") {
+            output = solver.run_state_estimation(se_input, 1e-10, 20, info, newton_raphson, y_bus_sym);
+        }
 
         CHECK(real(output.bus_injection[0]) == doctest::Approx(1.95));
         CHECK(real(output.source[0].s) == doctest::Approx(1.95));
@@ -1153,7 +1226,13 @@ TEST_CASE("Math solver, measurements") {
         auto topo_ptr = std::make_shared<MathModelTopology const>(topo);
         YBus<true> const y_bus_sym{topo_ptr, param_ptr};
         MathSolver<true> solver{topo_ptr};
-        output = solver.run_state_estimation(se_input, 1e-10, 20, info, iterative_linear, y_bus_sym);
+
+        SUBCASE("iterative linear") {
+            output = solver.run_state_estimation(se_input, 1e-10, 20, info, iterative_linear, y_bus_sym);
+        }
+        SUBCASE("Newton-Raphson") {
+            output = solver.run_state_estimation(se_input, 1e-10, 20, info, newton_raphson, y_bus_sym);
+        }
 
         CHECK(real(output.bus_injection[1]) == doctest::Approx(-1.95));
         CHECK(real(output.load_gen[0].s) == doctest::Approx(-1.95));
@@ -1181,7 +1260,13 @@ TEST_CASE("Math solver, measurements") {
         auto topo_ptr = std::make_shared<MathModelTopology const>(topo);
         YBus<true> const y_bus_sym{topo_ptr, param_ptr};
         MathSolver<true> solver{topo_ptr};
-        output = solver.run_state_estimation(se_input, 1e-10, 20, info, iterative_linear, y_bus_sym);
+
+        SUBCASE("iterative linear") {
+            output = solver.run_state_estimation(se_input, 1e-10, 20, info, iterative_linear, y_bus_sym);
+        }
+        SUBCASE("Newton-Raphson") {
+            output = solver.run_state_estimation(se_input, 1e-10, 20, info, newton_raphson, y_bus_sym);
+        }
 
         CHECK(real(output.bus_injection[0]) == doctest::Approx(2.0));
         CHECK(real(output.source[0].s) == doctest::Approx(2.0));
@@ -1209,7 +1294,13 @@ TEST_CASE("Math solver, measurements") {
         auto topo_ptr = std::make_shared<MathModelTopology const>(topo);
         YBus<true> const y_bus_sym{topo_ptr, param_ptr};
         MathSolver<true> solver{topo_ptr};
-        output = solver.run_state_estimation(se_input, 1e-10, 20, info, iterative_linear, y_bus_sym);
+
+        SUBCASE("iterative linear") {
+            output = solver.run_state_estimation(se_input, 1e-10, 20, info, iterative_linear, y_bus_sym);
+        }
+        SUBCASE("Newton-Raphson") {
+            output = solver.run_state_estimation(se_input, 1e-10, 20, info, newton_raphson, y_bus_sym);
+        }
 
         CHECK(real(output.bus_injection[0]) == doctest::Approx(2.0));
         CHECK(real(output.source[0].s) == doctest::Approx(2.0));
@@ -1237,7 +1328,13 @@ TEST_CASE("Math solver, measurements") {
         auto topo_ptr = std::make_shared<MathModelTopology const>(topo);
         YBus<true> const y_bus_sym{topo_ptr, param_ptr};
         MathSolver<true> solver{topo_ptr};
-        output = solver.run_state_estimation(se_input, 1e-10, 20, info, iterative_linear, y_bus_sym);
+
+        SUBCASE("iterative linear") {
+            output = solver.run_state_estimation(se_input, 1e-10, 20, info, iterative_linear, y_bus_sym);
+        }
+        SUBCASE("Newton-Raphson") {
+            output = solver.run_state_estimation(se_input, 1e-10, 20, info, newton_raphson, y_bus_sym);
+        }
 
         CHECK(real(output.bus_injection[1]) == doctest::Approx(-2.0));
         CHECK(real(output.load_gen[0].s) == doctest::Approx(-2.0));
@@ -1264,7 +1361,13 @@ TEST_CASE("Math solver, measurements") {
         auto topo_ptr = std::make_shared<MathModelTopology const>(topo);
         YBus<true> const y_bus_sym{topo_ptr, param_ptr};
         MathSolver<true> solver{topo_ptr};
-        output = solver.run_state_estimation(se_input, 1e-10, 20, info, iterative_linear, y_bus_sym);
+
+        SUBCASE("iterative linear") {
+            output = solver.run_state_estimation(se_input, 1e-10, 20, info, iterative_linear, y_bus_sym);
+        }
+        SUBCASE("Newton-Raphson") {
+            output = solver.run_state_estimation(se_input, 1e-10, 20, info, newton_raphson, y_bus_sym);
+        }
 
         CHECK(real(output.bus_injection[1]) == doctest::Approx(-2.0));
         CHECK(real(output.branch[0].s_t) == doctest::Approx(-2.0));
@@ -1294,7 +1397,13 @@ TEST_CASE("Math solver, measurements") {
         auto topo_ptr = std::make_shared<MathModelTopology const>(topo);
         YBus<true> const y_bus_sym{topo_ptr, param_ptr};
         MathSolver<true> solver{topo_ptr};
-        output = solver.run_state_estimation(se_input, 1e-10, 20, info, iterative_linear, y_bus_sym);
+
+        SUBCASE("iterative linear") {
+            output = solver.run_state_estimation(se_input, 1e-10, 20, info, iterative_linear, y_bus_sym);
+        }
+        SUBCASE("Newton-Raphson") {
+            output = solver.run_state_estimation(se_input, 1e-10, 20, info, newton_raphson, y_bus_sym);
+        }
 
         CHECK(real(output.bus_injection[1]) == doctest::Approx(-1.0));
         CHECK(real(output.load_gen[0].s) == doctest::Approx(-1.85));
@@ -1323,7 +1432,13 @@ TEST_CASE("Math solver, measurements") {
         auto topo_ptr = std::make_shared<MathModelTopology const>(topo);
         YBus<true> const y_bus_sym{topo_ptr, param_ptr};
         MathSolver<true> solver{topo_ptr};
-        output = solver.run_state_estimation(se_input, 1e-10, 20, info, iterative_linear, y_bus_sym);
+
+        SUBCASE("iterative linear") {
+            output = solver.run_state_estimation(se_input, 1e-10, 20, info, iterative_linear, y_bus_sym);
+        }
+        SUBCASE("Newton-Raphson") {
+            output = solver.run_state_estimation(se_input, 1e-10, 20, info, newton_raphson, y_bus_sym);
+        }
 
         // the different aggregation of the load gen's P and Q measurements cause differences compared to the case with
         // identical variances

tests/cpp_unit_tests/test_serializer.cpp

@@ -156,7 +156,6 @@ constexpr std::string_view batch_dataset_dict_indent =
 } // namespace
 
 TEST_CASE("Serializer") {
-
     std::vector<SymLoadGenUpdate> sym_load_gen(4);
     meta_data.get_dataset("update").get_component("sym_load").set_nan(sym_load_gen.data(), 0, 4);
     sym_load_gen[0].id = 9;