Sparse matrix solver: Trigger pivot perturbation based on each possibly ill-conditioned one

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

@@ -100,12 +100,14 @@ template <symmetry_tag sym_type> class IterativeLinearSESolver {
         MeasuredValues<sym> const measured_values{y_bus.shared_topology(), input};
         auto const observability_result =
             observability_check(measured_values, y_bus.math_topology(), y_bus.y_bus_structure());
+        std::pair<bool, std::vector<int8_t>> const possibly_ill_conditioned_pivots{
+            observability_result.use_perturbation(), observability_result.row_is_possibly_ill_conditioned};
 
         // prepare matrix
         sub_timer = Timer{log, LogEvent::prepare_matrix_including_prefactorization};
         prepare_matrix(y_bus, measured_values);
         // prefactorize
-        sparse_solver_.prefactorize(data_gain_, perm_, observability_result.use_perturbation());
+        sparse_solver_.prefactorize(data_gain_, perm_, possibly_ill_conditioned_pivots);
 
         // initialize voltage with initial angle
         sub_timer = Timer{log, LogEvent::initialize_voltages}; // TODO(mgovers): make scoped subtimers

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

@@ -166,8 +166,11 @@ template <symmetry_tag sym_type> class NewtonRaphsonSESolver {
         // preprocess measured value
         sub_timer = Timer{log, LogEvent::preprocess_measured_value};
         MeasuredValues<sym> const measured_values{y_bus.shared_topology(), input};
+
         auto const observability_result =
             observability_check(measured_values, y_bus.math_topology(), y_bus.y_bus_structure());
+        std::pair<bool, std::vector<int8_t>> const possibly_ill_conditioned_pivots{
+            observability_result.use_perturbation(), observability_result.row_is_possibly_ill_conditioned};
 
         // initialize voltage with initial angle
         sub_timer = Timer{log, LogEvent::initialize_voltages};
@@ -184,7 +187,7 @@ template <symmetry_tag sym_type> class NewtonRaphsonSESolver {
             // solve with prefactorization
             sub_timer = Timer{log, LogEvent::solve_sparse_linear_equation};
             sparse_solver_.prefactorize_and_solve(data_gain_, perm_, delta_x_rhs_, delta_x_rhs_,
-                                                  observability_result.use_perturbation());
+                                                  possibly_ill_conditioned_pivots);
             sub_timer = Timer{log, LogEvent::iterate_unknown};
             max_dev = iterate_unknown(output.u, measured_values);
         };

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp

@@ -15,7 +15,7 @@ namespace detail {
 struct ObservabilitySensorsResult {
     std::vector<int8_t> flow_sensors;
     std::vector<int8_t> voltage_phasor_sensors;
-    bool is_possibly_ill_conditioned{false};
+    std::vector<int8_t> row_is_possibly_ill_conditioned;
 };
 
 // count flow and voltage phasor sensors for the observability check
@@ -26,7 +26,7 @@ struct ObservabilitySensorsResult {
 // return a ObservabilitySensorsResult struct with
 //      a vector of flow sensor count
 //      a vector of voltage phasor sensor count
-//      a boolean indicating if the system is possibly ill-conditioned
+//      a vector of boolean indicating if each row is possibly ill-conditioned
 template <symmetry_tag sym>
 ObservabilitySensorsResult count_observability_sensors(MeasuredValues<sym> const& measured_values,
                                                        MathModelTopology const& topo,
@@ -35,7 +35,7 @@ ObservabilitySensorsResult count_observability_sensors(MeasuredValues<sym> const
 
     ObservabilitySensorsResult result{.flow_sensors = std::vector<int8_t>(y_bus_structure.row_indptr.back(), 0),
                                       .voltage_phasor_sensors = std::vector<int8_t>(n_bus, 0),
-                                      .is_possibly_ill_conditioned = false};
+                                      .row_is_possibly_ill_conditioned = std::vector<int8_t>(n_bus, 0)};
 
     auto has_flow_sensor = [&measured_values](Idx branch) {
         return measured_values.has_branch_from_power(branch) || measured_values.has_branch_to_power(branch) ||
@@ -81,7 +81,7 @@ ObservabilitySensorsResult count_observability_sensors(MeasuredValues<sym> const
         }
         // the system could be ill-conditioned if there is no flow sensor for one bus, except the last bus
         if (!has_at_least_one_sensor && row != n_bus - 1) {
-            result.is_possibly_ill_conditioned = true;
+            result.row_is_possibly_ill_conditioned[row] = 1;
         }
     }
     return result;
@@ -195,8 +195,11 @@ inline bool sufficient_observability_condition(YBusStructure const& y_bus_struct
 
 struct ObservabilityResult {
     bool is_observable{false};
-    bool is_possibly_ill_conditioned{false};
-    constexpr bool use_perturbation() const { return is_possibly_ill_conditioned && is_observable; }
+    std::vector<int8_t> row_is_possibly_ill_conditioned{};
+    bool use_perturbation() const {
+        return std::ranges::any_of(row_is_possibly_ill_conditioned, [](int8_t element) { return element == 1; }) &&
+               is_observable;
+    }
 };
 
 template <symmetry_tag sym>
@@ -227,7 +230,8 @@ inline ObservabilityResult observability_check(MeasuredValues<sym> const& measur
     }
     //  ToDo(JGuo): meshed network will require a different treatment
     return ObservabilityResult{.is_observable = is_necessary_condition_met && is_sufficient_condition_met,
-                               .is_possibly_ill_conditioned = observability_sensors.is_possibly_ill_conditioned};
+                               .row_is_possibly_ill_conditioned =
+                                   observability_sensors.row_is_possibly_ill_conditioned};
 }
 
 } // namespace power_grid_model::math_solver

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp

@@ -24,9 +24,9 @@ constexpr double epsilon_sqrt = 1.49011745e-8;      // sqrt(epsilon), sqrt does
 // it will get modified if perturbation happens
 // also has_pivot_perturbation will be updated if perturbation happens
 template <scalar_value Scalar>
-inline void perturb_pivot_if_needed(double perturb_threshold, Scalar& value, double& abs_value,
+inline void perturb_pivot_if_needed(int8_t pertub_pivot, double perturb_threshold, Scalar& value, double& abs_value,
                                     bool& has_pivot_perturbation) {
-    if (abs_value < perturb_threshold) {
+    if (pertub_pivot == 1 && abs_value < perturb_threshold) {
         Scalar const scale = (abs_value == 0.0) ? Scalar{1.0} : (value / abs_value);
         value = scale * perturb_threshold;
         has_pivot_perturbation = true;
@@ -60,8 +60,9 @@ template <rk2_tensor Matrix> class DenseLUFactor {
     // put pivot perturbation in has_pivot_perturbation in place
     template <class Derived>
     static void factorize_block_in_place(Eigen::MatrixBase<Derived>&& matrix, BlockPerm& block_perm,
-                                         double perturb_threshold, bool use_pivot_perturbation,
-                                         bool& has_pivot_perturbation)
+                                         double perturb_threshold,
+                                         std::pair<bool, std::vector<int8_t>> possibly_ill_conditioned_pivots,
+                                         Idx block_node_idx, bool& has_pivot_perturbation, double matrix_norm)
         requires(std::same_as<typename Derived::Scalar, Scalar> && rk2_tensor<Derived> &&
                  (Derived::RowsAtCompileTime == size) && (Derived::ColsAtCompileTime == size))
     {
@@ -84,7 +85,7 @@ template <rk2_tensor Matrix> class DenseLUFactor {
             assert(col_biggest_eigen + pivot < size);
 
             // check absolute singular matrix
-            if (biggest_score == 0.0 && !use_pivot_perturbation) {
+            if (biggest_score == 0.0 && !possibly_ill_conditioned_pivots.first) {
                 // pivot perturbation not possible, cannot proceed
                 // set identity permutation and break the loop
                 for (int8_t remaining_rows_cols = pivot; remaining_rows_cols != size; ++remaining_rows_cols) {
@@ -96,8 +97,11 @@ template <rk2_tensor Matrix> class DenseLUFactor {
 
             // perturb pivot if needed
             double abs_pivot = sqrt(biggest_score);
-            perturb_pivot_if_needed(perturb_threshold, matrix(row_biggest, col_biggest), abs_pivot,
-                                    has_pivot_perturbation);
+            int8_t const possibly_ill_conditioned_pivot = possibly_ill_conditioned_pivots.first
+                                                              ? possibly_ill_conditioned_pivots.second[block_node_idx]
+                                                              : int8_t{0};
+            perturb_pivot_if_needed(possibly_ill_conditioned_pivot, perturb_threshold, matrix(row_biggest, col_biggest),
+                                    abs_pivot, has_pivot_perturbation);
             max_pivot = std::max(max_pivot, abs_pivot);
 
             // swap rows and columns
@@ -132,9 +136,12 @@ template <rk2_tensor Matrix> class DenseLUFactor {
 
         // throw SparseMatrixError if the matrix is ill-conditioned
         // only check condition number if pivot perturbation is not used
-        double const pivot_threshold = has_pivot_perturbation ? 0.0 : epsilon * max_pivot;
+        double const absolute_threshold = epsilon;
+        double const relative_threshold = has_pivot_perturbation ? 0.0 : epsilon * matrix_norm;
         for (int8_t pivot = 0; pivot != size; ++pivot) {
-            if (cabs(matrix(pivot, pivot)) < pivot_threshold || !is_normal(matrix(pivot, pivot))) {
+            auto const pivot_abs_v = cabs(matrix(pivot, pivot));
+            if (pivot_abs_v < relative_threshold || pivot_abs_v < absolute_threshold ||
+                !is_normal(matrix(pivot, pivot))) {
                 throw SparseMatrixError{}; // can not specify error code
             }
         }
@@ -214,8 +221,8 @@ template <class Tensor, class RHSVector, class XVector> class SparseLUSolver {
     prefactorize_and_solve(std::vector<Tensor>& data,        // matrix data, factorize in-place
                            BlockPermArray& block_perm_array, // pre-allocated permutation array, will be overwritten
                            std::vector<RHSVector> const& rhs, std::vector<XVector>& x,
-                           bool use_pivot_perturbation = false) {
-        prefactorize(data, block_perm_array, use_pivot_perturbation);
+                           std::pair<bool, std::vector<int8_t>> possibly_ill_conditioned_pivots = {false, {}}) {
+        prefactorize(data, block_perm_array, possibly_ill_conditioned_pivots);
         // call solve with const method
         solve_with_prefactorized_matrix((std::vector<Tensor> const&)data, block_perm_array, rhs, x);
     }
@@ -239,9 +246,10 @@ template <class Tensor, class RHSVector, class XVector> class SparseLUSolver {
     // fill-ins should be pre-allocated with zero
     // block permutation array should be pre-allocated
     void prefactorize(std::vector<Tensor>& data, BlockPermArray& block_perm_array,
-                      bool use_pivot_perturbation = false) {
+                      std::pair<bool, std::vector<int8_t>> possibly_ill_conditioned_pivots = {false, {}}) {
         reset_matrix_cache();
-        if (use_pivot_perturbation) {
+        if (possibly_ill_conditioned_pivots.first) {
+            assert(static_cast<Idx>(possibly_ill_conditioned_pivots.second.size()) == size_);
             initialize_pivot_perturbation(data);
         }
         double const perturb_threshold = epsilon_perturbation * matrix_norm_;
@@ -268,15 +276,17 @@ template <class Tensor, class RHSVector, class XVector> class SparseLUSolver {
                     // use machine precision by default
                     // record block permutation
                     LUFactor::factorize_block_in_place(lu_matrix[pivot_idx].matrix(), block_perm_array[pivot_row_col],
-                                                       perturb_threshold, use_pivot_perturbation,
-                                                       has_pivot_perturbation_);
+                                                       perturb_threshold, possibly_ill_conditioned_pivots,
+                                                       pivot_row_col, // Pass the node index
+                                                       has_pivot_perturbation_, matrix_norm_);
                     return block_perm_array[pivot_row_col];
                 } else {
-                    if (use_pivot_perturbation) {
+                    if (possibly_ill_conditioned_pivots.first) {
                         // use machine precision by default
                         // record pivot perturbation
                         double abs_pivot = cabs(lu_matrix[pivot_idx]);
-                        perturb_pivot_if_needed(perturb_threshold, lu_matrix[pivot_idx], abs_pivot,
+                        perturb_pivot_if_needed(possibly_ill_conditioned_pivots.second[pivot_row_col],
+                                                perturb_threshold, lu_matrix[pivot_idx], abs_pivot,
                                                 has_pivot_perturbation_);
                     }
                     if (!is_normal(lu_matrix[pivot_idx])) {
@@ -539,19 +549,14 @@ template <class Tensor, class RHSVector, class XVector> class SparseLUSolver {
         // calculate the block-wise non-diagonal infinite norm of the matrix
         // that is:
         // 1. calculate the infinite norm of each individual block
-        // 2. sum all norms of the blocks per row, except the diagonal block
+        // 2. sum all norms of the blocks per row,
         // 3. take the maximum of all the sums
         matrix_norm_ = 0.0;
         auto const& row_indptr = *row_indptr_;
-        auto const& col_indices = *col_indices_;
         for (Idx row = 0; row != size_; ++row) {
             // calculate the sum of the norms of the blocks in the row
             double row_norm = 0.0;
             for (Idx idx = row_indptr[row]; idx != row_indptr[row + 1]; ++idx) {
-                // skip diagonal
-                if (col_indices[idx] == row) {
-                    continue;
-                }
                 if constexpr (is_block) {
                     row_norm += cabs(data[idx]).rowwise().sum().maxCoeff();
                 } else {

tests/cpp_unit_tests/test_sparse_lu_solver.cpp

@@ -7,6 +7,8 @@
 
 #include <doctest/doctest.h>
 
+#include <bitset>
+
 namespace power_grid_model::math_solver {
 namespace {
 using lu_trait_double = math_solver::sparse_lu_entry_trait<double, double, double>;
@@ -158,11 +160,32 @@ TEST_CASE("LU solver with ill-conditioned system") {
         Idx perm{};
 
         SUBCASE("Error without perturbation") {
-            CHECK_THROWS_AS(solver.prefactorize(data, perm, false), SparseMatrixError);
+            CHECK_THROWS_AS(solver.prefactorize(data, perm, {false, {}}), SparseMatrixError);
+        }
+
+        // first row is ill conditioned, others are fine: the cases where the first row is not perturbed should fail
+        SUBCASE("Error with wrong perturbation") {
+            // generate all possible combinations with 4 possibly ill conditioned rows (16 total)
+            for (int8_t combination = 0; combination < 16; ++combination) {
+                std::string binary_possibly_ill_conditioned_rows = std::bitset<4>(combination).to_string();
+                // we test combinations where the first row is not perturbed
+                if (binary_possibly_ill_conditioned_rows[0] == '1') {
+                    continue;
+                }
+                std::vector<int8_t> possibly_ill_conditioned_rows(4);
+                std::ranges::transform(binary_possibly_ill_conditioned_rows, possibly_ill_conditioned_rows.begin(),
+                                       [](char c) -> int8_t { return c == '1' ? 1 : 0; });
+                CAPTURE(combination);
+                CAPTURE(possibly_ill_conditioned_rows);
+                // first row is not perturbed: should fail
+                CHECK_THROWS_AS(solver.prefactorize(data, perm, {true, possibly_ill_conditioned_rows}),
+                                SparseMatrixError);
+            }
         }
 
-        SUBCASE("Success with perturbation") {
-            CHECK_NOTHROW(solver.prefactorize(data, perm, true));
+        SUBCASE("Success with correct perturbation") {
+            // Only the first row is ill conditioned: possibly trigger perturbation on the correct row
+            CHECK_NOTHROW(solver.prefactorize(data, perm, {true, {1, 0, 0, 0}}));
             solver.solve_with_prefactorized_matrix(data, perm, rhs, x);
             check_result(x, x_ref);
         }
@@ -186,11 +209,17 @@ TEST_CASE("LU solver with ill-conditioned system") {
         SparseLUSolver<Tensor, Array, Array> solver{row_indptr, col_indices, diag_lu};
 
         SUBCASE("Error without perturbation") {
-            CHECK_THROWS_AS(solver.prefactorize(data, block_perm, false), SparseMatrixError);
+            CHECK_THROWS_AS(solver.prefactorize(data, block_perm, {false, {}}), SparseMatrixError);
+        }
+
+        SUBCASE("Error with wrong perturbation") {
+            // first block is ill conditioned, second is fine: trigger perturbation on the wrong block
+            CHECK_THROWS_AS(solver.prefactorize(data, block_perm, {true, {0, 1}}), SparseMatrixError);
         }
 
-        SUBCASE("Success with perturbation") {
-            CHECK_NOTHROW(solver.prefactorize(data, block_perm, true));
+        SUBCASE("Success with correct perturbation") {
+            // first block is ill conditioned, second is fine: trigger perturbation on the correct block
+            CHECK_NOTHROW(solver.prefactorize(data, block_perm, {true, {1, 0}}));
             solver.solve_with_prefactorized_matrix(data, block_perm, rhs, x);
             check_result(x, x_ref);
         }

tests/data/state_estimation/ill-conditioned-by-line-meshed/params.json

@@ -3,9 +3,5 @@
   "rtol": 1e-8,
   "atol": {
     "default": 1e-8
-  },
-  "raises": {
-    "raises": "SparseMatrixError",
-    "reason": "Observability check for meshed grids is not yet implemented. Also, pivot perturbation wouldn't be triggered as Matrix norm is smaller than the max pivot. See https://github.com/PowerGridModel/power-grid-model/pull/1118#issuecomment-3269339573"
   }
 }

tests/data/state_estimation/ill-conditioned-by-link-meshed/params.json

@@ -6,9 +6,5 @@
   "rtol": 1e-8,
   "atol": {
     "default": 1e-8
-  },
-  "raises": {
-    "raises": "SparseMatrixError",
-    "reason": "Observability check for meshed grids is not yet implemented. Also, pivot perturbation wouldn't be triggered as Matrix norm is smaller than the max pivot. See https://github.com/PowerGridModel/power-grid-model/pull/1118#issuecomment-3269339573"
   }
 }