Incremental Y-Bus admittance matrix update

power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp

@@ -159,6 +159,11 @@ template <bool sym> struct MathModelParam {
     ComplexTensorVector<sym> source_param;
 };
 
+struct MathModelParamIncrement {
+    std::vector<Idx> branch_param_to_change; // indices of changed branch_param
+    std::vector<Idx> shunt_param_to_change;  // indices of changed shunt_param
+};
+
 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/container.hpp

@@ -140,13 +140,20 @@ class Container<RetrievableTypes<GettableTypes...>, StorageableTypes...> {
         return size_[get_cls_pos_v<Gettable, GettableTypes...>];
     }
 
+    // get sequence idx based on idx_2d
+    // E.g. when you know the idx_2d of the derived class but want to know the index of the base class getter
+    template <class Gettable> Idx get_seq(Idx2D idx_2d) const {
+        assert(construction_complete_);
+        std::array<Idx, num_storageable + 1> const& cum_size = cum_size_[get_cls_pos_v<Gettable, GettableTypes...>];
+        return cum_size[idx_2d.group] + idx_2d.pos;
+    }
+
     // get sequence idx based on id
     template <class Gettable> Idx get_seq(ID id) const {
         assert(construction_complete_);
-        std::array<Idx, num_storageable + 1> const& cum_size = cum_size_[get_cls_pos_v<Gettable, GettableTypes...>];
         auto const found = map_.find(id);
         assert(found != map_.end());
-        return cum_size[found->second.group] + found->second.pos;
+        return get_seq<Gettable>(found->second);
     }
 
     // get idx_2d based on sequence

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

