diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp index 1acdc400f..b81c3efb6 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/calculation_parameters.hpp @@ -159,6 +159,11 @@ template struct MathModelParam { ComplexTensorVector source_param; }; +struct MathModelParamIncrement { + std::vector branch_param_to_change; // indices of changed branch_param + std::vector shunt_param_to_change; // indices of changed shunt_param +}; + template struct PowerFlowInput { ComplexVector source; // Complex u_ref of each source ComplexValueVector s_injection; // Specified injection power of each load_gen diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/container.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/container.hpp index 781be626f..63d5dca2f 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/container.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/container.hpp @@ -140,13 +140,20 @@ class Container, StorageableTypes...> { return size_[get_cls_pos_v]; } + // 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 Idx get_seq(Idx2D idx_2d) const { + assert(construction_complete_); + std::array const& cum_size = cum_size_[get_cls_pos_v]; + return cum_size[idx_2d.group] + idx_2d.pos; + } + // get sequence idx based on id template Idx get_seq(ID id) const { assert(construction_complete_); - std::array const& cum_size = cum_size_[get_cls_pos_v]; auto const found = map_.find(id); assert(found != map_.end()); - return cum_size[found->second.group] + found->second.pos; + return get_seq(found->second); } // get idx_2d based on sequence diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/main_core/update.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/main_core/update.hpp index ce012aa50..85792bfac 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/main_core/update.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/main_core/update.hpp @@ -59,22 +59,30 @@ inline std::vector get_component_sequence(MainModelState Component, class ComponentContainer, std::forward_iterator ForwardIterator> +template Component, class ComponentContainer, std::forward_iterator ForwardIterator, + std::output_iterator OutputIterator> requires model_component_state inline UpdateChange update_component(MainModelState& state, ForwardIterator begin, - ForwardIterator end, std::vector const& sequence_idx = {}) { + ForwardIterator end, OutputIterator changed_it, + std::vector const& sequence_idx = {}) { using UpdateType = typename Component::UpdateType; - UpdateChange changed; + UpdateChange state_changed; detail::iterate_component_sequence( - [&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(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 const& state, Forw } template -inline void update_y_bus(YBus& y_bus, std::shared_ptr const> const& math_model_param) { - y_bus.update_admittance(math_model_param); +inline void update_y_bus(MathState& math_state, std::vector> 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(y_bus_vec.size()); ++i) { + y_bus_vec[i].update_admittance(std::make_shared const>(std::move(math_model_params[i]))); + } } template inline void update_y_bus(MathState& math_state, std::vector> const& math_model_params, - Idx n_math_solvers) { - for (Idx i = 0; i != n_math_solvers; ++i) { + std::vector 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 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 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(y_bus_vec.size()); ++i) { + y_bus_vec[i].update_admittance_increment( + std::make_shared const>(std::move(math_model_params[i])), + math_model_param_increments[i]); + } } } // namespace power_grid_model::main_core diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/main_model.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/main_model.hpp index 348f99b21..095fe0aac 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/main_model.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/main_model.hpp @@ -152,17 +152,18 @@ class MainModelImpl, ComponentLis // if sequence_idx is given, it will be used to load the object instead of using IDs via hash map. template void update_component(ForwardIterator begin, ForwardIterator end, std::vector const& sequence_idx) { + constexpr auto comp_index = AllComponents::template index_of(); + assert(construction_complete_); assert(static_cast(sequence_idx.size()) == std::distance(begin, end)); if constexpr (CacheType::value) { - constexpr auto comp_index = AllComponents::template index_of(); - main_core::update_inverse( state_, begin, end, std::back_inserter(std::get(cached_inverse_update_)), sequence_idx); } - UpdateChange const changed = main_core::update_component(state_, begin, end, sequence_idx); + UpdateChange const changed = main_core::update_component( + state_, begin, end, std::back_inserter(std::get(parameter_changed_components_)), sequence_idx); // update, get changed variable update_state(changed); @@ -748,7 +749,7 @@ class MainModelImpl, ComponentLis } template - ResIt output_result(std::vector const& math_output, ResIt res_it) { + ResIt output_result(std::vector const& math_output, ResIt res_it) const { assert(construction_complete_); return main_core::output_result(state_, math_output, res_it); } @@ -783,7 +784,7 @@ class MainModelImpl, 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, 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, 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, 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(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, ComponentLis state_.components.template get_item_by_seq(i).template calc_param(); } // loop all branch3 - for (Idx i = 0; i != (Idx)state_.comp_topo->branch3_node_idx.size(); ++i) { + for (Idx i = 0; i != static_cast(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, 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(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, ComponentLis state_.components.template get_item_by_seq(i).template calc_param(); } // loop all source - for (Idx i = 0; i != (Idx)state_.comp_topo->source_node_idx.size(); ++i) { + for (Idx i = 0; i != static_cast(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, ComponentLis } return math_param; } + template std::vector get_math_param_increment() { + using AddToIncrement = void (*)(std::vector&, MainModelState const&, Idx2D const&); + + static constexpr std::array add_to_increments{ + [](std::vector& increments, MainModelState const& state, + Idx2D const& changed_component_idx) { + if constexpr (std::derived_from) { + Idx2D const math_idx = + state.topo_comp_coup->branch[state.components.template get_seq(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) { + Idx2DBranch3 const math_idx = + state.topo_comp_coup + ->branch3[state.components.template get_seq(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(changed_component_idx).template calc_param(); + 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) { + Idx2D const math_idx = + state.topo_comp_coup->shunt[state.components.template get_seq(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 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, ComponentLis y_bus_vec.reserve(n_math_solvers_); auto math_params = get_math_param(); + // Check the branch and shunt indices + constexpr auto branch_param_in_seq_map = + std::array{AllComponents::template index_of(), AllComponents::template index_of(), + AllComponents::template index_of()}; + constexpr auto shunt_param_in_seq_map = std::array{AllComponents::template index_of()}; + 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, ComponentLis y_bus_vec.emplace_back(state_.math_topology[i], std::make_shared 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,33 @@ class MainModelImpl, ComponentLis rebuild_topology(); } prepare_y_bus(); - // if solvers do not exist, build them - if (n_math_solvers_ != (Idx)solvers.size()) { + + if (n_math_solvers_ != static_cast(solvers.size())) { assert(solvers.empty()); + assert(n_math_solvers_ == static_cast(state_.math_topology.size())); + assert(n_math_solvers_ == static_cast(get_y_bus().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]); + std::ranges::transform(state_.math_topology, std::back_inserter(solvers), + [](auto math_topo) { return MathSolver{std::move(math_topo)}; }); + for (Idx idx = 0; idx < n_math_solvers_; ++idx) { + get_y_bus()[idx].register_parameters_changed_callback( + [solver = std::ref(solvers[idx])](bool changed) { solver.get().parameters_changed(changed); }); } - } - // if parameters are not up to date, update them - else if (!is_parameter_up_to_date()) { - // get param, will be consumed + } else if (!is_parameter_up_to_date()) { std::vector> const math_params = get_math_param(); - main_core::update_y_bus(math_state_, math_params, n_math_solvers_); + std::vector const math_param_increments = get_math_param_increment(); + 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() = true; + std::ranges::for_each(parameter_changed_components_, [](auto& comps) { comps.clear(); }); + last_updated_calculation_symmetry_mode_ = sym; } }; diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_current_pf_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_current_pf_solver.hpp index 4d8939729..32f29159f 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_current_pf_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/iterative_current_pf_solver.hpp @@ -80,7 +80,6 @@ template class IterativeCurrentPFSolver : public IterativePFSolver const& y_bus, std::shared_ptr const& topo_ptr) : IterativePFSolver{y_bus, topo_ptr}, rhs_u_(y_bus.size()), - y_data_ptr_(nullptr), sparse_solver_{y_bus.shared_indptr_lu(), y_bus.shared_indices_lu(), y_bus.shared_diag_lu()} {} // Add source admittance to Y bus and set variable for prepared y bus to true @@ -89,7 +88,7 @@ template class IterativeCurrentPFSolver : public IterativePFSolver mat_data(y_bus.nnz_lu()); detail::copy_y_bus(y_bus, mat_data); @@ -106,9 +105,8 @@ template class IterativeCurrentPFSolver : public IterativePFSolver const>(std::move(mat_data)); perm_ = std::make_shared(std::move(perm)); - // cache pointer - y_data_ptr_ = &y_bus.admittance(); } + parameters_changed_ = false; } // Prepare matrix calculates injected current ie. RHS of solver for each iteration. @@ -146,13 +144,15 @@ template class IterativeCurrentPFSolver : public IterativePFSolver rhs_u_; std::shared_ptr const> mat_data_; - ComplexTensorVector const* y_data_ptr_; // sparse solver SparseLUSolver, ComplexValue, ComplexValue> sparse_solver_; std::shared_ptr perm_; + bool parameters_changed_ = true; void add_loads(boost::iterator_range const& load_gens, Idx bus_number, PowerFlowInput const& input, std::vector const& load_gen_type, ComplexValueVector const& u) { diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/math_solver.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/math_solver.hpp index aa6491667..4abbb46fb 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/math_solver.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/math_solver.hpp @@ -98,6 +98,12 @@ template class MathSolver { iterative_linear_se_solver_.reset(); } + void parameters_changed(bool changed) { + if (iterative_current_pf_solver_.has_value()) { + iterative_current_pf_solver_->parameters_changed(changed); + } + } + private: std::shared_ptr topo_ptr_; bool all_const_y_; // if all the load_gen is const element_admittance (impedance) type diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/y_bus.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/y_bus.hpp index c73d81c7f..e3ce35183 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/y_bus.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/math_solver/y_bus.hpp @@ -10,6 +10,9 @@ #include "../power_grid_model.hpp" #include "../three_phase_tensor.hpp" +#include +#include + namespace power_grid_model { // hide implementation in inside namespace @@ -71,6 +74,7 @@ struct YBusStructure { // element of ybus of all components std::vector y_bus_element; std::vector y_bus_entry_indptr; + std::vector y_bus_pos_in_entries; // sequence entry of bus data IdxVector bus_entry; // LU csr structure for the y bus with sparse fill-ins @@ -205,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; @@ -272,6 +277,8 @@ struct YBusStructure { // See also "Node Admittance Matrix" in "State Estimation Alliander" template class YBus { public: + using ParamChangedCallback = std::function; + YBus(std::shared_ptr const& topo_ptr, std::shared_ptr const> const& param, std::shared_ptr const& y_bus_struct = {}) @@ -300,7 +307,7 @@ template class YBus { MathModelTopology const& math_topology() const { return *math_topology_; } MathModelParam const& math_model_param() const { return *math_model_param_; } - ComplexTensorVector const& admittance() const { return *admittance_; } + ComplexTensorVector const& admittance() const { return admittance_; } IdxVector const& bus_entry() const { return y_bus_struct_->bus_entry; } IdxVector const& lu_diag() const { return y_bus_struct_->diag_lu; } IdxVector const& map_lu_y_bus() const { return y_bus_struct_->map_lu_y_bus; } @@ -318,11 +325,17 @@ template class YBus { constexpr auto& get_y_bus_structure() const { return y_bus_struct_; } + void set_branch_param_idx(IdxVector branch_param_idx) { branch_param_idx_ = std::move(branch_param_idx); } + void set_shunt_param_idx(IdxVector shunt_param_idx) { shunt_param_idx_ = std::move(shunt_param_idx); } + IdxVector const& get_branch_param_idx() const { return branch_param_idx_; } + IdxVector const& get_shunt_param_idx() const { return shunt_param_idx_; } + void update_admittance(std::shared_ptr const> const& math_model_param) { // overwrite the old cached parameters math_model_param_ = math_model_param; // construct admittance data - ComplexTensorVector admittance(nnz()); + admittance_ = ComplexTensorVector(nnz()); + auto const& branch_param = math_model_param_->branch_param; auto const& shunt_param = math_model_param_->shunt_param; auto const& y_bus_element = y_bus_struct_->y_bus_element; @@ -331,22 +344,85 @@ template class YBus { for (Idx entry = 0; entry != nnz(); ++entry) { // start admittance accumulation with zero ComplexTensor entry_admittance{0.0}; + IdxVector entry_param_shunt; + IdxVector entry_param_branch; + // loop over all entries of this position + for (Idx element = y_bus_entry_indptr[entry]; element != y_bus_entry_indptr[entry + 1]; ++element) { + auto param_idx = y_bus_element[element].idx; + if (y_bus_element[element].element_type == YBusElementType::shunt) { + entry_admittance += shunt_param[param_idx]; + entry_param_shunt.push_back(param_idx); + } else { + entry_admittance += + branch_param[param_idx].value[static_cast(y_bus_element[element].element_type)]; + entry_param_branch.push_back(param_idx); + } + } + // assign + admittance_[entry] = entry_admittance; + map_admittance_param_branch_.push_back(entry_param_branch); + map_admittance_param_shunt_.push_back(entry_param_shunt); + } + + parameters_changed(true); + } + + IdxVector increments_to_entries(auto const& math_model_param_incrmt) { + // construct affected entries + IdxVector affected_entries; + + auto query_params_in_map = [&affected_entries](auto const& params_to_change, auto const& mapping) { + for (size_t i = 0; i < mapping.size(); ++i) { + if (std::ranges::any_of(mapping[i], [&](Idx val) { + return std::ranges::find(params_to_change, val) != params_to_change.end(); + })) { + affected_entries.push_back(static_cast(i)); + } + } + }; + + query_params_in_map(math_model_param_incrmt.branch_param_to_change, map_admittance_param_branch_); + query_params_in_map(math_model_param_incrmt.shunt_param_to_change, map_admittance_param_shunt_); + return affected_entries; + } + + /** + * @brief Updates the admittance of the y_bus according to what's changed in math_model. + * + * @param math_model_param Shared pointer to the constant math_model parameters. + * @param math_model_param_incrmt Shared pointer to the constant mathematical model parameters . + */ + void update_admittance_increment(std::shared_ptr const> const& math_model_param, + MathModelParamIncrement const& math_model_param_incrmt) { + // swap the old cached parameters + math_model_param_ = math_model_param; + + auto const& y_bus_element = y_bus_struct_->y_bus_element; + auto const& y_bus_entry_indptr = y_bus_struct_->y_bus_entry_indptr; + auto const& math_param_shunt = math_model_param_->shunt_param; + auto const& math_param_branch = math_model_param_->branch_param; + + // process and update affected entries + for (auto const affected_entries = increments_to_entries(math_model_param_incrmt); + auto const entry : affected_entries) { + // start admittance accumulation with zero + ComplexTensor 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 += shunt_param[y_bus_element[element].idx]; + entry_admittance += math_param_shunt[y_bus_element[element].idx]; } else { // branch - entry_admittance += branch_param[y_bus_element[element].idx] + entry_admittance += math_param_branch[y_bus_element[element].idx] .value[static_cast(y_bus_element[element].element_type)]; } } // assign - admittance[entry] = entry_admittance; + admittance_[entry] = entry_admittance; } - // move to shared ownership - admittance_ = std::make_shared const>(std::move(admittance)); + + parameters_changed(true); } ComplexValue calculate_injection(ComplexValueVector const& u, Idx bus_number) const { @@ -415,18 +491,55 @@ template class YBus { return shunt_flow; } + /// @brief register a new callback to signal a parameter change + /// @param callback the callback to register + /// @return the unique key referencing this callback (used for unregistering) + uint64_t register_parameters_changed_callback(ParamChangedCallback callback) { + static uint64_t num_added = 0; + + auto const new_key = num_added; + ++num_added; + + assert(!parameters_changed_callbacks_.contains(new_key)); + parameters_changed_callbacks_.emplace_hint(parameters_changed_callbacks_.cend(), new_key, std::move(callback)); + return new_key; + } + + /// @brief unregister a callback to signal a parameter change + /// @param key the unique key referencing the callback (returned by register_parameters_changed_callback) + void unregister_parameters_changed_callback(uint64_t key) { + assert(parameters_changed_callbacks_.contains(key)); + parameters_changed_callbacks_.erase(key); + } + private: // csr structure std::shared_ptr y_bus_struct_; // admittance - std::shared_ptr const> admittance_; + ComplexTensorVector admittance_; // cache math topology std::shared_ptr math_topology_; // cache the math parameters std::shared_ptr const> math_model_param_; + + // cache the branch and shunt parameters in sequence_idx_map + IdxVector branch_param_idx_{}; + IdxVector shunt_param_idx_{}; + + // map index between admittance entries and parameter entries + std::vector map_admittance_param_branch_{}; + std::vector map_admittance_param_shunt_{}; + + std::unordered_map parameters_changed_callbacks_; + + void parameters_changed(bool param_changed) const { + std::ranges::for_each(parameters_changed_callbacks_, [param_changed](auto const& key_and_callback) { + key_and_callback.second(param_changed); + }); + } }; template class YBus; @@ -438,4 +551,4 @@ using math_solver::YBus; } // namespace power_grid_model -#endif \ No newline at end of file +#endif diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/three_phase_tensor.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/three_phase_tensor.hpp index e36476ad4..d7aeca0c9 100644 --- a/power_grid_model_c/power_grid_model/include/power_grid_model/three_phase_tensor.hpp +++ b/power_grid_model_c/power_grid_model/include/power_grid_model/three_phase_tensor.hpp @@ -100,7 +100,6 @@ template using RealDiagonalTensor = std::conditional_t>; template using ComplexDiagonalTensor = std::conditional_t>; - template using RealValue = std::conditional_t>; template using ComplexValue = std::conditional_t>; diff --git a/tests/cpp_unit_tests/test_container.cpp b/tests/cpp_unit_tests/test_container.cpp index 087ef480c..34830bfc7 100644 --- a/tests/cpp_unit_tests/test_container.cpp +++ b/tests/cpp_unit_tests/test_container.cpp @@ -26,11 +26,9 @@ struct C2 : C { TEST_CASE("Test component container") { using CompContainer = Container; - CompContainer container; - container.reserve(3); - container.reserve(2); - container.reserve(1); using CompContainer2 = Container, C1, C2>; + + CompContainer container; CompContainer2 container2; container.emplace(1, 5); @@ -112,6 +110,18 @@ TEST_CASE("Test component container") { CHECK(const_container.size() == 1); } + SUBCASE("Test get sequence based on idx_2d") { + CHECK(const_container.get_seq(Idx2D{0, 0}) == 0); + CHECK(const_container.get_seq(Idx2D{0, 1}) == 1); + CHECK(const_container.get_seq(Idx2D{0, 2}) == 2); + CHECK(const_container.get_seq(Idx2D{1, 0}) == 3); + CHECK(const_container.get_seq(Idx2D{1, 1}) == 4); + CHECK(const_container.get_seq(Idx2D{2, 0}) == 5); + CHECK(const_container.get_seq(Idx2D{1, 0}) == 0); + CHECK(const_container.get_seq(Idx2D{1, 1}) == 1); + CHECK(const_container.get_seq(Idx2D{2, 0}) == 0); + } + SUBCASE("Test get sequence based on id") { CHECK(const_container.get_seq(1) == 0); CHECK(const_container.get_seq(11) == 1); diff --git a/tests/cpp_unit_tests/test_main_model.cpp b/tests/cpp_unit_tests/test_main_model.cpp index b6fc2b0c3..8f5253f05 100644 --- a/tests/cpp_unit_tests/test_main_model.cpp +++ b/tests/cpp_unit_tests/test_main_model.cpp @@ -149,6 +149,7 @@ struct State { std::vector sym_load_update{{7, 1, 1.0e6, nan}}; std::vector asym_load_update{{8, 0, RealValue{nan}, RealValue{nan}}}; std::vector shunt_update{{9, 0, nan, 0.02, nan, 0.02}}; + std::vector shunt_update_2{{6, 0, nan, 0.01, nan, 0.01}}; // used for test case alternate compute mode std::vector source_update{{10, 1, test::u1, nan}}; std::vector link_update{{5, 1, 0}}; std::vector fault_update{{30, 1, FaultType::three_phase, FaultPhase::abc, 1, nan, nan}}; @@ -960,6 +961,82 @@ TEST_CASE_TEMPLATE("Test main model - restore components", settings, regular_upd } } +TEST_CASE_TEMPLATE("Test main model - updates w/ alternating compute mode", settings, regular_update, cached_update) { + constexpr auto check_sym = [](MainModel const& model_, std::vector> const& math_output_) { + State state_; + model_.output_result(math_output_, state_.sym_node.begin()); + model_.output_result(math_output_, state_.sym_branch.begin()); + model_.output_result(math_output_, state_.sym_appliance.begin()); + + CHECK(state_.sym_node[0].u_pu == doctest::Approx(1.05)); + CHECK(state_.sym_node[1].u_pu == doctest::Approx(test::u1)); + CHECK(state_.sym_node[2].u_pu == doctest::Approx(test::u1)); + CHECK(state_.sym_branch[0].i_from == doctest::Approx(test::i)); + CHECK(state_.sym_appliance[0].i == doctest::Approx(test::i)); + CHECK(state_.sym_appliance[1].i == doctest::Approx(0.0)); + CHECK(state_.sym_appliance[2].i == doctest::Approx(test::i_load * 2 + test::i_shunt)); + CHECK(state_.sym_appliance[3].i == doctest::Approx(0.0)); + CHECK(state_.sym_appliance[4].i == doctest::Approx(0.0)); + }; + constexpr auto check_asym = [](MainModel const& model_, std::vector> const& math_output_) { + State state_; + model_.output_result(math_output_, state_.asym_node.begin()); + model_.output_result(math_output_, state_.asym_branch.begin()); + model_.output_result(math_output_, state_.asym_appliance.begin()); + CHECK(state_.asym_node[0].u_pu(0) == doctest::Approx(1.05)); + CHECK(state_.asym_node[1].u_pu(1) == doctest::Approx(test::u1)); + CHECK(state_.asym_node[2].u_pu(2) == doctest::Approx(test::u1)); + CHECK(state_.asym_branch[0].i_from(0) == doctest::Approx(test::i)); + CHECK(state_.asym_appliance[0].i(1) == doctest::Approx(test::i)); + CHECK(state_.asym_appliance[1].i(2) == doctest::Approx(0.0)); + CHECK(state_.asym_appliance[2].i(0) == doctest::Approx(test::i_load * 2 + test::i_shunt)); + CHECK(state_.asym_appliance[3].i(1) == doctest::Approx(0.0)); + CHECK(state_.asym_appliance[4].i(2) == doctest::Approx(0.0)); + }; + + State state; + auto main_model = default_model(state); + + state.sym_load_update[0].p_specified = 2.5e6; + + ConstDataset const update_data{ + {"sym_load", ConstDataPointer{state.sym_load_update.data(), static_cast(state.sym_load_update.size())}}, + {"asym_load", ConstDataPointer{state.asym_load_update.data(), static_cast(state.asym_load_update.size())}}, + {"shunt", ConstDataPointer{state.shunt_update.data(), static_cast(state.shunt_update.size())}}}; + + // This will lead to no topo change but param change + main_model.update_component(update_data); + + auto const math_output_sym_1 = main_model.calculate_power_flow(1e-8, 20, CalculationMethod::linear); + check_sym(main_model, math_output_sym_1); + + auto const math_output_asym_1 = main_model.calculate_power_flow(1e-8, 20, CalculationMethod::linear); + check_asym(main_model, math_output_asym_1); + + SUBCASE("No new update") { + // Math state may be fully cached + } + if constexpr (std::same_as) { + SUBCASE("No new parameter change") { + // Math state may be fully cached due to no change + main_model.update_component(update_data); + } + } + SUBCASE("With parameter change") { + // Restore to original state and re-apply same update: causes param change for cached update + main_model.restore_components(main_model.get_sequence_idx_map(update_data)); + main_model.update_component(update_data); + } + + auto const math_output_asym_2 = main_model.calculate_power_flow(1e-8, 20, CalculationMethod::linear); + check_asym(main_model, math_output_asym_2); + + auto const math_output_sym_2 = main_model.calculate_power_flow(1e-8, 20, CalculationMethod::linear); + check_sym(main_model, math_output_sym_2); + + main_model.restore_components(main_model.get_sequence_idx_map(update_data)); +} + TEST_CASE("Test main model - runtime dispatch") { State state; auto main_model = default_model(state); diff --git a/tests/cpp_unit_tests/test_y_bus.cpp b/tests/cpp_unit_tests/test_y_bus.cpp index 60a5328cd..06efe4373 100644 --- a/tests/cpp_unit_tests/test_y_bus.cpp +++ b/tests/cpp_unit_tests/test_y_bus.cpp @@ -7,6 +7,11 @@ #include +#include +#include + +#include + namespace power_grid_model { namespace { @@ -14,31 +19,28 @@ using math_solver::YBusStructure; } // namespace TEST_CASE("Test 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 - */ - + // 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 param_sym; topo.phase_shift.resize(4, 0.0); @@ -62,13 +64,12 @@ TEST_CASE("Test y bus") { // output IdxVector row_indptr = {0, 2, 5, 8, 10}; - /* Use col_indices to find the location in Y bus - * e.g. col_indices = {0, 1, 0} results in Y bus: - * [ - * x, x - * x, 0 - * ] - */ + // Use col_indices to find the location in Y bus + // e.g. col_indices = {0, 1, 0} results in Y bus: + // [ + // x, x + // x, 0 + // ] IdxVector col_indices = {// Culumn col_indices for each non-zero element in Y bus. 0, 1, 0, 1, 2, 1, 2, 3, 2, 3}; Idx nnz = 10; // Number of non-zero elements in Y bus @@ -229,18 +230,14 @@ TEST_CASE("Test one bus system") { } TEST_CASE("Test fill-in y bus") { - /* - * struct - * [1] --0--> [0] --1--> [2] - * extra fill-in: (1, 2) by removing node 0 - * - * [ - * 0, 1, 2 - * 3, 4, f - * 5, f, 6 - * ] - */ - + // [1] --0--> [0] --1--> [2] + // extra fill-in: (1, 2) by removing node 0 + // + // [ + // 0, 1, 2 + // 3, 4, f + // 5, f, 6 + // ] MathModelTopology topo{}; topo.phase_shift.resize(3, 0.0); topo.branch_bus_idx = { @@ -276,9 +273,162 @@ TEST_CASE("Test fill-in y bus") { CHECK(ybus.map_lu_y_bus == map_lu_y_bus); } -/* -TODO: -- test counting_sort_element() -*/ +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 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(topo); + + const ComplexTensorVector 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 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 param_sym_update; + // Swap based params + param_sym_update.branch_param = { + // ff, ft, tf, tt + {2.0i, 2.0i, 3.0i, 4.0i}, // (1, 0) + {5.0, 7.0, 7.0, 8.0}, // (1, 2) + {9.0i, 10.0i, 11.0i, 14.0i}, // (2, 3) + {13.0, 14.0, 17.0, 16.0}, // (3, 2) + {17.0, 18.0, 19.0, 20.0}, // (0, 1) + {1001.0i, 0.0, 0.0, 0.0} // (2,-1) + }; + param_sym_update.shunt_param = {101.0i, 200.0i}; + + const ComplexTensorVector 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 + }; + + auto verify_admittance = [](ComplexTensorVector const& admittance, + ComplexTensorVector const& admittance_ref) { + CHECK(admittance.size() == admittance_ref.size()); + for (size_t i = 0; i < admittance.size(); i++) { + CHECK(cabs(admittance[i] - admittance_ref[i]) < numerical_tolerance); + } + }; + SUBCASE("Test whole scale update") { + YBus ybus{topo_ptr, std::make_shared const>(param_sym)}; + verify_admittance(ybus.admittance(), admittance_sym); + + ybus.update_admittance(std::make_shared const>(param_sym)); + verify_admittance(ybus.admittance(), admittance_sym); + } + + SUBCASE("Test progressive update") { + YBus ybus{topo_ptr, std::make_shared const>(param_sym)}; + verify_admittance(ybus.admittance(), admittance_sym); + + auto branch_param_to_change_views = + boost::irange(0, static_cast(param_sym_update.branch_param.size())) | + boost::adaptors::filtered([¶m_sym_update](Idx i) { + return param_sym_update.branch_param[i].yff() != ComplexTensor{0.0} || + param_sym_update.branch_param[i].yft() != ComplexTensor{0.0} || + param_sym_update.branch_param[i].ytf() != ComplexTensor{0.0} || + param_sym_update.branch_param[i].ytt() != ComplexTensor{0.0}; + }); + auto shunt_param_to_change_views = boost::irange(0, static_cast(param_sym_update.shunt_param.size())) | + boost::adaptors::filtered([¶m_sym_update](Idx i) { + return param_sym_update.shunt_param[i] != ComplexTensor{0.0}; + }); + + MathModelParamIncrement math_model_param_incrmt; + 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 param_update_ptr = std::make_shared const>(param_sym_update); + + ybus.update_admittance_increment(param_update_ptr, math_model_param_incrmt); + verify_admittance(ybus.admittance(), admittance_sym_2); + } +} + +// TODO: +// - test counting_sort_element() } // namespace power_grid_model diff --git a/tests/data/power_flow/dummy-test-batch-shunt/asym_output_batch.json b/tests/data/power_flow/dummy-test-batch-shunt/asym_output_batch.json index 9082e4c58..e884d61fa 100644 --- a/tests/data/power_flow/dummy-test-batch-shunt/asym_output_batch.json +++ b/tests/data/power_flow/dummy-test-batch-shunt/asym_output_batch.json @@ -36,6 +36,38 @@ } ] }, + { + "shunt": [ + { + "id": 9, + "p": [ + 352799.852441, + 352799.852441, + 352799.852441 + ], + "q": [ + 0.0, + 0.0, + 0.0 + ], + "i": [ + 72.74613391789285, + 72.74613391789285, + 72.74613391789285 + ], + "s": [ + 352799.852441, + 352799.852441, + 352799.852441 + ], + "pf": [ + 1.0000000000000000, + 1.0000000000000000, + 1.0000000000000000 + ] + } + ] + }, { "shunt": [ { diff --git a/tests/data/power_flow/dummy-test-batch-shunt/sym_output_batch.json b/tests/data/power_flow/dummy-test-batch-shunt/sym_output_batch.json index 0dee9f9e1..0f5c97a06 100644 --- a/tests/data/power_flow/dummy-test-batch-shunt/sym_output_batch.json +++ b/tests/data/power_flow/dummy-test-batch-shunt/sym_output_batch.json @@ -16,6 +16,18 @@ } ] }, + { + "shunt": [ + { + "id": 9, + "p": 1058399.5573238870, + "q": 0.0, + "i": 72.74613391789285, + "s": 1058399.5573238870, + "pf": 1.0000000000000000 + } + ] + }, { "shunt": [ { diff --git a/tests/data/power_flow/dummy-test-batch-shunt/update_batch.json b/tests/data/power_flow/dummy-test-batch-shunt/update_batch.json index 2444f2a0c..d96dcfb3e 100644 --- a/tests/data/power_flow/dummy-test-batch-shunt/update_batch.json +++ b/tests/data/power_flow/dummy-test-batch-shunt/update_batch.json @@ -5,6 +5,14 @@ "attributes": {}, "data": [ {}, + { + "shunt": [ + { + "id": 9, + "g1": 0.015 + } + ] + }, { "shunt": [ {