Incremental Y-Bus admittance matrix update

power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp

@@ -159,6 +159,14 @@ template <bool sym> struct MathModelParam {
     ComplexTensorVector<sym> source_param;
 };
 
+template <bool sym> struct MathModelParamIncrement {
+    std::vector<BranchCalcParam<sym>> branch_param;
+    ComplexTensorVector<sym> shunt_param;
+    ComplexTensorVector<sym> source_param;
+    std::vector<Idx> branch_param_to_change;
+    std::vector<Idx> shunt_param_to_change;
+};
+
 template <bool sym> struct PowerFlowInput {
     ComplexVector source;                // Complex u_ref of each source
     ComplexValueVector<sym> s_injection; // Specified injection power of each load_gen

power_grid_model_c/power_grid_model/include/power_grid_model/main_core/update.hpp

@@ -10,6 +10,8 @@
 
 #include "../all_components.hpp"
 
+#include <ranges>
+
 namespace power_grid_model::main_core {
 
 namespace detail {
@@ -82,6 +84,33 @@ inline void update_y_bus(YBus<sym>& y_bus, std::shared_ptr<MathModelParam<sym> c
     y_bus.update_admittance(math_model_param);
 }
 
+template <bool sym>
+inline void update_y_bus_increment(YBus<sym>& y_bus, std::shared_ptr<MathModelParam<sym> const> const& math_model_param,
+                                   bool increment) {
+    auto branch_param_to_change_views =
+        std::views::iota(0, math_model_param->branch_param.size()) | std::views::filter([&math_model_param](Idx i) {
+            return math_model_param->branch_param[i].yff() != ComplexTensor<sym>{0.0} ||
+                   math_model_param->branch_param[i].yft() != ComplexTensor<sym>{0.0} ||
+                   math_model_param->branch_param[i].ytf() != ComplexTensor<sym>{0.0} ||
+                   math_model_param->branch_param[i].ytt() != ComplexTensor<sym>{0.0};
+        });
+    auto shunt_param_to_change_views =
+        std::views::iota(0, math_model_param->shunt_param.size()) | std::views::filter([&math_model_param](Idx i) {
+            return math_model_param->shunt_param[i] != ComplexTensor<sym>{0.0};
+        });
+
+    MathModelParamIncrement<sym> math_model_param_incrmt;
+    math_model_param_incrmt.branch_param = math_model_param->branch_param;
+    math_model_param_incrmt.shunt_param = math_model_param->shunt_param;
+    math_model_param_incrmt.source_param = math_model_param->source_param; // not sure if we actually need this
+    math_model_param_incrmt.branch_param_to_change = branch_param_to_change_views;
+    math_model_param_incrmt.shunt_param_to_change = shunt_param_to_change_views;
+
+    auto math_model_param_incrmt_ptr = std::make_shared<MathModelParamIncrement<sym> const>(math_model_param_incrmt);
+
+    y_bus.update_admittance_increment(math_model_param_incrmt_ptr, !increment);
+}
+
 template <bool sym>
 inline void update_y_bus(MathState& math_state, std::vector<MathModelParam<sym>> const& math_model_params,
                          Idx n_math_solvers) {

power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/y_bus.hpp

@@ -11,6 +11,8 @@
 #include "../sparse_mapping.hpp"
 #include "../three_phase_tensor.hpp"
 
+#include <ranges>
+
 namespace power_grid_model {
 
 // hide implementation in inside namespace
@@ -72,6 +74,7 @@ struct YBusStructure {
     // element of ybus of all components
     std::vector<YBusElement> y_bus_element;
     std::vector<Idx> y_bus_entry_indptr;
+    std::vector<MatrixPos> y_bus_pos_in_entries;
     // sequence entry of bus data
     IdxVector bus_entry;
     // LU csr structure for the y bus with sparse fill-ins
@@ -206,6 +209,7 @@ struct YBusStructure {
                 // all entries in the same position are looped, append indptr
                 // need to be offset by fill-in
                 y_bus_entry_indptr.push_back((it_element - vec_map_element.cbegin()) - fill_in_counter);
+                y_bus_pos_in_entries.push_back(pos);
                 // iterate linear nnz
                 ++nnz_counter;
                 ++nnz_counter_lu;
@@ -350,6 +354,79 @@ template <bool sym> class YBus {
         admittance_ = std::make_shared<ComplexTensorVector<sym> const>(std::move(admittance));
     }
 
+    // void increments_to_entries(std::shared_ptr<IdxVector> affected_entries) {
+    IdxVector increments_to_entries() {
+        // construct affected entries
+        IdxVector affected_entries = {};
+        auto const& y_bus_pos_in_entries = y_bus_struct_->y_bus_pos_in_entries;
+        auto process_params = [&affected_entries, &y_bus_pos_in_entries, this](auto params_to_change, bool is_branch) {
+            for (auto param_to_change : params_to_change) {
+                Idx id_x{0};
+                Idx id_y{0};
+                if (is_branch) {
+                    id_x = math_topology_->branch_bus_idx[param_to_change][0];
+                    id_y = math_topology_->branch_bus_idx[param_to_change][1];
+                } else {
+                    id_x = math_topology_->shunts_per_bus.get_group(param_to_change);
+                    id_y = math_topology_->shunts_per_bus.get_group(param_to_change);
+                }
+                auto it =
+                    std::ranges::find(y_bus_pos_in_entries.begin(), y_bus_pos_in_entries.end(), MatrixPos{id_x, id_y});
+                if (it != y_bus_pos_in_entries.end()) {
+                    affected_entries.push_back(std::distance(y_bus_pos_in_entries.begin(), it));
+                }
+            }
+        };
+
+        process_params(math_model_param_incrmt_->branch_param_to_change, true);
+        process_params(math_model_param_incrmt_->shunt_param_to_change, false);
+        return affected_entries;
+    }
+
+    void update_admittance_increment(std::shared_ptr<MathModelParamIncrement<sym> const> const& math_model_param_incrmt,
+                                     bool is_decrement = false) {
+        const double mode = is_decrement ? -1.0 : 1.0;
+        if (!is_decrement) {
+            // overwrite the old cached parameters in increment mode (update)
+            math_model_param_incrmt_ = math_model_param_incrmt;
+        } else if (!decremented_) {
+            math_model_param_incrmt_ = {}; // reset recorded increment that is already decremented
+            return;
+        }
+
+        // construct admittance data
+        ComplexTensorVector<sym> admittance(nnz());
+        assert(Idx(admittance_->size()) == nnz());
+        std::ranges::copy(*admittance_, admittance.begin());
+        auto const& y_bus_element = y_bus_struct_->y_bus_element;
+        auto const& y_bus_entry_indptr = y_bus_struct_->y_bus_entry_indptr;
+
+        // construct affected entries
+        auto const affected_entries = increments_to_entries();
+
+        // process and update affected entries
+        for (auto const entry : affected_entries) {
+            // start admittance accumulation with zero
+            ComplexTensor<sym> entry_admittance{0.0};
+            // loop over all entries of this position
+            for (Idx element = y_bus_entry_indptr[entry]; element != y_bus_entry_indptr[entry + 1]; ++element) {
+                if (y_bus_element[element].element_type == YBusElementType::shunt) {
+                    // shunt
+                    entry_admittance += math_model_param_incrmt_->shunt_param[y_bus_element[element].idx];
+                } else {
+                    // branch
+                    entry_admittance += math_model_param_incrmt_->branch_param[y_bus_element[element].idx]
+                                            .value[static_cast<Idx>(y_bus_element[element].element_type)];
+                }
+            }
+            // assign
+            admittance[entry] += mode * entry_admittance;
+        }
+        // move to shared ownership
+        admittance_ = std::make_shared<ComplexTensorVector<sym> const>(std::move(admittance));
+        decremented_ = is_decrement;
+    }
+
     ComplexValue<sym> calculate_injection(ComplexValueVector<sym> const& u, Idx bus_number) const {
         Idx const begin = row_indptr()[bus_number];
         Idx const end = row_indptr()[bus_number + 1];
@@ -428,6 +505,10 @@ template <bool sym> class YBus {
 
     // cache the math parameters
     std::shared_ptr<MathModelParam<sym> const> math_model_param_;
+
+    // cache the increment math parameters
+    std::shared_ptr<MathModelParamIncrement<sym> const> math_model_param_incrmt_;
+    bool decremented_ = false;
 };
 
 template class YBus<true>;

tests/cpp_unit_tests/test_y_bus.cpp

@@ -7,6 +7,8 @@
 
 #include <doctest/doctest.h>
 
+#include <ranges>
+
 namespace power_grid_model {
 
 TEST_CASE("Test y bus") {
@@ -32,7 +34,7 @@ TEST_CASE("Test y bus") {
  ^         |             |
  |         |             5
   --- 0 ---              |
-                             X
+                         X
     */
 
     MathModelTopology topo{};
@@ -272,6 +274,181 @@ TEST_CASE("Test fill-in y bus") {
     CHECK(ybus.map_lu_y_bus == map_lu_y_bus);
 }
 
+TEST_CASE("Incremental update y-bus") {
+    /*
+    test Y bus struct
+    [
+            x, x, 0, 0
+            x, x, x, 0
+            0, x, x, x
+            0, 0, x, x
+    ]
+
+     [0]   = Node
+--0--> = Branch (from --id--> to)
+ -X-   = Open switch / not connected
+
+    Topology:
+
+  --- 4 ---               ----- 3 -----
+ |         |             |             |
+ |         v             v             |
+[0]       [1] --- 1 --> [2] --- 2 --> [3]
+ ^         |             |
+ |         |             5
+  --- 0 ---              |
+                         X
+    */
+
+    MathModelTopology topo{};
+    MathModelParam<true> param_sym;
+    topo.phase_shift.resize(4, 0.0);
+    topo.branch_bus_idx = {
+        {1, 0}, // branch 0 from node 1 to 0
+        {1, 2}, // branch 1 from node 1 to 2
+        {2, 3}, // branch 2 from node 2 to 3
+        {3, 2}, // branch 3 from node 3 to 2
+        {0, 1}, // branch 4 from node 0 to 1
+        {2, -1} // branch 5 from node 2 to "not connected"
+    };
+    param_sym.branch_param = {                             //  ff,    ft,    tf,   tt
+                              {1.0i, 2.0i, 3.0i, 4.0i},    // { 1,  0 }
+                              {5.0, 6.0, 7.0, 8.0},        // { 1,  2 }
+                              {9.0i, 10.0i, 11.0i, 12.0i}, // { 2,  3 }
+                              {13.0, 14.0, 15.0, 16.0},    // { 3,  2 }
+                              {17.0, 18.0, 19.0, 20.0},    // { 0,  1 }
+                              {1000i, 0.0, 0.0, 0.0}};     // { 2, -1 }
+    topo.shunts_per_bus = {from_sparse, {0, 1, 1, 1, 2}};  // 4 buses, 2 shunts -> shunt connected to bus 0 and bus 3
+    param_sym.shunt_param = {100.0i, 200.0i};
+
+    // get shared ptr
+    auto topo_ptr = std::make_shared<MathModelTopology const>(topo);
+
+    const ComplexTensorVector<true> admittance_sym = {
+        17.0 + 104.0i,  // 0, 0 -> {1, 0}tt + {0, 1}ff + shunt(0) = 4.0i + 17.0 + 100.0i
+        18.0 + 3.0i,    // 0, 1 -> {0, 1}ft + {1, 0}tf = 18.0 + 3.0i
+        19.0 + 2.0i,    // 1, 0 -> {0, 1}tf + {1, 0}ft = 19.0 + 2.0i
+        25.0 + 1.0i,    // 1, 1 -> {0, 1}tt + {1, 0}ff + {1,2}ff = 20.0 + 1.0i + 5.0
+        6.0,            // 1, 2 -> {1,2}ft = 6.0
+        7.0,            // 2, 1 -> {1,2}tf = 7.0
+        24.0 + 1009.0i, // 2, 2 -> {1,2}tt + {2,3}ff + {3, 2}tt + {2,-1}ff = 8.0 + 9.0i + 16.0 + 1000.0i = 24.0 + 1009i
+        15.0 + 10.0i,   // 2, 3 -> {2,3}ft + {3,2}tf = 10.0i + 15.0
+        14.0 + 11.0i,   // 3, 2 -> {2,3}tf + {3,2}ft = 11.0i + 14.0
+        13.0 + 212.0i   // 3, 3 -> {2,3}tt + {3,2}ff + shunt(1) = 12.0i + 13.0 + 200.0i
+    };
+
+    const ComplexTensorVector<true> admittance_sym_state_1 = {34.0 + 208.0i, 36.0 + 6.0i,  38.0 + 4.0i,    50.0 + 2.0i,
+                                                              12.0,          14.0,         48.0 + 2018.0i, 30.0 + 20.0i,
+                                                              14.0 + 22.0i,  26.0 + 424.0i};
+    topo.branch_bus_idx = {
+        {1, 0}, // branch 0 from node 1 to 0
+        {1, 2}, // branch 1 from node 1 to 2
+        {2, 3}, // branch 2 from node 2 to 3
+        {3, 2}, // branch 3 from node 3 to 2
+        {0, 1}, // branch 4 from node 0 to 1
+        {2, -1} // branch 5 from node 2 to "not connected"
+    };
+    MathModelParam<true> param_sym_update;
+    param_sym_update.branch_param = {                        //   ff,    ft,   tf,   tt
+                                     {1.0i, 0.0, 0.0, 0.0},  // {1,  0}
+                                     {0.0, 1.0, 0.0, 0.0},   // {1,  2}
+                                     {0.0, 0.0, 0.0, 2.0i},  // {2,  3}
+                                     {0.0, 0.0, 2.0, 0.0},   // {3,  2}
+                                     {0.0, 0.0, 0.0, 0.0},   // {0,  1}
+                                     {1.0i, 0.0, 0.0, 0.0}}; // {2, -1}
+    param_sym_update.shunt_param = {1.0i, 0.0i};
+
+    const ComplexTensorVector<true> admittance_sym_2 = {
+        // 17.0 + 104.0i,                                                                        [v]
+        17.0 + 105.0i, // 0, 0 -> += {1, 0}tt + {0, 1}ff + shunt(0) = 0.0 + 0.0 + 1.0i
+        // 18.0 + 3.0i,                                                                          [v]
+        18.0 + 3.0i, // 0, 1 -> += {0, 1}ft + {1, 0}tf = 0.0 + 0.0
+        // 19.0 + 2.0i,                                                                          [v]
+        19.0 + 2.0i, // 1, 0 -> += {0, 1}tf + {1, 0}ft = 0.0 + 0.0
+        // 25.0 + 1.0i,                                                                          [v]
+        25.0 + 2.0i, // 1, 1 -> += {0, 1}tt + {1, 0}ff + {1,2}ff = 0.0 + 1.0i + 0.0
+        // 6.0,                                                                                  [v]
+        7.0, // 1, 2 -> += {1,2}ft = 1.0
+        // 7.0,                                                                                  [v]
+        7.0, // 2, 1 -> += {1,2}tf = 0.0
+        // 24.0 + 1009.0i,                                                                       [v]
+        24.0 + 1010.0i, // 2, 2 -> += {1,2}tt + {2,3}ff + {3, 2}tt + {2,-1}ff = 0.0 + 0.0 + 0.0 + 1.0i
+        // 15.0 + 10.0i,                                                                         [v]
+        17.0 + 10.0i, // 2, 3 -> += {2,3}ft + {3,2}tf = 0.0 + 2.0
+        // 14.0 + 11.0i,                                                                         [v]
+        14.0 + 11.0i, // 3, 2 -> += {2,3}tf + {3,2}ft = 0.0 + 0.0
+        // 13.0 + 212.0i                                                                         [v]
+        13.0 + 214.0i // 3, 3 -> += {2,3}tt + {3,2}ff + shunt(1) = 2.0i + 0.0 + 0.0
+    };
+
+    SUBCASE("Test whole scale update") {
+        YBus<true> ybus{topo_ptr, std::make_shared<MathModelParam<true> const>(param_sym)};
+        CHECK(ybus.admittance().size() == admittance_sym.size());
+        for (size_t i = 0; i < admittance_sym.size(); i++) {
+            CHECK(cabs(ybus.admittance()[i] - admittance_sym[i]) < numerical_tolerance);
+        }
+        ybus.update_admittance(std::make_shared<MathModelParam<true> const>(param_sym));
+
+        CHECK(ybus.admittance().size() == admittance_sym.size());
+        for (size_t i = 0; i < admittance_sym.size(); i++) {
+            CHECK(cabs(ybus.admittance()[i] - admittance_sym[i]) < numerical_tolerance);
+        }
+    }
+
+    SUBCASE("Test progressive update") {
+        YBus<true> ybus{topo_ptr, std::make_shared<MathModelParam<true> const>(param_sym)};
+        CHECK(ybus.admittance().size() == admittance_sym.size());
+        for (size_t i = 0; i < admittance_sym.size(); i++) {
+            CHECK(cabs(ybus.admittance()[i] - admittance_sym[i]) < numerical_tolerance);
+        }
+        auto branch_param_to_change_views =
+            std::views::iota(0, static_cast<int>(param_sym_update.branch_param.size())) |
+            std::views::filter([&param_sym_update](Idx i) {
+                return param_sym_update.branch_param[i].yff() != ComplexTensor<true>{0.0} ||
+                       param_sym_update.branch_param[i].yft() != ComplexTensor<true>{0.0} ||
+                       param_sym_update.branch_param[i].ytf() != ComplexTensor<true>{0.0} ||
+                       param_sym_update.branch_param[i].ytt() != ComplexTensor<true>{0.0};
+            });
+        auto shunt_param_to_change_views = std::views::iota(0, static_cast<int>(param_sym_update.shunt_param.size())) |
+                                           std::views::filter([&param_sym_update](Idx i) {
+                                               return param_sym_update.shunt_param[i] != ComplexTensor<true>{0.0};
+                                           });
+
+        MathModelParamIncrement<true> math_model_param_incrmt;
+        math_model_param_incrmt.branch_param = param_sym_update.branch_param;
+        math_model_param_incrmt.shunt_param = param_sym_update.shunt_param;
+        math_model_param_incrmt.source_param = param_sym_update.source_param; // not sure if we actually need this
+        math_model_param_incrmt.branch_param_to_change = {branch_param_to_change_views.begin(),
+                                                          branch_param_to_change_views.end()};
+        math_model_param_incrmt.shunt_param_to_change = {shunt_param_to_change_views.begin(),
+                                                         shunt_param_to_change_views.end()};
+
+        auto math_model_param_incrmt_ptr =
+            std::make_shared<MathModelParamIncrement<true> const>(math_model_param_incrmt);
+
+        ybus.update_admittance_increment(math_model_param_incrmt_ptr, false);
+        // check increment
+        CHECK(ybus.admittance().size() == admittance_sym_2.size());
+        for (size_t i = 0; i < admittance_sym_2.size(); i++) {
+            CHECK(cabs(ybus.admittance()[i] - admittance_sym_2[i]) < numerical_tolerance);
+        }
+
+        ybus.update_admittance_increment(math_model_param_incrmt_ptr, true);
+        // check decrement
+        CHECK(ybus.admittance().size() == admittance_sym.size());
+        for (size_t i = 0; i < admittance_sym.size(); i++) {
+            CHECK(cabs(ybus.admittance()[i] - admittance_sym[i]) < numerical_tolerance);
+        }
+
+        ybus.update_admittance_increment(math_model_param_incrmt_ptr, true);
+        // check second decrement, which should not change anything
+        CHECK(ybus.admittance().size() == admittance_sym.size());
+        for (size_t i = 0; i < admittance_sym.size(); i++) {
+            CHECK(cabs(ybus.admittance()[i] - admittance_sym[i]) < numerical_tolerance);
+        }
+    }
+}
+
 /*
 TODO:
 - test counting_sort_element()