diff --git a/CHANGELOG.md b/CHANGELOG.md index 2999c732..60eba3d9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -37,6 +37,7 @@ - Added support for DependencyTracking::Variable in PowerElectronics models. - Updated Jacobian value storage from `ScalarT` to `RealT`. - Added a header file defining constants to be used throughout the code. +- Added Jacobian sparsity pattern computation into `PowerElectronics`. ## v0.1 diff --git a/GridKit/LinearAlgebra/SparseMatrix/COO_Matrix.hpp b/GridKit/LinearAlgebra/SparseMatrix/COO_Matrix.hpp index 80f80ecb..1cd2f6e9 100644 --- a/GridKit/LinearAlgebra/SparseMatrix/COO_Matrix.hpp +++ b/GridKit/LinearAlgebra/SparseMatrix/COO_Matrix.hpp @@ -292,11 +292,6 @@ namespace GridKit template inline void COO_Matrix::axpy(RealT alpha, COO_Matrix& a) { - if (alpha == 0) - { - return; - } - if (!this->sorted_) { this->sortSparse(); @@ -372,9 +367,6 @@ namespace GridKit template inline void COO_Matrix::axpy(RealT alpha, std::vector r, std::vector c, std::vector v) { - if (alpha == 0) - return; - if (!this->sorted_) { this->sortSparse(); diff --git a/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.hpp b/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.hpp index 975db353..cd7e18ff 100644 --- a/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.hpp +++ b/GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.hpp @@ -76,6 +76,8 @@ namespace GridKit virtual int setValuesPointer(RealT* new_vals, memory::MemorySpace memspace); + int resize(IdxT n, IdxT m); + private: IdxT n_{0}; ///< number of rows IdxT m_{0}; ///< number of columns diff --git a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp index 5e1eee4d..88ed090e 100644 --- a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp +++ b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.cpp @@ -39,9 +39,9 @@ namespace GridKit { // internals [\delta_i, Pi, Qi, phi_di, phi_qi, gamma_di, gamma_qi, il_di, il_qi, vo_di, vo_qi, io_di, io_qi] // externals [\omega_ref, vba_out, vbb_out] - size_ = 16; - n_intern_ = 13; - n_extern_ = 3; + size_ = NUM_INTERNALS + NUM_EXTERNALS; + n_intern_ = NUM_INTERNALS; + n_extern_ = NUM_EXTERNALS; extern_indices_ = {0, 1, 2}; idc_ = id; } diff --git a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp index d2c930ac..8b54c9aa 100644 --- a/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp +++ b/GridKit/Model/PowerElectronics/DistributedGenerator/DistributedGenerator.hpp @@ -64,6 +64,9 @@ namespace GridKit using CircuitComponent::n_intern_; public: + static constexpr size_t NUM_INTERNALS = 13; + static constexpr size_t NUM_EXTERNALS = 3; + DistributedGenerator(IdxT id, DistributedGeneratorParameters parm, bool reference_frame); diff --git a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp index fd5dbe8f..2f159448 100644 --- a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp +++ b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.cpp @@ -24,9 +24,9 @@ namespace GridKit : RN_(RN) { // externals [vbus_d, vbus_q] - size_ = 2; - n_intern_ = 0; - n_extern_ = 2; + size_ = NUM_INTERNALS + NUM_EXTERNALS; + n_intern_ = NUM_INTERNALS; + n_extern_ = NUM_EXTERNALS; extern_indices_ = {0, 1}; idc_ = id; } diff --git a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp index 8c0d8dad..646e56fd 100644 --- a/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp +++ b/GridKit/Model/PowerElectronics/MicrogridBusDQ/MicrogridBusDQ.hpp @@ -43,6 +43,9 @@ namespace GridKit using CircuitComponent::n_intern_; public: + static constexpr size_t NUM_INTERNALS = 0; + static constexpr size_t NUM_EXTERNALS = 2; + MicrogridBusDQ(IdxT id, RealT RN); virtual ~MicrogridBusDQ(); diff --git a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp index 44f9073d..b5a8f6b8 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp +++ b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.cpp @@ -28,9 +28,9 @@ namespace GridKit { // internals [id, iq] // externals [\omegaref, vbd_in, vbq_in, vbd_out, vbq_out] - size_ = 7; - n_intern_ = 2; - n_extern_ = 5; + size_ = NUM_INTERNALS + NUM_EXTERNALS; + n_intern_ = NUM_INTERNALS; + n_extern_ = NUM_EXTERNALS; extern_indices_ = {0, 1, 2, 3, 4}; idc_ = id; } diff --git a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp index 64992c5c..1f646251 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp +++ b/GridKit/Model/PowerElectronics/MicrogridLine/MicrogridLine.hpp @@ -44,6 +44,9 @@ namespace GridKit using CircuitComponent::n_intern_; public: + static constexpr size_t NUM_INTERNALS = 2; + static constexpr size_t NUM_EXTERNALS = 5; + MicrogridLine(IdxT id, RealT R, RealT L); virtual ~MicrogridLine(); diff --git a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp index e16b70b6..2ca50f90 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp +++ b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.cpp @@ -26,9 +26,9 @@ namespace GridKit { // internals [id, iq] // externals [\omegaref, vbd_out, vbq_out] - size_ = 5; - n_intern_ = 2; - n_extern_ = 3; + size_ = NUM_INTERNALS + NUM_EXTERNALS; + n_intern_ = NUM_INTERNALS; + n_extern_ = NUM_EXTERNALS; extern_indices_ = {0, 1, 2}; idc_ = id; } diff --git a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp index 0a8f34bb..b8e0432c 100644 --- a/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp +++ b/GridKit/Model/PowerElectronics/MicrogridLoad/MicrogridLoad.hpp @@ -44,6 +44,9 @@ namespace GridKit using CircuitComponent::n_intern_; public: + static constexpr size_t NUM_INTERNALS = 2; + static constexpr size_t NUM_EXTERNALS = 3; + MicrogridLoad(IdxT id, RealT R, RealT L); virtual ~MicrogridLoad(); diff --git a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp index 736f9d00..98e92243 100644 --- a/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp +++ b/GridKit/Model/PowerElectronics/SystemModelPowerElectronics.hpp @@ -2,11 +2,15 @@ #pragma once +#include #include +#include #include #include #include +#include +#include #include #include #include @@ -60,8 +64,11 @@ namespace GridKit { using RealT = typename CircuitComponent::RealT; using component_type = CircuitComponent; + using MatrixT = CircuitComponent::MatrixT; using CircuitComponent::size_; + using CircuitComponent::n_intern_; + using CircuitComponent::n_extern_; using CircuitComponent::nnz_; using CircuitComponent::time_; using CircuitComponent::alpha_; @@ -71,6 +78,7 @@ namespace GridKit using CircuitComponent::jac_; using CircuitComponent::rel_tol_; using CircuitComponent::abs_tol_; + using CircuitComponent::max_steps_; public: /** @@ -81,11 +89,11 @@ namespace GridKit PowerElectronicsModel() { // Set system model parameters as default - rel_tol_ = 1e-4; - abs_tol_ = 1e-4; - this->max_steps_ = 2000; + rel_tol_ = 1e-4; + abs_tol_ = 1e-4; + max_steps_ = 2000; // By default don't use the jacobian - use_jac_ = false; + use_jac_ = false; } /** @@ -104,11 +112,12 @@ namespace GridKit IdxT max_steps = 2000) { // Set system model tolerances from input - rel_tol_ = rel_tol; - abs_tol_ = abs_tol; - this->max_steps_ = max_steps; + rel_tol_ = rel_tol; + abs_tol_ = abs_tol; + max_steps_ = max_steps; // Can choose if to use jacobian - use_jac_ = use_jac; + use_jac_ = use_jac; + csr_jac_ = nullptr; } /** @@ -121,20 +130,186 @@ namespace GridKit */ virtual ~PowerElectronicsModel() { - for (auto comp : this->components_) + for (auto comp : components_) delete comp; + + if (csr_jac_ != nullptr) + delete csr_jac_; } /** - * @brief allocator default - * - * @todo this should throw an exception as no allocation without a graph is allowed. - * Or needs to be removed from the base class - * - * @return int + * @brief Allocates and constructs sparsity patterns for system jacobian. + * To do this, sparsity patterns of components are needed, so + * each component's jacobian is evaluated, and the sparsity pattern + * of that component's jacobian is expected to match the sparsity pattern + * for the rest of the simulation. + * @pre `size_` must be correctly initialized */ int allocate() { + // Helper struct for identifying a particular component's local system variable + struct ComponentContribution + { + // The index of the component in `components_`. + size_t comp_idx_; + // The local system variable index + size_t local_row_idx_; + }; + + // A reverse mapping from external system variables -> component variables + std::forward_list* reverse_extern_map = new std::forward_list[n_extern_]; + size_t component_nnz = 0; + + // Loop over all components, evaluate their jacobians, save their sparsity information, + // and construct the reverse variable mapping. + for (size_t comp_idx = 0; comp_idx < components_.size(); comp_idx++) + { + component_type* component = components_[comp_idx]; + component->evaluateJacobian(); + + MatrixT& comp_jacobian = component->getJacobian(); + + for (IdxT local_external_row : component->getExternIndices()) + { + IdxT global_row = component->getNodeConnection(local_external_row); + + // Not all local variables map to a global variable + if (global_row != neg1_) + { + // global_row >= n_intern_ (see allocate(s)) + reverse_extern_map[global_row - n_intern_].push_front({.comp_idx_ = comp_idx, .local_row_idx_ = local_external_row}); + } + } + + component_nnz += comp_jacobian.nnz(); + } + + // Allocate the final sparsity pattern info. Not deleted, because ownership is stolen by the jacobian + IdxT* global_row_indices = new IdxT[size_ + 1]; + IdxT* global_col_indices = new IdxT[component_nnz]; // Use component_nnz as an upper bound on nnz + global_row_indices[0] = 0; + + // Construct Jacobian sparsity pattern + + // Start with internal rows, which must be before external rows, are grouped by component, + // strictly increasing wrt component internals, and must be sorted by component (see allocate(s)) + size_t curr_internal_row = 0; + for (size_t comp_idx = 0; comp_idx < components_.size(); comp_idx++) + { + component_type* component = components_[comp_idx]; + std::set comp_externals = component->getExternIndices(); + std::tuple&, std::vector&, std::vector&> entries = component->getJacobian().getEntries(); + const std::vector& rows = std::get<0>(entries); + const std::vector& columns = std::get<1>(entries); + + if (rows.empty()) + continue; + + // Loop through all elements of the component jacobian, adding them to the system jacobian if needed + for (size_t elem_idx = 0; elem_idx < rows.size(); elem_idx++) + { + IdxT local_row = rows[elem_idx]; + IdxT global_row = component->getNodeConnection(local_row); + + // Skip variables which aren't system variables + if (comp_externals.contains(local_row) || global_row == neg1_) + continue; + + if (global_row > curr_internal_row) + { + curr_internal_row++; + global_row_indices[curr_internal_row + 1] = global_row_indices[curr_internal_row]; + } + + assert(global_row == curr_internal_row); + + IdxT local_col_idx = columns[elem_idx]; + IdxT global_col_idx = component->getNodeConnection(local_col_idx); + + if (global_col_idx == neg1_) + continue; + + global_col_indices[global_row_indices[curr_internal_row + 1]] = global_col_idx; + global_row_indices[curr_internal_row + 1]++; + } + } + + // One last row after + curr_internal_row++; + global_row_indices[curr_internal_row + 1] = global_row_indices[curr_internal_row]; + + // Need to sort internal rows because even though the mapping from component internal to system internal + // is monotonic, the mapping from component external to system external may not be, and internal rows + // may contain external columns. + for (size_t row_idx = 0; row_idx < curr_internal_row; row_idx++) + { + IdxT* global_row_start = global_col_indices + global_row_indices[row_idx]; + IdxT* global_row_end = global_col_indices + global_row_indices[row_idx + 1]; + std::sort(global_row_start, global_row_end); + } + + assert(curr_internal_row == n_intern_); + + // Then move on to external rows, which come after internal rows + for (size_t row = n_intern_; row < size_; row++) + { + global_row_indices[row + 1] = global_row_indices[row]; + + // Collect columns from each component which has a row which contributes to this row + for (ComponentContribution contrib : reverse_extern_map[row - n_intern_]) + { + component_type* component = components_[contrib.comp_idx_]; + std::tuple&, std::vector&, std::vector&> entries = component->getJacobian().getEntries(); + const std::vector& rows = std::get<0>(entries); + const std::vector& columns = std::get<1>(entries); + + typename std::vector::const_iterator row_start = std::ranges::lower_bound(rows, contrib.local_row_idx_); + + // It can happen where external contributions are only constants, and do not appear in the jacobian. + // If that is the case, we won't be able to find local_row_idx_ and must skip this contribution + if (row_start == rows.end() || *row_start != contrib.local_row_idx_) + continue; + + typename std::vector::const_iterator row_end = std::upper_bound(row_start, rows.end(), contrib.local_row_idx_); + + for (size_t local_elem_idx = std::distance(rows.begin(), row_start); + local_elem_idx < static_cast(std::distance(rows.begin(), row_end)); + local_elem_idx++) + { + IdxT local_col = columns[local_elem_idx]; + IdxT global_col = component->getNodeConnection(local_col); + + // Not all columns map to columns in the system jacobian + if (global_col != neg1_) + { + global_col_indices[global_row_indices[row + 1]] = global_col; + global_row_indices[row + 1]++; + } + } + } + + // Sort the row by column indices. Since the mapping from local indices to global indices isn't monotonically increasing, + // this is necessary. + IdxT* start = global_col_indices + global_row_indices[row]; + IdxT* end = global_col_indices + global_row_indices[row + 1]; + std::sort(start, end); + + // De-duplicate the columns + IdxT* new_end = std::unique(start, end); + global_row_indices[row + 1] = global_row_indices[row] + static_cast(std::distance(start, new_end)); + } + // Allocate new sparsity buffers + IdxT nnz = global_row_indices[size_]; + + if (csr_jac_ != nullptr) + delete csr_jac_; + + RealT* value_buffer = new RealT[nnz]; + csr_jac_ = new LinearAlgebra::CsrMatrix( + size_, size_, nnz, &global_row_indices, &global_col_indices, &value_buffer, LinearAlgebra::memory::HOST, LinearAlgebra::memory::HOST); + + delete[] reverse_extern_map; + return 1; } @@ -147,7 +322,7 @@ namespace GridKit */ bool hasJacobian() { - if (!this->use_jac_) + if (!use_jac_) return false; for (const auto& component : components_) @@ -173,10 +348,74 @@ namespace GridKit int allocate(IdxT s) { // Allocate all components - size_ = s; + size_ = s; + n_intern_ = 0; for (const auto& component : components_) { component->allocate(); + + // Count up the amount of internal variables which get mapped to system variables + std::set extern_indices = component->getExternIndices(); + for (IdxT comp_var_idx = 0; comp_var_idx < component->size(); comp_var_idx++) + { + IdxT sys_var_idx = component->getNodeConnection(comp_var_idx); + if (!extern_indices.contains(comp_var_idx) && sys_var_idx != neg1_) + { + n_intern_++; + } + } + } + n_extern_ = size_ - n_intern_; + + // Ensure that all component locals are mapped to system locals + // and all component externals are mapped to system externals. + // System locals are stored first, in 0..n_intern_ and externals + // are stored after, in n_intern_.. + // Additionally, ensure that components locals are grouped in the system vectors - no other + // variables are between locals from a single component. As well, ensure that components are + // sorted by these groupings, so the first component is the first block and so on. + for (size_t comp_idx = 0; comp_idx < components_.size(); comp_idx++) + { + component_type* component = components_[comp_idx]; + std::set extern_indices = component->getExternIndices(); + + // Whether or not we've seen a local variable yet + bool has_seen_local = false; + // The system index of the last local we've seen + IdxT last_local_sys; + + for (IdxT comp_var_idx = 0; comp_var_idx < component->size(); comp_var_idx++) + { + IdxT sys_var_idx = component->getNodeConnection(comp_var_idx); + + // Ignore local variables which aren't mapped to sytem variables + if (sys_var_idx == neg1_) + continue; + + if (extern_indices.contains(comp_var_idx)) + { + assert(sys_var_idx >= n_intern_); + } + else + { + assert(sys_var_idx < n_intern_); + + // Make sure that all of the locals for a particular component are grouped + // together in the sytem vector, and have increasing indices. + if (has_seen_local) + { + assert(sys_var_idx == last_local_sys + 1); + } + else if (comp_idx > 0) + { + // Ensure increasing indices in components - so no need to sort components + assert(sys_var_idx > last_local_sys); + } + + has_seen_local = true; + last_local_sys = sys_var_idx; + } + } } // Allocate global vectors @@ -184,6 +423,8 @@ namespace GridKit yp_.resize(size_); f_.resize(size_); + allocate(); + return 0; } @@ -199,7 +440,7 @@ namespace GridKit { component->initialize(); } - this->distributeVectors(); + distributeVectors(); return 0; } @@ -244,12 +485,12 @@ namespace GridKit */ int evaluateResidual() { - for (IdxT i = 0; i < this->f_.size(); i++) + for (IdxT i = 0; i < f_.size(); i++) { f_[i] = 0.0; } - this->distributeVectors(); + distributeVectors(); // Update system residual vector @@ -401,9 +642,20 @@ namespace GridKit components_.push_back(component); } + /** + * @brief Returns a reference to the CSR Jacobian. + * @todo Can be removed once `getJacobian()` returns a `CsrMatrix` + */ + LinearAlgebra::CsrMatrix& getCsrJac() + { + return *csr_jac_; + } + private: static constexpr IdxT neg1_ = static_cast(-1); + LinearAlgebra::CsrMatrix* csr_jac_; + std::vector components_; int jac_call_count_{0}; diff --git a/examples/PowerElectronics/DistributedGeneratorTest/CMakeLists.txt b/examples/PowerElectronics/DistributedGeneratorTest/CMakeLists.txt index 225dd73f..d17aa846 100644 --- a/examples/PowerElectronics/DistributedGeneratorTest/CMakeLists.txt +++ b/examples/PowerElectronics/DistributedGeneratorTest/CMakeLists.txt @@ -5,7 +5,8 @@ add_executable(dgtest DGTest.cpp) target_link_libraries(dgtest GridKit::power_elec_disgen GridKit::power_elec_microline - GridKit::power_elec_microload) + GridKit::power_elec_microload + GridKit::CsrMatrix) add_test(NAME DistributedGeneratorTest COMMAND $) install(TARGETS dgtest RUNTIME DESTINATION bin) diff --git a/examples/PowerElectronics/Microgrid/CMakeLists.txt b/examples/PowerElectronics/Microgrid/CMakeLists.txt index 1c67f74f..ef82d2e4 100644 --- a/examples/PowerElectronics/Microgrid/CMakeLists.txt +++ b/examples/PowerElectronics/Microgrid/CMakeLists.txt @@ -7,7 +7,8 @@ target_link_libraries(microgrid GridKit::power_elec_disgen GridKit::power_elec_microline GridKit::power_elec_microload GridKit::solvers_dyn - GridKit::power_elec_microbusdq) + GridKit::power_elec_microbusdq + GridKit::CsrMatrix) add_test(NAME Microgrid COMMAND $) install(TARGETS microgrid RUNTIME DESTINATION bin) diff --git a/examples/PowerElectronics/RLCircuit/CMakeLists.txt b/examples/PowerElectronics/RLCircuit/CMakeLists.txt index 42ffce87..52af279f 100644 --- a/examples/PowerElectronics/RLCircuit/CMakeLists.txt +++ b/examples/PowerElectronics/RLCircuit/CMakeLists.txt @@ -7,7 +7,8 @@ target_link_libraries(rlcircuit GridKit::power_elec_capacitor GridKit::power_elec_inductor GridKit::power_elec_resistor GridKit::power_elec_voltagesource - GridKit::solvers_dyn) + GridKit::solvers_dyn + GridKit::CsrMatrix) add_test(NAME RLCircuit COMMAND $) install(TARGETS rlcircuit RUNTIME DESTINATION bin) diff --git a/examples/PowerElectronics/RLCircuit/RLCircuit.cpp b/examples/PowerElectronics/RLCircuit/RLCircuit.cpp index dce13668..733aabbe 100644 --- a/examples/PowerElectronics/RLCircuit/RLCircuit.cpp +++ b/examples/PowerElectronics/RLCircuit/RLCircuit.cpp @@ -31,15 +31,21 @@ int main(int /* argc */, char const** /* argv */) double linit = 1.0; double vinit = 1.0; + // System vector map + // 0 -> Inductor internal + // 1 -> VSource internal + // 2 -> Resistor/VSource external + // 3 -> Inductor/Resistor external + // inductor GridKit::Inductor* induct = new GridKit::Inductor(idoff, linit); // Form index to node uid realations // input - induct->setExternalConnectionNodes(0, 1); + induct->setExternalConnectionNodes(0, 3); // output induct->setExternalConnectionNodes(1, static_cast(-1)); // internal - induct->setExternalConnectionNodes(2, 2); + induct->setExternalConnectionNodes(2, 0); // add component sysmodel.addComponent(induct); @@ -48,9 +54,9 @@ int main(int /* argc */, char const** /* argv */) GridKit::Resistor* resis = new GridKit::Resistor(idoff, rinit); // Form index to node uid realations // input - resis->setExternalConnectionNodes(0, 0); + resis->setExternalConnectionNodes(0, 2); // output - resis->setExternalConnectionNodes(1, 1); + resis->setExternalConnectionNodes(1, 3); // add sysmodel.addComponent(resis); @@ -61,9 +67,9 @@ int main(int /* argc */, char const** /* argv */) // input vsource->setExternalConnectionNodes(0, static_cast(-1)); // output - vsource->setExternalConnectionNodes(1, 0); + vsource->setExternalConnectionNodes(1, 2); // internal - vsource->setExternalConnectionNodes(2, 3); + vsource->setExternalConnectionNodes(2, 1); sysmodel.addComponent(vsource); @@ -74,15 +80,15 @@ int main(int /* argc */, char const** /* argv */) // Grounding for IDA. If no grounding then circuit is \mu > 1 // v_0 (grounded) // Create Intial points - sysmodel.y()[0] = vinit; // v_1 - sysmodel.y()[1] = vinit; // v_2 - sysmodel.y()[2] = 0.0; // i_L - sysmodel.y()[3] = 0.0; // i_s + sysmodel.y()[0] = 0.0; // i_L + sysmodel.y()[1] = 0.0; // i_s + sysmodel.y()[2] = vinit; // v_1 + sysmodel.y()[3] = vinit; // v_2 - sysmodel.yp()[0] = 0.0; // v'_1 - sysmodel.yp()[1] = 0.0; // v'_2 - sysmodel.yp()[2] = -vinit / linit; // i'_s - sysmodel.yp()[3] = -vinit / linit; // i'_L + sysmodel.yp()[0] = -vinit / linit; // i'_s + sysmodel.yp()[1] = -vinit / linit; // i'_L + sysmodel.yp()[2] = 0.0; // v'_1 + sysmodel.yp()[3] = 0.0; // v'_2 sysmodel.initialize(); sysmodel.evaluateResidual(); @@ -123,10 +129,10 @@ int main(int /* argc */, char const** /* argv */) std::vector yexact(4); // analytical solution to the circuit - yexact[0] = vinit; - yexact[2] = (vinit / rinit) * (exp(-(rinit / linit) * t_final) - 1.0); - yexact[3] = yexact[2]; - yexact[1] = vinit + rinit * yexact[2]; + yexact[2] = vinit; + yexact[0] = (vinit / rinit) * (exp(-(rinit / linit) * t_final) - 1.0); + yexact[1] = yexact[0]; + yexact[3] = vinit + rinit * yexact[0]; std::cout << "Element-wise Relative error at t=" << t_final << "\n"; for (size_t i = 0; i < yfinial.size(); i++) diff --git a/examples/PowerElectronics/ScaleMicrogrid/CMakeLists.txt b/examples/PowerElectronics/ScaleMicrogrid/CMakeLists.txt index aeea5f01..c711d0d6 100644 --- a/examples/PowerElectronics/ScaleMicrogrid/CMakeLists.txt +++ b/examples/PowerElectronics/ScaleMicrogrid/CMakeLists.txt @@ -8,13 +8,15 @@ target_link_libraries(scalemicrogrid GridKit::power_elec_disgen GridKit::power_elec_microline GridKit::power_elec_microload GridKit::solvers_dyn - GridKit::power_elec_microbusdq) + GridKit::power_elec_microbusdq + GridKit::CsrMatrix) target_link_libraries(scalemicrogridarbitrary GridKit::power_elec_disgen GridKit::power_elec_microline GridKit::power_elec_microload GridKit::solvers_dyn - GridKit::power_elec_microbusdq) + GridKit::power_elec_microbusdq + GridKit::CsrMatrix) add_test(NAME ScaleMicrogrid COMMAND $) install(TARGETS scalemicrogrid RUNTIME DESTINATION bin) diff --git a/tests/UnitTests/CMakeLists.txt b/tests/UnitTests/CMakeLists.txt index c002660b..3c7e0892 100644 --- a/tests/UnitTests/CMakeLists.txt +++ b/tests/UnitTests/CMakeLists.txt @@ -2,5 +2,6 @@ add_subdirectory(AutomaticDifferentiation) add_subdirectory(LinearAlgebra) add_subdirectory(PhasorDynamics) +add_subdirectory(PowerElectronics) add_subdirectory(Solver) add_subdirectory(Utilities) diff --git a/tests/UnitTests/PowerElectronics/CMakeLists.txt b/tests/UnitTests/PowerElectronics/CMakeLists.txt new file mode 100644 index 00000000..24e1905c --- /dev/null +++ b/tests/UnitTests/PowerElectronics/CMakeLists.txt @@ -0,0 +1,11 @@ +add_executable(test_power_elec_system runSystemTests.cpp) +target_link_libraries(test_power_elec_system GridKit::power_elec_disgen + GridKit::power_elec_microline + GridKit::power_elec_microload + GridKit::power_elec_microbusdq + GridKit::CsrMatrix) + +add_test(NAME PowerElectronicsSystemTest COMMAND test_power_elec_system) + +install(TARGETS test_power_elec_system + RUNTIME DESTINATION bin) diff --git a/tests/UnitTests/PowerElectronics/SystemTests.hpp b/tests/UnitTests/PowerElectronics/SystemTests.hpp new file mode 100644 index 00000000..44780a46 --- /dev/null +++ b/tests/UnitTests/PowerElectronics/SystemTests.hpp @@ -0,0 +1,312 @@ +#pragma once + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + + template + class SystemTests + { + public: + SystemTests() = default; + ~SystemTests() = default; + + TestOutcome jacobianSparsity(size_t Nsize) + { + TestStatus success = true; + + double abs_tol = 1.0e-8; + double rel_tol = 1.0e-8; + size_t max_step_number = 3000; + bool use_jac = true; + + // Create model + auto* sysmodel = new GridKit::PowerElectronicsModel(rel_tol, abs_tol, use_jac, max_step_number); + + // Modeled after the problem in the paper + double RN = 1.0e4; + + // DG Params Vector + // All DGs have the same set of parameters except for the first two. + GridKit::DistributedGeneratorParameters DG_parms1; + DG_parms1.wb_ = 2.0 * M_PI * 50.0; + DG_parms1.wc_ = 31.41; + DG_parms1.mp_ = 9.4e-5; + DG_parms1.Vn_ = 380.0; + DG_parms1.nq_ = 1.3e-3; + DG_parms1.F_ = 0.75; + DG_parms1.Kiv_ = 420.0; + DG_parms1.Kpv_ = 0.1; + DG_parms1.Kic_ = 2.0e4; + DG_parms1.Kpc_ = 15.0; + DG_parms1.Cf_ = 5.0e-5; + DG_parms1.rLf_ = 0.1; + DG_parms1.Lf_ = 1.35e-3; + DG_parms1.rLc_ = 0.03; + DG_parms1.Lc_ = 0.35e-3; + + GridKit::DistributedGeneratorParameters DG_parms2; + DG_parms2.wb_ = 2.0 * M_PI * 50.0; + DG_parms2.wc_ = 31.41; + DG_parms2.mp_ = 12.5e-5; + DG_parms2.Vn_ = 380.0; + DG_parms2.nq_ = 1.5e-3; + DG_parms2.F_ = 0.75; + DG_parms2.Kiv_ = 390.0; + DG_parms2.Kpv_ = 0.05; + DG_parms2.Kic_ = 16.0e3; + DG_parms2.Kpc_ = 10.5; + DG_parms2.Cf_ = 50.0e-6; + DG_parms2.rLf_ = 0.1; + DG_parms2.Lf_ = 1.35e-3; + DG_parms2.rLc_ = 0.03; + DG_parms2.Lc_ = 0.35e-3; + + std::vector> DGParams_list(2 * Nsize, DG_parms2); + + DGParams_list[0] = DG_parms1; + DGParams_list[1] = DG_parms1; + + // line vector params + // Every odd line has the same parameters and every even line has the same parameters + double rline1 = 0.23; + double Lline1 = 0.1 / (2.0 * M_PI * 50.0); + double rline2 = 0.35; + double Lline2 = 0.58 / (2.0 * M_PI * 50.0); + std::vector rline_list(2 * Nsize - 1, 0.0); + std::vector Lline_list(2 * Nsize - 1, 0.0); + for (size_t i = 0; i < rline_list.size(); i++) + { + rline_list[i] = (i % 2) ? rline2 : rline1; + Lline_list[i] = (i % 2) ? Lline2 : Lline1; + } + + // load parms + // Only the first load has the same paramaters. + double rload1 = 3.0; + double Lload1 = 2.0 / (2.0 * M_PI * 50.0); + double rload2 = 2.0; + double Lload2 = 1.0 / (2.0 * M_PI * 50.0); + + std::vector rload_list(Nsize, rload2); + std::vector Lload_list(Nsize, Lload2); + rload_list[0] = rload1; + Lload_list[0] = Lload1; + + size_t num_dgs = 2 * Nsize; + size_t num_lines = 2 * Nsize - 1; + size_t num_loads = Nsize; + size_t num_buses = 2 * Nsize; + + size_t vec_size_internals = + // - refframe + num_dgs * DistributedGenerator::NUM_INTERNALS - 1 + + num_lines * MicrogridLine::NUM_INTERNALS + + num_loads * MicrogridLoad::NUM_INTERNALS + + num_buses * MicrogridBusDQ::NUM_INTERNALS; + + size_t vec_size_externals = + // + refframe + num_buses * MicrogridBusDQ::NUM_EXTERNALS + 1; + + std::vector vdqbus_index(2 * Nsize, 0); + vdqbus_index[0] = vec_size_internals + 1; + for (size_t i = 1; i < vdqbus_index.size(); i++) + { + vdqbus_index[i] = vdqbus_index[i - 1] + 2; + } + + // Total size of the vector setup + size_t vec_size_total = vec_size_internals + vec_size_externals; + + // Create the reference DG + auto* dg_ref = new DistributedGenerator(0, + DGParams_list[0], + true); + // ref motor + dg_ref->setExternalConnectionNodes(0, vec_size_internals); + // outputs + dg_ref->setExternalConnectionNodes(1, vdqbus_index[0]); + dg_ref->setExternalConnectionNodes(2, vdqbus_index[0] + 1); + //"grounding" of the difference + dg_ref->setExternalConnectionNodes(3, static_cast(-1)); + // internal connections + for (size_t i = 0; i < dg_ref->NUM_INTERNALS - 1; i++) + { + dg_ref->setExternalConnectionNodes(4 + i, i); + } + sysmodel->addComponent(dg_ref); + + // Keep track of models and index location + size_t indexv = dg_ref->NUM_INTERNALS - 1; + size_t model_id = 1; + // Add all other DGs + for (size_t i = 1; i < 2 * Nsize; i++) + { + // current DG to add + auto* dg = new DistributedGenerator(model_id++, + DGParams_list[i], + false); + // ref motor + dg->setExternalConnectionNodes(0, vec_size_internals); + // outputs + dg->setExternalConnectionNodes(1, vdqbus_index[i]); + dg->setExternalConnectionNodes(2, vdqbus_index[i] + 1); + // internal connections + for (size_t j = 0; j < dg_ref->NUM_INTERNALS; j++) + { + dg->setExternalConnectionNodes(3 + j, indexv + j); + } + indexv += dg_ref->NUM_INTERNALS; + sysmodel->addComponent(dg); + } + + // Load all the Line compoenents + for (size_t i = 0; i < 2 * Nsize - 1; i++) + { + // line + auto* line_model = new MicrogridLine(model_id++, + rline_list[i], + Lline_list[i]); + // ref motor + line_model->setExternalConnectionNodes(0, vec_size_internals); + // input connections + line_model->setExternalConnectionNodes(1, vdqbus_index[i]); + line_model->setExternalConnectionNodes(2, vdqbus_index[i] + 1); + // output connections + line_model->setExternalConnectionNodes(3, vdqbus_index[i + 1]); + line_model->setExternalConnectionNodes(4, vdqbus_index[i + 1] + 1); + // internal connections + for (size_t j = 0; j < 2; j++) + { + line_model->setExternalConnectionNodes(5 + j, indexv + j); + } + indexv += 2; + sysmodel->addComponent(line_model); + } + + // Load all the Load components + for (size_t i = 0; i < Nsize; i++) + { + auto* load_model = new MicrogridLoad(model_id++, + rload_list[i], + Lload_list[i]); + // ref motor + load_model->setExternalConnectionNodes(0, vec_size_internals); + // input connections + load_model->setExternalConnectionNodes(1, vdqbus_index[2 * i]); + load_model->setExternalConnectionNodes(2, vdqbus_index[2 * i] + 1); + // internal connections + for (size_t j = 0; j < 2; j++) + { + load_model->setExternalConnectionNodes(3 + j, indexv + j); + } + indexv += 2; + sysmodel->addComponent(load_model); + } + + // Add all the microgrid Virtual DQ Buses + for (size_t i = 0; i < 2 * Nsize; i++) + { + auto* virDQbus_model = new MicrogridBusDQ(model_id++, RN); + + virDQbus_model->setExternalConnectionNodes(0, vdqbus_index[i]); + virDQbus_model->setExternalConnectionNodes(1, vdqbus_index[i] + 1); + sysmodel->addComponent(virDQbus_model); + } + + sysmodel->allocate(vec_size_total); + + // Create Intial points for states + for (size_t i = 0; i < vec_size_total; i++) + { + sysmodel->y()[i] = 0.0; + sysmodel->yp()[i] = 0.0; + sysmodel->y()[i].setVariableNumber(i); + } + + // Create Intial derivatives specifics generated in MATLAB + for (size_t i = 0; i < 2 * Nsize; i++) + { + sysmodel->yp()[dg_ref->NUM_INTERNALS * i - 1 + 3] = DGParams_list[i].Vn_; + sysmodel->yp()[dg_ref->NUM_INTERNALS * i - 1 + 5] = DGParams_list[i].Kpv_ * DGParams_list[i].Vn_; + sysmodel->yp()[dg_ref->NUM_INTERNALS * i - 1 + 7] = (DGParams_list[i].Kpc_ * DGParams_list[i].Kpv_ * DGParams_list[i].Vn_) / DGParams_list[i].Lf_; + } + + // since the intial P_com = 0, the set the intial vector to the reference frame + sysmodel->y()[vec_size_internals] = DG_parms1.wb_; + sysmodel->y()[vec_size_internals].setVariableNumber(vec_size_internals); + + for (size_t i = 0; i < vec_size_total; i++) + { + sysmodel->yp()[i].setVariableNumber(i); + } + + sysmodel->initialize(); + sysmodel->evaluateResidual(); + + std::vector& residuals = sysmodel->getResidual(); + + size_t* row_indices = sysmodel->getCsrJac().getRowData(); + size_t* col_indices = sysmodel->getCsrJac().getColData(); + + for (size_t row = 0; row < residuals.size(); row++) + { + bool row_correct = true; + + DependencyTracking::Variable& row_residual = residuals[row]; + auto& dependency_map = row_residual.getDependencies(); + auto dependencies = std::vector(dependency_map.begin(), dependency_map.end()); + + // Verify both rows are the same size + row_correct &= (row_indices[row + 1] - row_indices[row]) == dependencies.size(); + + if (row_correct) + { + for (size_t col_idx = 0; col_idx < dependencies.size(); col_idx++) + { + row_correct &= col_indices[col_idx + row_indices[row]] == std::get<0>(dependencies[col_idx]); + } + } + + if (!row_correct) + { + std::cerr << "Row " << row << " is incorrect:\n"; + std::cerr << " CSR Row: {"; + for (size_t col_idx = row_indices[row]; col_idx < row_indices[row + 1]; col_idx++) + { + std::cerr << std::format("{:3}{}", col_indices[col_idx], col_idx == row_indices[row + 1] - 1 ? "" : ", "); + } + std::cerr << "}\n" + << " Dep Tracker Row: {"; + for (size_t col_idx = 0; col_idx < dependencies.size(); col_idx++) + { + std::cerr << std::format("{:3}{}", std::get<0>(dependencies[col_idx]), col_idx == dependencies.size() - 1 ? "" : ", "); + } + std::cerr << "}\n\n"; + } + + success *= row_correct; + } + + delete sysmodel; + + std::string test_name = __func__ + std::format(" (Nsize = {})", Nsize); + return success.report(test_name.c_str()); + } + }; + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/PowerElectronics/runSystemTests.cpp b/tests/UnitTests/PowerElectronics/runSystemTests.cpp new file mode 100644 index 00000000..de40ecfd --- /dev/null +++ b/tests/UnitTests/PowerElectronics/runSystemTests.cpp @@ -0,0 +1,15 @@ +#include "SystemTests.hpp" + +int main() +{ + using namespace GridKit; + using namespace GridKit::Testing; + + GridKit::Testing::TestingResults result; + GridKit::Testing::SystemTests test; + + result += test.jacobianSparsity(2); + result += test.jacobianSparsity(128); + + return result.summary(); +}