@@ -59,22 +59,30 @@ inline std::vector<Idx2D> get_component_sequence(MainModelState<ComponentContain
 // using forward interators
 // different selection based on component type
 // if sequence_idx is given, it will be used to load the object instead of using IDs via hash map.
-template <std::derived_from<Base> Component, class ComponentContainer, std::forward_iterator ForwardIterator>
+template <std::derived_from<Base> Component, class ComponentContainer, std::forward_iterator ForwardIterator,
+          std::output_iterator<Idx2D> OutputIterator>
     requires model_component_state<MainModelState, ComponentContainer, Component>
 inline UpdateChange update_component(MainModelState<ComponentContainer>& state, ForwardIterator begin,
-                                     ForwardIterator end, std::vector<Idx2D> const& sequence_idx = {}) {
+                                     ForwardIterator end, OutputIterator changed_it,
+                                     std::vector<Idx2D> const& sequence_idx = {}) {
     using UpdateType = typename Component::UpdateType;
 
-    UpdateChange changed;
+    UpdateChange state_changed;
 
     detail::iterate_component_sequence<Component>(
-        [&changed, &state](UpdateType const& update_data, Idx2D const& sequence_single) {
+        [&state_changed, &changed_it, &state](UpdateType const& update_data, Idx2D const& sequence_single) {
             auto& comp = state.components.template get_item<Component>(sequence_single);
-            changed = changed || comp.update(update_data);
+            auto const comp_changed = comp.update(update_data);
+
+            state_changed = state_changed || comp_changed;
+
+            if (comp_changed.param || comp_changed.topo) {
+                *changed_it++ = sequence_single;
+            }
         },
         begin, end, sequence_idx);
 
-    return changed;
+    return state_changed;
 }
 
 // template to get the inverse update for components
@@ -97,23 +105,42 @@ inline void update_inverse(MainModelState<ComponentContainer> const& state, Forw
 }
 
 template <bool sym>
-inline void update_y_bus(YBus<sym>& y_bus, std::shared_ptr<MathModelParam<sym> const> const& math_model_param) {
-    y_bus.update_admittance(math_model_param);
+inline void update_y_bus(MathState& math_state, std::vector<MathModelParam<sym>> const& math_model_params) {
+    auto& y_bus_vec = [&math_state]() -> auto& {
+        if constexpr (sym) {
+            return math_state.y_bus_vec_sym;
+        } else {
+            return math_state.y_bus_vec_asym;
+        }
+    }
+    ();
+
+    assert(y_bus_vec.size() == math_model_params.size());
+
+    for (Idx i = 0; i != static_cast<Idx>(y_bus_vec.size()); ++i) {
+        y_bus_vec[i].update_admittance(std::make_shared<MathModelParam<sym> const>(std::move(math_model_params[i])));
+    }
 }
 
 template <bool sym>
 inline void update_y_bus(MathState& math_state, std::vector<MathModelParam<sym>> const& math_model_params,
-                         Idx n_math_solvers) {
-    for (Idx i = 0; i != n_math_solvers; ++i) {
+                         std::vector<MathModelParamIncrement> const& math_model_param_increments) {
+    auto& y_bus_vec = [&math_state]() -> auto& {
         if constexpr (sym) {
-            update_y_bus(math_state.y_bus_vec_sym[i],
-                         std::make_shared<MathModelParam<sym> const>(std::move(math_model_params[i])));
-
+            return math_state.y_bus_vec_sym;
         } else {
-            update_y_bus(math_state.y_bus_vec_asym[i],
-                         std::make_shared<MathModelParam<sym> const>(std::move(math_model_params[i])));
+            return math_state.y_bus_vec_asym;
         }
     }
+    ();
+
+    assert(y_bus_vec.size() == math_model_params.size());
+
+    for (Idx i = 0; i != static_cast<Idx>(y_bus_vec.size()); ++i) {
+        y_bus_vec[i].update_admittance_increment(
+            std::make_shared<MathModelParam<sym> const>(std::move(math_model_params[i])),
+            math_model_param_increments[i]);
+    }
 }
 
 } // namespace power_grid_model::main_core

power_grid_model_c/power_grid_model/include/power_grid_model/main_model.hpp

@@ -152,17 +152,18 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
     // if sequence_idx is given, it will be used to load the object instead of using IDs via hash map.
     template <class CompType, class CacheType, std::forward_iterator ForwardIterator>
     void update_component(ForwardIterator begin, ForwardIterator end, std::vector<Idx2D> const& sequence_idx) {
+        constexpr auto comp_index = AllComponents::template index_of<CompType>();
+
         assert(construction_complete_);
         assert(static_cast<ptrdiff_t>(sequence_idx.size()) == std::distance(begin, end));
 
         if constexpr (CacheType::value) {
-            constexpr auto comp_index = AllComponents::template index_of<CompType>();
-
             main_core::update_inverse<CompType>(
                 state_, begin, end, std::back_inserter(std::get<comp_index>(cached_inverse_update_)), sequence_idx);
         }
 
-        UpdateChange const changed = main_core::update_component<CompType>(state_, begin, end, sequence_idx);
+        UpdateChange const changed = main_core::update_component<CompType>(
+            state_, begin, end, std::back_inserter(std::get<comp_index>(parameter_changed_components_)), sequence_idx);
 
         // update, get changed variable
         update_state(changed);
@@ -748,7 +749,7 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
     }
 
     template <typename Component, math_output_type MathOutputType, std::forward_iterator ResIt>
-    ResIt output_result(std::vector<MathOutputType> const& math_output, ResIt res_it) {
+    ResIt output_result(std::vector<MathOutputType> const& math_output, ResIt res_it) const {
         assert(construction_complete_);
         return main_core::output_result<Component, ComponentContainer>(state_, math_output, res_it);
     }
@@ -783,7 +784,7 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
         }
     }
 
-    CalculationInfo calculation_info() { return calculation_info_; }
+    CalculationInfo calculation_info() const { return calculation_info_; }
 
   private:
     CalculationInfo calculation_info_; // needs to be first due to padding override
@@ -797,9 +798,12 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
     bool is_topology_up_to_date_{false};
     bool is_sym_parameter_up_to_date_{false};
     bool is_asym_parameter_up_to_date_{false};
+    bool is_accumulated_component_updated_{true};
+    bool last_updated_calculation_symmetry_mode_{false};
 
     OwnedUpdateDataset cached_inverse_update_{};
     UpdateChange cached_state_changes_{};
+    std::array<std::vector<Idx2D>, n_types> parameter_changed_components_{};
 #ifndef NDEBUG
     // construction_complete is used for debug assertions only
     bool construction_complete_{false};
@@ -877,7 +881,7 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
             math_param[i].source_param.resize(state_.math_topology[i]->n_source());
         }
         // loop all branch
-        for (Idx i = 0; i != (Idx)state_.comp_topo->branch_node_idx.size(); ++i) {
+        for (Idx i = 0; i != static_cast<Idx>(state_.comp_topo->branch_node_idx.size()); ++i) {
             Idx2D const math_idx = state_.topo_comp_coup->branch[i];
             if (math_idx.group == -1) {
                 continue;
@@ -887,7 +891,7 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
                 state_.components.template get_item_by_seq<Branch>(i).template calc_param<sym>();
         }
         // loop all branch3
-        for (Idx i = 0; i != (Idx)state_.comp_topo->branch3_node_idx.size(); ++i) {
+        for (Idx i = 0; i != static_cast<Idx>(state_.comp_topo->branch3_node_idx.size()); ++i) {
             Idx2DBranch3 const math_idx = state_.topo_comp_coup->branch3[i];
             if (math_idx.group == -1) {
                 continue;
@@ -900,7 +904,7 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
             }
         }
         // loop all shunt
-        for (Idx i = 0; i != (Idx)state_.comp_topo->shunt_node_idx.size(); ++i) {
+        for (Idx i = 0; i != static_cast<Idx>(state_.comp_topo->shunt_node_idx.size()); ++i) {
             Idx2D const math_idx = state_.topo_comp_coup->shunt[i];
             if (math_idx.group == -1) {
                 continue;
@@ -910,7 +914,7 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
                 state_.components.template get_item_by_seq<Shunt>(i).template calc_param<sym>();
         }
         // loop all source
-        for (Idx i = 0; i != (Idx)state_.comp_topo->source_node_idx.size(); ++i) {
+        for (Idx i = 0; i != static_cast<Idx>(state_.comp_topo->source_node_idx.size()); ++i) {
             Idx2D const math_idx = state_.topo_comp_coup->source[i];
             if (math_idx.group == -1) {
                 continue;
@@ -921,6 +925,56 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
         }
         return math_param;
     }
+    template <bool sym> std::vector<MathModelParamIncrement> get_math_param_increment() {
+        using AddToIncrement = void (*)(std::vector<MathModelParamIncrement>&, MainModelState const&, Idx2D const&);
+
+        static constexpr std::array<AddToIncrement, n_types> add_to_increments{
+            [](std::vector<MathModelParamIncrement>& increments, MainModelState const& state,
+               Idx2D const& changed_component_idx) {
+                if constexpr (std::derived_from<ComponentType, Branch>) {
+                    Idx2D const math_idx =
+                        state.topo_comp_coup->branch[state.components.template get_seq<Branch>(changed_component_idx)];
+                    if (math_idx.group == -1) {
+                        return;
+                    }
+                    // assign parameters
+                    increments[math_idx.group].branch_param_to_change.push_back(math_idx.pos);
+                } else if constexpr (std::derived_from<ComponentType, Branch3>) {
+                    Idx2DBranch3 const math_idx =
+                        state.topo_comp_coup
+                            ->branch3[state.components.template get_seq<Branch3>(changed_component_idx)];
+                    if (math_idx.group == -1) {
+                        return;
+                    }
+                    // assign parameters, branch3 param consists of three branch parameters
+                    // auto const branch3_param =
+                    //   state.components.template get_item<Branch3>(changed_component_idx).template calc_param<sym>();
+                    for (size_t branch2 = 0; branch2 < 3; ++branch2) {
+                        increments[math_idx.group].branch_param_to_change.push_back(math_idx.pos[branch2]);
+                    }
+                } else if constexpr (std::same_as<ComponentType, Shunt>) {
+                    Idx2D const math_idx =
+                        state.topo_comp_coup->shunt[state.components.template get_seq<Shunt>(changed_component_idx)];
+                    if (math_idx.group == -1) {
+                        return;
+                    }
+                    // assign parameters
+                    increments[math_idx.group].shunt_param_to_change.push_back(math_idx.pos);
+                }
+            }...};
+
+        std::vector<MathModelParamIncrement> math_param_increment(n_math_solvers_);
+
+        for (size_t i = 0; i < n_types; ++i) {
+            auto const& changed_type_components = parameter_changed_components_[i];
+            auto const& add_type_to_increment = add_to_increments[i];
+            for (auto const& changed_component : changed_type_components) {
+                add_type_to_increment(math_param_increment, state_, changed_component);
+            }
+        }
+
+        return math_param_increment;
+    }
 
     static constexpr auto include_all = [](Idx) { return true; };
 
@@ -1181,6 +1235,12 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
             y_bus_vec.reserve(n_math_solvers_);
             auto math_params = get_math_param<sym>();
 
+            // Check the branch and shunt indices
+            constexpr auto branch_param_in_seq_map =
+                std::array{AllComponents::template index_of<Line>(), AllComponents::template index_of<Link>(),
+                           AllComponents::template index_of<Transformer>()};
+            constexpr auto shunt_param_in_seq_map = std::array{AllComponents::template index_of<Shunt>()};
+
             for (Idx i = 0; i != n_math_solvers_; ++i) {
                 // construct from existing Y_bus structure if possible
                 if (other_y_bus_exist) {
@@ -1191,6 +1251,11 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
                     y_bus_vec.emplace_back(state_.math_topology[i],
                                            std::make_shared<MathModelParam<sym> const>(std::move(math_params[i])));
                 }
+
+                y_bus_vec.back().set_branch_param_idx(
+                    IdxVector{branch_param_in_seq_map.begin(), branch_param_in_seq_map.end()});
+                y_bus_vec.back().set_shunt_param_idx(
+                    IdxVector{shunt_param_in_seq_map.begin(), shunt_param_in_seq_map.end()});
             }
         }
     }
@@ -1202,23 +1267,28 @@ class MainModelImpl<ExtraRetrievableTypes<ExtraRetrievableType...>, ComponentLis
             rebuild_topology();
         }
         prepare_y_bus<sym>();
-        // if solvers do not exist, build them
-        if (n_math_solvers_ != (Idx)solvers.size()) {
+
+        if (n_math_solvers_ != static_cast<Idx>(solvers.size())) {
             assert(solvers.empty());
+            assert(n_math_solvers_ == static_cast<Idx>(state_.math_topology.size()));
+
+            solvers.clear();
             solvers.reserve(n_math_solvers_);
-            // loop to build
-            for (Idx i = 0; i != n_math_solvers_; ++i) {
-                solvers.emplace_back(state_.math_topology[i]);
-            }
-        }
-        // if parameters are not up to date, update them
-        else if (!is_parameter_up_to_date<sym>()) {
-            // get param, will be consumed
+            std::ranges::transform(state_.math_topology, std::back_inserter(solvers),
+                                   [](auto math_topo) { return MathSolver<sym>{std::move(math_topo)}; });
+        } else if (!is_parameter_up_to_date<sym>()) {
             std::vector<MathModelParam<sym>> const math_params = get_math_param<sym>();
-            main_core::update_y_bus(math_state_, math_params, n_math_solvers_);
+            std::vector<MathModelParamIncrement> const math_param_increments = get_math_param_increment<sym>();
+            if (last_updated_calculation_symmetry_mode_ == sym) {
+                main_core::update_y_bus(math_state_, math_params, math_param_increments);
+            } else {
+                main_core::update_y_bus(math_state_, math_params);
+            }
         }
         // else do nothing, set everything up to date
         is_parameter_up_to_date<sym>() = true;
+        std::ranges::for_each(parameter_changed_components_, [](auto& comps) { comps.clear(); });
+        last_updated_calculation_symmetry_mode_ = sym;
     }
 };