Merge pull request #1030 from PowerGridModel/feature/nrse-global-current-sensor

NRSE current sensor: Global angle current sensor

Author: nitbharambe

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

@@ -122,6 +122,8 @@ template <symmetry_tag sym_type> class NewtonRaphsonSESolver {
         auto const& abs_u_psi_inv(Order ij_voltage_order) const {
             return ij_voltage_order == Order::row_major ? abs_uj_inv : abs_ui_inv;
         }
+        auto const& u_chi(Order ij_voltage_order) const { return ij_voltage_order == Order::row_major ? ui : uj; }
+        auto const& u_psi(Order ij_voltage_order) const { return ij_voltage_order == Order::row_major ? uj : ui; }
     };
 
     struct NRSEJacobian {
@@ -336,7 +338,10 @@ template <symmetry_tag sym_type> class NewtonRaphsonSESolver {
                                                                          measurement);
                                 break;
                             case AngleMeasurementType::global_angle:
-                                throw SparseMatrixError{};
+                                process_branch_global_current_measurement(block, diag_block, rhs_block, y_branch.yff(),
+                                                                          y_branch.yft(), u_state, ij_voltage_order,
+                                                                          measurement);
+                                break;
                             default:
                                 assert(measurement.angle_measurement_type == AngleMeasurementType::local_angle ||
                                        measurement.angle_measurement_type == AngleMeasurementType::global_angle);
@@ -353,7 +358,10 @@ template <symmetry_tag sym_type> class NewtonRaphsonSESolver {
                                                                          measurement);
                                 break;
                             case AngleMeasurementType::global_angle:
-                                throw SparseMatrixError{};
+                                process_branch_global_current_measurement(block, diag_block, rhs_block, y_branch.ytt(),
+                                                                          y_branch.ytf(), u_state, ij_voltage_order,
+                                                                          measurement);
+                                break;
                             default:
                                 assert(measurement.angle_measurement_type == AngleMeasurementType::local_angle ||
                                        measurement.angle_measurement_type == AngleMeasurementType::global_angle);
@@ -503,8 +511,8 @@ template <symmetry_tag sym_type> class NewtonRaphsonSESolver {
         }
     }
 
-    /// @brief Adds contribution of branch current measurements. Logic is similar to powermeasurements, but with an
-    /// additional voltage division as per mathematical workout.
+    /// @brief Adds contribution of local angle branch current measurements. Logic is similar to power measurements, but
+    /// with an additional voltage division as per mathematical workout.
     void process_branch_local_current_measurement(NRSEGainBlock<sym>& block, NRSEGainBlock<sym>& diag_block,
                                                   NRSERhs<sym>& rhs_block, auto const& y_xi_xi, auto const& y_xi_mu,
                                                   auto const& u_state, Order const order,
@@ -533,6 +541,29 @@ template <symmetry_tag sym_type> class NewtonRaphsonSESolver {
         }
     }
 
+    /// @brief Adds contribution of global angle branch current measurements.
+    void process_branch_global_current_measurement(NRSEGainBlock<sym>& block, NRSEGainBlock<sym>& diag_block,
+                                                   NRSERhs<sym>& rhs_block, auto const& y_xi_xi, auto const& y_xi_mu,
+                                                   auto const& u_state, Order const order,
+                                                   CurrentSensorCalcParam<sym> const& current_sensor) {
+        ComplexTensor<sym> const current_chi_chi = dot(y_xi_xi, ComplexDiagonalTensor<sym>{u_state.u_chi(order)});
+        ComplexTensor<sym> const current_chi_psi = dot(y_xi_mu, ComplexDiagonalTensor<sym>{u_state.u_psi(order)});
+        ComplexTensor<sym> const current_chi_chi_v_inv = dot(current_chi_chi, u_state.abs_u_chi_inv(order));
+        ComplexTensor<sym> const current_chi_psi_v_inv = dot(current_chi_psi, u_state.abs_u_psi_inv(order));
+        ComplexValue<sym> const f_x_complex = sum_row(current_chi_chi + current_chi_psi);
+
+        auto const block_rr_or_cc = calculate_jacobian(-current_chi_chi, current_chi_chi_v_inv);
+        auto const block_rc_or_cr = calculate_jacobian(-current_chi_psi, current_chi_psi_v_inv);
+
+        if (order == Order::row_major) {
+            multiply_add_branch_blocks(block, diag_block, rhs_block, block_rr_or_cc, block_rc_or_cr,
+                                       current_sensor.measurement, f_x_complex);
+        } else {
+            multiply_add_branch_blocks(block, diag_block, rhs_block, block_rc_or_cr, block_rr_or_cc,
+                                       current_sensor.measurement, f_x_complex);
+        }
+    }
+
     void multiply_add_branch_blocks(NRSEGainBlock<sym>& block, NRSEGainBlock<sym>& diag_block, NRSERhs<sym>& rhs_block,
                                     auto& left_block, auto const& right_block,
                                     DecomposedComplexRandVar<sym> const& measured_flow, auto const& f_x_complex) {
@@ -709,6 +740,10 @@ template <symmetry_tag sym_type> class NewtonRaphsonSESolver {
     /// @brief Construct the F_k(u1, u2, y12) block using helper function of hnml complex form
     /// The 4 members are H, N, M, L in the order.
     ///
+    /// Function is also reused for global current sensor.
+    /// The notation in that case is as per current and current * V^-1 instead of HNML.
+    ///
+    ///
     /// @param hm_complex hm_complex
     /// @param nl_complex hm_complex / abs(u2)
     /// @return  F_k(u1, u2, y12)