diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp index c80406693..53e908d66 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_linear_se_solver.hpp @@ -100,12 +100,14 @@ template class IterativeLinearSESolver { MeasuredValues 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> 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 diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp index 54b9ba9c4..e2993fd22 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/newton_raphson_se_solver.hpp @@ -166,8 +166,11 @@ template class NewtonRaphsonSESolver { // preprocess measured value sub_timer = Timer{log, LogEvent::preprocess_measured_value}; MeasuredValues 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> 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 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); }; diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp index 0890f3962..d7068401a 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/observability.hpp +++ b/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 flow_sensors; std::vector voltage_phasor_sensors; - bool is_possibly_ill_conditioned{false}; + std::vector 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 ObservabilitySensorsResult count_observability_sensors(MeasuredValues const& measured_values, MathModelTopology const& topo, @@ -35,7 +35,7 @@ ObservabilitySensorsResult count_observability_sensors(MeasuredValues const ObservabilitySensorsResult result{.flow_sensors = std::vector(y_bus_structure.row_indptr.back(), 0), .voltage_phasor_sensors = std::vector(n_bus, 0), - .is_possibly_ill_conditioned = false}; + .row_is_possibly_ill_conditioned = std::vector(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 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 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 @@ -227,7 +230,8 @@ inline ObservabilityResult observability_check(MeasuredValues 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 diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp index 15b9d8303..b39a80972 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/sparse_lu_solver.hpp +++ b/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 -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 class DenseLUFactor { // put pivot perturbation in has_pivot_perturbation in place template static void factorize_block_in_place(Eigen::MatrixBase&& matrix, BlockPerm& block_perm, - double perturb_threshold, bool use_pivot_perturbation, - bool& has_pivot_perturbation) + double perturb_threshold, + std::pair> possibly_ill_conditioned_pivots, + Idx block_node_idx, bool& has_pivot_perturbation, double matrix_norm) requires(std::same_as && rk2_tensor && (Derived::RowsAtCompileTime == size) && (Derived::ColsAtCompileTime == size)) { @@ -84,7 +85,7 @@ template 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 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 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 SparseLUSolver { prefactorize_and_solve(std::vector& data, // matrix data, factorize in-place BlockPermArray& block_perm_array, // pre-allocated permutation array, will be overwritten std::vector const& rhs, std::vector& x, - bool use_pivot_perturbation = false) { - prefactorize(data, block_perm_array, use_pivot_perturbation); + std::pair> 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 const&)data, block_perm_array, rhs, x); } @@ -239,9 +246,10 @@ template class SparseLUSolver { // fill-ins should be pre-allocated with zero // block permutation array should be pre-allocated void prefactorize(std::vector& data, BlockPermArray& block_perm_array, - bool use_pivot_perturbation = false) { + std::pair> possibly_ill_conditioned_pivots = {false, {}}) { reset_matrix_cache(); - if (use_pivot_perturbation) { + if (possibly_ill_conditioned_pivots.first) { + assert(static_cast(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 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 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 { diff --git a/tests/cpp_unit_tests/test_sparse_lu_solver.cpp b/tests/cpp_unit_tests/test_sparse_lu_solver.cpp index d1915de71..978a74abd 100644 --- a/tests/cpp_unit_tests/test_sparse_lu_solver.cpp +++ b/tests/cpp_unit_tests/test_sparse_lu_solver.cpp @@ -7,6 +7,8 @@ #include +#include + namespace power_grid_model::math_solver { namespace { using lu_trait_double = math_solver::sparse_lu_entry_trait; @@ -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 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 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); } diff --git a/tests/data/state_estimation/ill-conditioned-by-line-meshed/params.json b/tests/data/state_estimation/ill-conditioned-by-line-meshed/params.json index 3358a01ad..9e558925e 100644 --- a/tests/data/state_estimation/ill-conditioned-by-line-meshed/params.json +++ b/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" } } \ No newline at end of file diff --git a/tests/data/state_estimation/ill-conditioned-by-link-meshed/params.json b/tests/data/state_estimation/ill-conditioned-by-link-meshed/params.json index 6dd5b9a8d..88cb11fee 100644 --- a/tests/data/state_estimation/ill-conditioned-by-link-meshed/params.json +++ b/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" } } \ No newline at end of file