diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp index 52355f576..79f1b886a 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/link_solver.hpp @@ -319,4 +319,54 @@ inline std::vector> set_projection_system(Idx free_in return projection_system; }; +inline void naive_gauss_elimination(std::vector>& system) { + auto const system_size = narrow_cast(std::ssize(system)); + + // we skip pivoting since the matrix system is mostly diagonally dominant + // and given the real power systems topology of super nodes, this should not introduce + // any numerical instabilities in practice + + // forward elimination + for (Idx const column : IdxRange{system_size}) { + auto const pivot = system[column][column]; + for (Idx const row : IdxRange{column + Idx{1}, system_size}) { + auto& row_col_value = system[row][column]; + row_col_value /= -pivot; + for (Idx const column_part : IdxRange{column + Idx{1}, system_size + Idx{1}}) { + system[row][column_part] += row_col_value * system[column][column_part]; + } + } + } + + // backward substitution + system[system_size - 1][system_size] /= system[system_size - 1][system_size - 1]; + for (Idx const row : IdxRange{system_size - 1} | std::views::reverse) { + auto element_sum = DoubleComplex{}; + for (Idx const column : IdxRange{row + Idx{1}, system_size}) { + element_sum -= system[row][column] * system[column][system_size]; + } + system[row][system_size] = (system[row][system_size] + element_sum) / system[row][row]; + } +}; + +inline std::vector compute_internal_loads(SolutionSet const& solution_set, + std::vector> const& system) { + auto const number_of_rows = narrow_cast(solution_set.extended_rhs.size()); + auto const number_of_columns = narrow_cast(system.size()); + std::vector internal_loads(number_of_rows); + + for (auto const row : IdxRange{number_of_rows}) { + internal_loads[row] = solution_set.extended_rhs[row]; + auto sum_value = DoubleComplex{}; + for (auto const column : IdxRange{number_of_columns}) { + solution_set.dfs_matrix.get_value(row, column).transform([&sum_value, &system, column](IntS value) { + sum_value += static_cast(value) * system[column].back(); + return value; + }); + } + internal_loads[row] -= sum_value; + } + + return internal_loads; +}; } // namespace power_grid_model::link_solver::detail diff --git a/tests/cpp_unit_tests/test_link_solver.cpp b/tests/cpp_unit_tests/test_link_solver.cpp index f38e9f963..a0ceba0e1 100644 --- a/tests/cpp_unit_tests/test_link_solver.cpp +++ b/tests/cpp_unit_tests/test_link_solver.cpp @@ -4,10 +4,13 @@ #include #include +#include +#include #include #include +#include #include #include #include @@ -519,5 +522,363 @@ TEST_CASE("Test the link solver algorithm") { CHECK(projection_system == test_system); } } + + SUBCASE("Testing the gauss elimination routine") { + auto compare_systems = [](std::vector> const& solution, + std::vector> const& reference, Idx col_number, + Idx row_number, double tolerance) { + for (Idx const col : IdxRange{col_number}) { + for (Idx const row : IdxRange{row_number}) { + CHECK(solution[row][col].real() == doctest::Approx(reference[row][col].real()).epsilon(tolerance)); + CHECK(solution[row][col].imag() == doctest::Approx(reference[row][col].imag()).epsilon(tolerance)); + } + } + }; + + auto compare_vectors = [](std::vector>& solution, + std::vector& reference, double tolerance) { + auto const solution_size = narrow_cast(solution.size()); + // we only test against the last column of the solution as this lambda is intended for large edge test cases + // this makes our lifes easier + for (Idx const row : IdxRange{solution_size}) { + CHECK(solution[row][solution_size].real() == doctest::Approx(reference[row].real()).epsilon(tolerance)); + CHECK(solution[row][solution_size].imag() == doctest::Approx(reference[row].imag()).epsilon(tolerance)); + } + }; + + SUBCASE("Linear system of the complex case") { + std::vector> system = {{{3, 0}, {1, 0}, {1, 0}, {2, 2}}, + {{1, 0}, {3, 0}, {-1, 0}, {3, 3}}, + {{1, 0}, {-1, 0}, {4, 0}, {-1, -1}}}; + naive_gauss_elimination(system); + std::vector> const test_system = { + {{3, 0}, {1, 0}, {1, 0}, {0.458333, 0.458333}}, + {{-0.333333, 0}, {2.66667, 0}, {-1.33333, 0}, {0.791667, 0.791667}}, + {{-0.333333, 0}, {0.5, -0}, {3, 0}, {-0.166667, -0.166667}}}; + + // the tolerance is set to 1e-5 because that's the test system tolerance + compare_systems(system, test_system, Idx{4}, Idx{3}, 1e-5); + } + + SUBCASE("A system that consisng of a 15 X 15 matrix with externally randomly generated elements") { + std::vector> system = {{{7, 0}, + {4, 0}, + {7, 0}, + {2, 0}, + {14, 0}, + {5, 0}, + {2, 0}, + {4, 0}, + {14, 0}, + {2, 0}, + {6, 0}, + {7, 0}, + {7, 0}, + {0, 0}, + {3, 0}, + {1, 0}}, + {{7, 0}, + {11, 0}, + {9, 0}, + {7, 0}, + {0, 0}, + {12, 0}, + {1, 0}, + {14, 0}, + {7, 0}, + {12, 0}, + {0, 0}, + {4, 0}, + {4, 0}, + {14, 0}, + {13, 0}, + {1, 0}}, + {{6, 0}, + {0, 0}, + {13, 0}, + {0, 0}, + {0, 0}, + {2, 0}, + {4, 0}, + {8, 0}, + {8, 0}, + {14, 0}, + {9, 0}, + {8, 0}, + {0, 0}, + {3, 0}, + {11, 0}, + {1, 0}}, + {{9, 0}, + {2, 0}, + {2, 0}, + {14, 0}, + {5, 0}, + {4, 0}, + {14, 0}, + {7, 0}, + {4, 0}, + {8, 0}, + {4, 0}, + {5, 0}, + {11, 0}, + {10, 0}, + {4, 0}, + {1, 0}}, + {{10, 0}, + {7, 0}, + {12, 0}, + {12, 0}, + {12, 0}, + {7, 0}, + {13, 0}, + {7, 0}, + {14, 0}, + {2, 0}, + {14, 0}, + {5, 0}, + {2, 0}, + {1, 0}, + {0, 0}, + {1, 0}}, + {{1, 0}, + {3, 0}, + {0, 0}, + {7, 0}, + {3, 0}, + {14, 0}, + {11, 0}, + {5, 0}, + {6, 0}, + {11, 0}, + {3, 0}, + {7, 0}, + {0, 0}, + {12, 0}, + {1, 0}, + {1, 0}}, + {{2, 0}, + {11, 0}, + {9, 0}, + {2, 0}, + {0, 0}, + {3, 0}, + {0, 0}, + {8, 0}, + {0, 0}, + {12, 0}, + {8, 0}, + {5, 0}, + {14, 0}, + {10, 0}, + {4, 0}, + {1, 0}}, + {{9, 0}, + {12, 0}, + {2, 0}, + {13, 0}, + {0, 0}, + {11, 0}, + {8, 0}, + {1, 0}, + {1, 0}, + {13, 0}, + {2, 0}, + {8, 0}, + {10, 0}, + {2, 0}, + {13, 0}, + {1, 0}}, + {{14, 0}, + {9, 0}, + {13, 0}, + {13, 0}, + {13, 0}, + {12, 0}, + {11, 0}, + {2, 0}, + {0, 0}, + {3, 0}, + {11, 0}, + {3, 0}, + {6, 0}, + {6, 0}, + {13, 0}, + {1, 0}}, + {{3, 0}, + {14, 0}, + {4, 0}, + {7, 0}, + {10, 0}, + {14, 0}, + {6, 0}, + {13, 0}, + {11, 0}, + {12, 0}, + {6, 0}, + {7, 0}, + {14, 0}, + {12, 0}, + {0, 0}, + {1, 0}}, + {{13, 0}, + {11, 0}, + {9, 0}, + {14, 0}, + {14, 0}, + {14, 0}, + {12, 0}, + {13, 0}, + {1, 0}, + {10, 0}, + {8, 0}, + {8, 0}, + {11, 0}, + {14, 0}, + {14, 0}, + {1, 0}}, + {{10, 0}, + {14, 0}, + {9, 0}, + {3, 0}, + {7, 0}, + {11, 0}, + {8, 0}, + {8, 0}, + {11, 0}, + {3, 0}, + {0, 0}, + {6, 0}, + {12, 0}, + {5, 0}, + {2, 0}, + {1, 0}}, + {{7, 0}, + {13, 0}, + {8, 0}, + {7, 0}, + {0, 0}, + {12, 0}, + {4, 0}, + {9, 0}, + {8, 0}, + {13, 0}, + {3, 0}, + {8, 0}, + {0, 0}, + {5, 0}, + {2, 0}, + {1, 0}}, + {{2, 0}, + {2, 0}, + {3, 0}, + {10, 0}, + {10, 0}, + {14, 0}, + {0, 0}, + {5, 0}, + {7, 0}, + {5, 0}, + {6, 0}, + {10, 0}, + {8, 0}, + {11, 0}, + {6, 0}, + {1, 0}}, + {{9, 0}, + {6, 0}, + {0, 0}, + {2, 0}, + {2, 0}, + {12, 0}, + {13, 0}, + {13, 0}, + {12, 0}, + {9, 0}, + {8, 0}, + {5, 0}, + {12, 0}, + {9, 0}, + {7, 0}, + {1, 0}}}; + + std::vector test_solution = { + {0.05461404, 0}, {0.03584441, 0}, {-0.00895461, 0}, {-0.00979037, 0}, {-0.01083266, 0}, + {-0.03845678, 0}, {-0.00652489, 0}, {-0.08356931, 0}, {0.05730963, 0}, {0.01390954, 0}, + {0.026622, 0}, {0.04469859, 0}, {-0.00946348, 0}, {0.08945877, 0}, {0.00377452, 0}}; + + naive_gauss_elimination(system); + + // error tolerance is increased as this is a stress test + // the test system solution now includes for significant digits for this reason + // this is important because we are using naive_gauss_elimination, which skips pivoting + // so this test makes sure that assumption holds + compare_vectors(system, test_solution, 1e-7); + } + } + + SUBCASE("Testing the compute_internal_loads_routine") { + auto generate_input_solution_set = [](std::span data, std::span rows, + std::span cols, Idx row_number, Idx col_number) { + SolutionSet solution_set{}; + solution_set.dfs_matrix.prepare(row_number, col_number); + for (auto idx = size_t{0}; idx < data.size(); ++idx) { + solution_set.dfs_matrix.set_value(data[idx], rows[idx], cols[idx]); + } + return solution_set; + }; + + auto compare_vectors = [](std::vector& loads, std::vector& test_loads, + double tolerance) { + auto const loads_size = narrow_cast(loads.size()); + auto const test_loads_size = narrow_cast(test_loads.size()); + REQUIRE(loads_size == test_loads_size); + + for (Idx const idx : IdxRange(loads_size)) { + CHECK(loads[idx].real() == doctest::Approx(test_loads[idx].real()).epsilon(tolerance)); + CHECK(loads[idx].imag() == doctest::Approx(test_loads[idx].imag()).epsilon(tolerance)); + } + }; + + SUBCASE("Complex case with complex loads") { + std::vector data = {1, 1, 1, -1, -1, -1, -1, -1, 1, -1}; + std::vector row = {0, 0, 1, 1, 2, 2, 3, 4, 5, 6}; + std::vector col = {0, 2, 1, 2, 0, 1, 0, 1, 2, 2}; + + auto solution_set = generate_input_solution_set(data, row, col, Idx{7}, Idx{3}); + + solution_set.extended_rhs = {{0, 0}, {1, 1}, {-2, -2}, {0, 0}, {0, 0}, {0, 0}, {0, 0}}; + + std::vector> const test_system = { + {{3, 0}, {1, 0}, {1, 0}, {0.458333, 0.458333}}, + {{-0.333333, 0}, {2.66667, 0}, {-1.33333, 0}, {0.791667, 0.791667}}, + {{-0.333333, 0}, {0.5, -0}, {3, 0}, {-0.166667, -0.166667}}}; + + std::vector internal_loads = compute_internal_loads(solution_set, test_system); + + std::vector test_loads = { + {-0.291667, -0.291667}, {0.0416667, 0.0416667}, {-0.75, -0.75}, {0.458333, 0.458333}, + {0.791667, 0.791667}, {0.166667, 0.166667}, {-0.166667, -0.166667}}; + + compare_vectors(internal_loads, test_loads, 1e-5); + } + + SUBCASE("Four edges, four nodes, two real loads") { + std::vector data = {1, 1, -1}; + std::vector row = {1, 2, 3}; + std::vector col = {0, 0, 0}; + + auto solution_set = generate_input_solution_set(data, row, col, Idx{4}, Idx{1}); + solution_set.extended_rhs = {{1, 0}, {-1, -0}, {-1, -0}, {0, 0}}; + + std::vector> const test_system = {{{3, 0}, {-0.666667, 0}}}; + + std::vector internal_loads = compute_internal_loads(solution_set, test_system); + + std::vector test_loads = {{1, 0}, {-0.333333, -0}, {-0.333333, -0}, {-0.666667, 0}}; + + compare_vectors(internal_loads, test_loads, 1.e-5); + } + } } } // namespace power_grid_model::link_solver