Fix NRSE line measurement implementation

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

@@ -189,7 +189,7 @@ template <bool sym> class IterativeLinearSESolver {
                     // shunt
                     if (type == YBusElementType::shunt) {
                         if (measured_value.has_shunt(obj)) {
-                            // G += Ys^H * (variance^-1) * Ys
+                            // G += (-Ys)^H * (variance^-1) * (-Ys)
                             auto const& shunt_power = measured_value.shunt_power(obj);
                             block.g() += dot(hermitian_transpose(param.shunt_param[obj]),
                                              diagonal_inverse(shunt_power.p_variance + shunt_power.q_variance),
@@ -282,7 +282,7 @@ template <bool sym> class IterativeLinearSESolver {
                 if (type == YBusElementType::shunt) {
                     if (measured_value.has_shunt(obj)) {
                         PowerSensorCalcParam<sym> const& m = measured_value.shunt_power(obj);
-                        // eta -= Ys^H * (variance^-1) * i_shunt
+                        // eta += (-Ys)^H * (variance^-1) * i_shunt
                         rhs_block.eta() -= dot(hermitian_transpose(param.shunt_param[obj]),
                                                diagonal_inverse(m.p_variance + m.q_variance), conj(m.value / u[bus]));
                     }

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 / static_cast<RealValue<sym>>(n_voltage_angle_measurements_);
+            mean_angle_shift_ = angle_cum / static_cast<RealValue<sym> const>(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/math_solver/newton_raphson_se_solver.hpp

@@ -111,13 +111,6 @@ template <bool sym> class NewtonRaphsonSESolver {
             this->dQ_dv += other.dQ_dv;
             return *this;
         }
-        NRSEJacobian& operator-=(NRSEJacobian const& other) {
-            this->dP_dt -= other.dP_dt;
-            this->dP_dv -= other.dP_dv;
-            this->dQ_dt -= other.dQ_dt;
-            this->dQ_dv -= other.dQ_dv;
-            return *this;
-        }
     };
 
   public:
@@ -221,7 +214,7 @@ template <bool sym> class NewtonRaphsonSESolver {
         std::ranges::fill(x_, default_unknown);
     }
 
-    void prepare_matrix_and_rhs(YBus<sym> const& y_bus, MeasuredValues<sym> const& measured_value,
+    void prepare_matrix_and_rhs(YBus<sym> const& y_bus, MeasuredValues<sym> const& measured_values,
                                 ComplexValueVector<sym> const& current_u) {
         MathModelParam<sym> const& param = y_bus.math_model_param();
         IdxVector const& row_indptr = y_bus.row_indptr_lu();
@@ -261,78 +254,89 @@ template <bool sym> class NewtonRaphsonSESolver {
 
                 // get a reference and reset block to zero
                 NRSEGainBlock<sym>& block = data_gain_[data_idx_lu];
-                // Diagonal block is being cleared outside this loop
-                if (row != col) {
-                    block.clear();
-                }
-                // fill block with voltage measurement, only diagonal
                 if (row == col) {
-                    process_voltage_measurements(block, rhs_block, measured_value, row);
+                    // fill block with voltage measurement, only diagonal
+                    process_voltage_measurements(diag_block, rhs_block, measured_values, row);
+                } else {
+                    // Diagonal block is being cleared outside this loop
+                    block.clear();
                 }
+
                 // fill block with branch, shunt measurement
                 for (Idx element_idx = y_bus.y_bus_entry_indptr()[data_idx];
                      element_idx != y_bus.y_bus_entry_indptr()[data_idx + 1]; ++element_idx) {
                     Idx const obj = y_bus.y_bus_element()[element_idx].idx;
                     YBusElementType const type = y_bus.y_bus_element()[element_idx].element_type;
-                    if (type == YBusElementType::shunt) {
-                        if (measured_value.has_shunt(obj)) {
-                            auto const& yii = param.shunt_param[obj];
-                            auto const& measured_power = measured_value.shunt_power(obj);
-                            process_shunt_measurement(block, rhs_block, yii, u_state, measured_power);
+
+                    switch (type) {
+                    case YBusElementType::shunt:
+                        if (measured_values.has_shunt(obj)) {
+                            auto const& ys = param.shunt_param[obj];
+                            auto const& measured_power = measured_values.shunt_power(obj);
+                            process_shunt_measurement(block, rhs_block, ys, u_state, measured_power);
                         }
-                    } else if (type == YBusElementType::bft || type == YBusElementType::btf) {
+                        break;
+                    case YBusElementType::bft:
+                        [[fallthrough]];
+                    case YBusElementType::btf: {
                         auto const& y_branch = param.branch_param[obj];
-                        if (measured_value.has_branch_from(obj)) {
+                        if (measured_values.has_branch_from(obj)) {
                             auto const ij_voltage_order = (type == YBusElementType::bft);
                             process_branch_measurement(block, diag_block, rhs_block, y_branch.yff(), y_branch.yft(),
                                                        u_state, ij_voltage_order,
-                                                       measured_value.branch_from_power(obj));
+                                                       measured_values.branch_from_power(obj));
                         }
-                        if (measured_value.has_branch_to(obj)) {
+                        if (measured_values.has_branch_to(obj)) {
                             auto const ij_voltage_order = (type == YBusElementType::btf);
                             process_branch_measurement(block, diag_block, rhs_block, y_branch.ytt(), y_branch.ytf(),
-                                                       u_state, ij_voltage_order, measured_value.branch_to_power(obj));
+                                                       u_state, ij_voltage_order, measured_values.branch_to_power(obj));
                         }
-                    } else {
+                        break;
+                    }
+                    default:
                         assert(type == YBusElementType::bff || type == YBusElementType::btt);
+                        break;
                     }
                 }
 
                 // fill block with injection measurement constraints
-                if (measured_value.has_bus_injection(row)) {
+                if (measured_values.has_bus_injection(row)) {
                     auto const& yij = y_bus.admittance()[data_idx];
                     process_injection_row(block, diag_block, rhs_block, yij, u_state);
-
-                    // R_ii = -variance, only diagonal
                     if (row == col) {
-                        auto const& injection = measured_value.bus_injection(row);
-                        // assign variance to diagonal of 3x3 tensor, for asym
-                        rhs_block.tau_p() += injection.value.real();
-                        rhs_block.tau_q() += injection.value.imag();
-                        block.r_P_theta() = RealTensor<sym>{RealValue<sym>{-injection.p_variance}};
-                        block.r_Q_v() = RealTensor<sym>{RealValue<sym>{-injection.q_variance}};
+                        process_injection_diagonal(block, rhs_block, measured_values.bus_injection(row));
                     }
                 } else {
-                    // virtually remove constraints from equation
-                    // Q_ij = 0
-                    // R_ii = -1.0, only diagonal
-                    // assign -1.0 to diagonal of 3x3 tensor, for asym
                     if (row == col) {
-                        block.r_P_theta() = RealTensor<sym>{-1.0};
-                        block.r_Q_v() = RealTensor<sym>{-1.0};
+                        virtually_remove_constraints(block);
                     }
                 }
             }
         }
 
-        // loop all transpose entry for QT
-        // assign the transpose of the transpose entry of Q
         fill_qt_process_lagrange_multiplier(y_bus);
 
         // prefactorize
         sparse_solver_.prefactorize(data_gain_, perm_);
     }
 
+    void virtually_remove_constraints(NRSEGainBlock<sym>& block) {
+        // Q_ij = 0
+        // R_ii = -1.0, only diagonal
+        // assign -1.0 to diagonal of 3x3 tensor, for asym
+        block.r_P_theta() = RealTensor<sym>{-1.0};
+        block.r_Q_v() = RealTensor<sym>{-1.0};
+    }
+
+    void process_injection_diagonal(NRSEGainBlock<sym>& block, NRSERhs<sym>& rhs_block, auto const& injection) {
+        // R_ii = -variance, only diagonal
+        // assign variance to diagonal of 3x3 tensor, for asym
+        rhs_block.tau_p() += injection.value.real();
+        rhs_block.tau_q() += injection.value.imag();
+        block.r_P_theta() = RealTensor<sym>{RealValue<sym>{-injection.p_variance}};
+        block.r_Q_v() = RealTensor<sym>{RealValue<sym>{-injection.q_variance}};
+    }
+
     /**
      * @brief Processes common part of all elements to fill from an injection measurement.
      * Also includes zero injection constraint.
@@ -342,8 +346,8 @@ template <bool sym> class NewtonRaphsonSESolver {
      * @param block LHS(r, c)
      * @param diag_block LHS(r, r)
      * @param rhs_block RHS(r)
-     * @param yij
-     * @param u_state
+     * @param yij admittance of (row with injection, c)
+     * @param u_state Voltage state of iteration
      */
     void process_injection_row(NRSEGainBlock<sym>& block, NRSEGainBlock<sym>& diag_block, NRSERhs<sym>& rhs_block,
                                auto const& yij, auto const& u_state) {
@@ -361,15 +365,25 @@ template <bool sym> class NewtonRaphsonSESolver {
         rhs_block.tau_q() -= imag(f_x_complex_row);
     }
 
-    void process_shunt_measurement(NRSEGainBlock<sym>& block, NRSERhs<sym>& rhs_block, auto const& yii,
+    /**
+     * @brief Adds contribution of G(i, i) form the shunt.
+     * Note: The sign of Y_s is inverted per injection direction.
+     *
+     * @param block LHS(i, i)
+     * @param rhs_block RHS(i)
+     * @param ys shunt admittance related to the shunt measurement
+     * @param u_state Voltage state of iteration Voltage state of iteration
+     * @param measured_power measured shunt power
+     */
+    void process_shunt_measurement(NRSEGainBlock<sym>& block, NRSERhs<sym>& rhs_block, auto const& ys,
                                    auto const& u_state, auto const& measured_power) {
-        auto const hm_ui_ui_yii = hm_complex_form(yii, u_state.ui_ui_conj);
-        auto const nl_ui_ui_yii = dot(hm_ui_ui_yii, u_state.abs_ui_inv);
-        auto const f_x_complex = sum_row(hm_ui_ui_yii);
-        auto const f_x_complex_abs_ui_inv = sum_row(nl_ui_ui_yii);
+        auto const hm_ui_ui_ys = hm_complex_form(-ys, u_state.ui_ui_conj);
+        auto const nl_ui_ui_ys = dot(hm_ui_ui_ys, u_state.abs_ui_inv);
+        auto const f_x_complex = sum_row(hm_ui_ui_ys);
+        auto const f_x_complex_abs_ui_inv = sum_row(nl_ui_ui_ys);
 
-        auto jac_block = calculate_jacobian(hm_ui_ui_yii, nl_ui_ui_yii);
-        jac_block -= jacobian_diagonal_component(f_x_complex_abs_ui_inv, f_x_complex);
+        auto jac_block = calculate_jacobian(hm_ui_ui_ys, nl_ui_ui_ys);
+        jac_block += jacobian_diagonal_component(f_x_complex_abs_ui_inv, f_x_complex);
         auto const& block_F_T_k_w = transpose_multiply_weight(jac_block, measured_power);
         multiply_add_jacobian_blocks_lhs(block, block_F_T_k_w, jac_block);
         multiply_add_jacobian_blocks_rhs(rhs_block, block_F_T_k_w, measured_power, f_x_complex);
@@ -398,9 +412,9 @@ template <bool sym> class NewtonRaphsonSESolver {
      * @param block G_(r, c)
      * @param diag_block G_(r, r)
      * @param rhs_block RHS(r)
-     * @param y_xi_xi
-     * @param y_xi_mu
-     * @param u_state voltage state vector
+     * @param y_xi_xi shunt admittance near to branch measurement
+     * @param y_xi_mu admittance from the branch measurement to other bus
+     * @param u_state Voltage state of iteration voltage state vector
      * @param order bool to determine if (chi, psi) = (row, col) or (col, row)
      * @param measured_power
      */
@@ -481,31 +495,31 @@ template <bool sym> class NewtonRaphsonSESolver {
      *
      * @param block LHS(row, col), ie. LHS(row, row)
      * @param rhs_block RHS(row)
-     * @param measured_value
-     * @param bus
+     * @param measured_values
+     * @param bus bus with voltage measurement
      */
     void process_voltage_measurements(NRSEGainBlock<sym>& block, NRSERhs<sym>& rhs_block,
-                                      MeasuredValues<sym> const& measured_value, Idx const& bus) {
-        if (!measured_value.has_voltage(bus)) {
+                                      MeasuredValues<sym> const& measured_values, Idx const& bus) {
+        if (!measured_values.has_voltage(bus)) {
             return;
         }
 
         // G += 1.0 / variance
         // for 3x3 tensor, fill diagonal
-        auto const w_v = RealTensor<sym>{1.0 / measured_value.voltage_var(bus)};
-        auto const abs_measured_v = detail::cabs_or_real<sym>(measured_value.voltage(bus));
+        auto const w_v = RealTensor<sym>{1.0 / measured_values.voltage_var(bus)};
+        auto const abs_measured_v = detail::cabs_or_real<sym>(measured_values.voltage(bus));
         auto const delta_v = abs_measured_v - x_[bus].v();
 
-        auto const virtual_angle_measurement_bus = measured_value.has_voltage(math_topo_->slack_bus)
+        auto const virtual_angle_measurement_bus = measured_values.has_voltage(math_topo_->slack_bus)
                                                        ? math_topo_->slack_bus
-                                                       : measured_value.first_voltage_measurement();
+                                                       : measured_values.first_voltage_measurement();
 
         RealTensor<sym> w_theta{};
         RealValue<sym> delta_theta{};
-        if (measured_value.has_angle_measurement(bus)) {
-            delta_theta = RealValue<sym>{arg(measured_value.voltage(bus))} - RealValue<sym>{x_[bus].theta()};
+        if (measured_values.has_angle_measurement(bus)) {
+            delta_theta = RealValue<sym>{arg(measured_values.voltage(bus))} - RealValue<sym>{x_[bus].theta()};
             w_theta = RealTensor<sym>{1.0};
-        } else if (bus == virtual_angle_measurement_bus && !measured_value.has_angle()) {
+        } else if (bus == virtual_angle_measurement_bus && !measured_values.has_angle()) {
             delta_theta = arg(ComplexValue<sym>{1.0}) - RealValue<sym>{x_[bus].theta()};
             w_theta = RealTensor<sym>{1.0};
         }
@@ -643,10 +657,10 @@ template <bool sym> class NewtonRaphsonSESolver {
             // accumulate the unknown variable
             x_[bus].theta() += delta_x_rhs_[bus].theta() - RealValue<sym>{angle_offset};
             x_[bus].v() += delta_x_rhs_[bus].v();
-            if (measured_values.has_bus_injection(bus) && any_zero(measured_values.bus_injection(bus).p_variance)) {
+            if (measured_values.has_bus_injection(bus) && all_zero(measured_values.bus_injection(bus).p_variance)) {
                 x_[bus].phi_p() += delta_x_rhs_[bus].phi_p();
             }
-            if (measured_values.has_bus_injection(bus) && any_zero(measured_values.bus_injection(bus).q_variance)) {
+            if (measured_values.has_bus_injection(bus) && all_zero(measured_values.bus_injection(bus).q_variance)) {
                 x_[bus].phi_q() += delta_x_rhs_[bus].phi_q();
             }
 

power_grid_model_c/power_grid_model/include/power_grid_model/three_phase_tensor.hpp

@@ -305,6 +305,10 @@ inline auto is_inf(RealValue<false> const& value) { return is_inf(value(0)) || i
 inline auto any_zero(std::floating_point auto value) { return value == 0.0; }
 inline auto any_zero(RealValue<false> const& value) { return (value == RealValue<false>{0.0}).any(); }
 
+// all_zero
+inline auto all_zero(std::floating_point auto value) { return value == 0.0; }
+inline auto all_zero(RealValue<false> const& value) { return (value == RealValue<false>{0.0}).all(); }
+
 // update real values
 //
 // RealValue is only updated when the update value is not nan

tests/cpp_unit_tests/test_main_model_se.cpp

@@ -236,20 +236,24 @@ TEST_CASE_TEMPLATE("Test main model - state estimation", CalculationMethod, Iter
                     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
+                    CHECK(shunt_output[0].p == doctest::Approx(1800.0).epsilon(0.01).scale(1e3));
+                    CHECK(shunt_output[0].q == doctest::Approx(180.0).epsilon(0.01).scale(1e3));
+
+                    CHECK(line_output[0].p_from == doctest::Approx(1800.0).epsilon(0.01).scale(1e3));
+                    CHECK(line_output[0].q_from == doctest::Approx(180.0).epsilon(0.01).scale(1e3));
+                    CHECK(line_output[0].p_to == doctest::Approx(-1800.0).epsilon(0.01).scale(1e3));
+                    CHECK(line_output[0].q_to == doctest::Approx(-180.0).epsilon(0.01).scale(1e3));
+
+                    CHECK(power_sensor_output[0].p_residual == doctest::Approx(0.0).epsilon(0.01).scale(1e3)); // shunt
+                    CHECK(power_sensor_output[0].q_residual == doctest::Approx(0.0).epsilon(0.01).scale(1e3)); // shunt
+                    CHECK(power_sensor_output[1].p_residual ==
+                          doctest::Approx(0.0).epsilon(0.01).scale(1e3)); // branch_from
+                    CHECK(power_sensor_output[1].q_residual ==
+                          doctest::Approx(0.0).epsilon(0.01).scale(1e3)); // branch_from
+                    CHECK(power_sensor_output[2].p_residual ==
+                          doctest::Approx(0.0).epsilon(0.01).scale(1e3)); // branch_to
+                    CHECK(power_sensor_output[2].q_residual ==
+                          doctest::Approx(0.0).epsilon(0.01).scale(1e3)); // branch_to
                 }
             }
         }

tests/data/state_estimation/1os2msr-no-angle/params.json

@@ -1,14 +1,5 @@
 {
   "calculation_method":  ["iterative_linear", "newton_raphson"],
   "rtol": 1e-8,
-  "atol": 1e-8,
-  "extra_params": {
-    "newton_raphson": {
-      "experimental_features": "enabled",
-      "fail": {
-        "reason": "Experimental feature under development",
-        "raises": "PowerGridError"
-      }
-    }
-  }
+  "atol": 1e-8
 }

tests/data/state_estimation/1os2msr/params.json

@@ -4,14 +4,5 @@
   "atol": {
     "default": 1e-8,
     ".+_residual": 5e-4
-  },
-  "extra_params": {
-    "newton_raphson": {
-      "experimental_features": "enabled",
-      "fail": {
-        "reason": "Experimental feature under development",
-        "raises": "PowerGridError"
-      }
-    }
   }
 }