diff --git a/CMake/FindSuiteSparse.cmake b/CMake/FindSuiteSparse.cmake index 2d73feb7..5a102693 100644 --- a/CMake/FindSuiteSparse.cmake +++ b/CMake/FindSuiteSparse.cmake @@ -69,11 +69,11 @@ Author(s): set(SUITESPARSE_MODULES amd colamd - klu) + klu + suitesparseconfig) find_library(SUITESPARSE_LIBRARY NAMES - suitesparseconfig ${SUITESPARSE_MODULES} PATHS ${SUITESPARSE_DIR} $ENV{SUITESPARSE_DIR} ${SUITESPARSE_ROOT_DIR} @@ -97,6 +97,7 @@ find_path(SUITESPARSE_INCLUDE_DIR amd.h colamd.h klu.h + SuiteSparse_config.h PATHS ${SUITESPARSE_DIR} $ENV{SUITESPARSE_DIR} ${SUITESPARSE_ROOT_DIR} ${SUITESPARSE_LIBRARY_DIR}/.. PATH_SUFFIXES diff --git a/CMakeLists.txt b/CMakeLists.txt index 19048d37..badea53a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -81,7 +81,7 @@ option(GRIDKIT_ENABLE_IPOPT "Enable Ipopt support" ON) option(GRIDKIT_ENABLE_SUNDIALS "Enable SUNDIALS support" ON) # Enable KLU -option(GRIDKIT_ENABLE_SUNDIALS_SPARSE "Enable SUNDIALS sparse linear solvers" OFF) +option(GRIDKIT_ENABLE_SUNDIALS_SPARSE "Enable SUNDIALS sparse linear solvers" ON) set(CMAKE_MACOSX_RPATH 1) @@ -110,7 +110,7 @@ endif("${isSystemDir}" STREQUAL "-1") # TODO: Probably beter to set a debug interface target -set(CMAKE_CXX_FLAGS_DEBUG "-Wall -O0 -g") +set(CMAKE_CXX_FLAGS_DEBUG "-Wall -O0 -g -DDEBUG") set(CMAKE_CXX_STANDARD 17) @@ -143,6 +143,9 @@ add_subdirectory(ComponentLib) # General Utilities and File IO add_subdirectory(Utilities) +#Local Sparse matrix operations +add_subdirectory(SparseMatrix) + # Create solvers add_subdirectory(Solver) diff --git a/CircuitGraph.hpp b/CircuitGraph.hpp new file mode 100644 index 00000000..070d838f --- /dev/null +++ b/CircuitGraph.hpp @@ -0,0 +1,122 @@ + + +#include +#include +#include +#include +#include + +/** + * @brief A very basic hypergraph setup for circuit representation. + * This forms the hypergraph as a bipartite graph. Doesn't allow + * removing. Can only grab sets of connections to nodes + * + * @todo should replace with something better and more efficent. + * Should replace with a libraries setup instead. This would allow + * fast and easy partitioning of circuits + * + * @todo This is to replace inserting vector size for allocating PowerElectronicsModel + * + * @todo should replace N and E with Node and Component classes respectively. + * + * @note Tested but currently not used in the rest of the code. + * + * @tparam IdxT + * @tparam Label + */ +template +class CircuitGraph +{ +private: + std::set hypernodes; + std::set hyperedges; + std::map> edgestonodes; + +public: + CircuitGraph(); + ~CircuitGraph(); + bool addHyperEdge(E he); + bool addHyperNode(N hn); + bool addConnection(N hn, E he); + std::set getHyperEdgeConnections(E he); + size_t amountHyperNodes(); + size_t amountHyperEdges(); + void printBiPartiteGraph(bool verbose = false); +}; + +template +CircuitGraph::CircuitGraph() +{ +} + +template +CircuitGraph::~CircuitGraph() +{ +} + +template +bool CircuitGraph::addHyperNode(N hn) +{ + return this->hypernodes.insert(hn).second; +} + +template +bool CircuitGraph::addHyperEdge(E he) +{ + return this->hyperedges.insert(he).second; +} + +template +bool CircuitGraph::addConnection(N hn, E he) +{ + if (this->hyperedges.count(he) == 0 || this->hypernodes.count(hn) == 0) + { + return false; + } + return this->edgestonodes[he].insert(hn).second; +} + +template +std::set CircuitGraph::getHyperEdgeConnections(E he) +{ + return this->edgestonodes[he]; +} + +template +size_t CircuitGraph::amountHyperNodes() +{ + return this->hypernodes.size(); +} + +template +size_t CircuitGraph::amountHyperEdges() +{ + return this->hyperedges.size(); +} + +/** + * @brief Printing + * + * @todo need to add verbose printing for connections display + * + * @tparam IdxT + * @param verbose + */ + +template +void CircuitGraph::printBiPartiteGraph(bool verbose) +{ + + std::cout << "Amount of HyperNodes: " << this->amountHyperNodes() << std::endl; + std::cout << "Amount of HyperEdges: " << this->amountHyperEdges() << std::endl; + std::cout << "Connections per Edge:" << std::endl; + for (auto i : this->edgestonodes) + { + std::cout << i.first << " : {"; + for (auto j : i.second) + { + std::cout << j << ", "; + } + std::cout << "}\n"; + } +} diff --git a/ComponentLib/Bus/BusPV.cpp b/ComponentLib/Bus/BusPV.cpp index 5daafa10..d5a05dc9 100644 --- a/ComponentLib/Bus/BusPV.cpp +++ b/ComponentLib/Bus/BusPV.cpp @@ -74,7 +74,7 @@ namespace ModelLib { */ template BusPV::BusPV() - : BaseBus(0), V_(0.0), theta0_(0.0), Pg_(0.0) + : BaseBus(0), V_(0.0), theta0_(0.0) { //std::cout << "Create BusPV..." << std::endl; //std::cout << "Number of equations is " << size_ << std::endl; @@ -83,9 +83,7 @@ BusPV::BusPV() } /*! - * @brief BusPV constructor. - * - * This constructor sets initial values for voltage and phase angle. + * @brief Constructor for a PV bus * * @todo Arguments that should be passed to ModelEvaluatorImpl constructor: * - Number of equations = 1 (size_) @@ -94,10 +92,10 @@ BusPV::BusPV() * - Number of optimization parameters = 0 */ template -BusPV::BusPV(ScalarT V, ScalarT theta0, ScalarT Pg) - : BaseBus(0), V_(V), theta0_(theta0), Pg_(Pg) +BusPV::BusPV(ScalarT V, ScalarT theta0) + : BaseBus(0), V_(V), theta0_(theta0) { - //std::cout << "Create BusPV ..." << std::endl; + //std::cout << "Create BusPV..." << std::endl; //std::cout << "Number of equations is " << size_ << std::endl; size_ = 1; @@ -171,8 +169,8 @@ template int BusPV::evaluateResidual() { // std::cout << "Evaluating residual of a PV bus ...\n"; - P() = Pg_; - Q() = 0.0; + P() = 0.0; // <-- Residual P + Q() = 0.0; // <-- Output Qg, the reactive power generator needs to supply return 0; } diff --git a/ComponentLib/Bus/BusPV.hpp b/ComponentLib/Bus/BusPV.hpp index fd0291b6..cf645d59 100644 --- a/ComponentLib/Bus/BusPV.hpp +++ b/ComponentLib/Bus/BusPV.hpp @@ -91,7 +91,7 @@ namespace ModelLib using BusData = GridKit::PowerSystemData::BusData; BusPV(); - BusPV(ScalarT V, ScalarT theta0, ScalarT P); + BusPV(ScalarT V, ScalarT theta0); BusPV(BusData& data); virtual ~BusPV(); @@ -197,10 +197,9 @@ namespace ModelLib } private: - ScalarT V_; - ScalarT theta0_; ///< Default initial value for phase - ScalarT Pg_; ///< Generator injection - ScalarT Q_; + ScalarT V_; ///< Bus voltage magnitude + ScalarT theta0_; ///< Default initial value for phase + ScalarT Q_; ///< Reactive power that generator needs to provide ScalarT VB_; ScalarT thetaB_; diff --git a/ComponentLib/Bus/BusSlack.cpp b/ComponentLib/Bus/BusSlack.cpp index 80a5ae09..54318de7 100644 --- a/ComponentLib/Bus/BusSlack.cpp +++ b/ComponentLib/Bus/BusSlack.cpp @@ -97,6 +97,8 @@ BusSlack::BusSlack(ScalarT V, ScalarT theta) { //std::cout << "Create BusSlack..." << std::endl; //std::cout << "Number of equations is " << size_ << std::endl; + P() = 0.0; + Q() = 0.0; size_ = 0; } @@ -106,6 +108,8 @@ BusSlack::BusSlack(BusData& data) { //std::cout << "Create BusSlack..." << std::endl; //std::cout << "Number of equations is " << size_ << std::endl; + P() = 0.0; + Q() = 0.0; size_ = 0; } diff --git a/ComponentLib/CMakeLists.txt b/ComponentLib/CMakeLists.txt index 0c422d5b..29a5a798 100644 --- a/ComponentLib/CMakeLists.txt +++ b/ComponentLib/CMakeLists.txt @@ -69,3 +69,5 @@ add_subdirectory(Generator4Governor) add_subdirectory(Generator4Param) add_subdirectory(Load) add_subdirectory(MiniGrid) +add_subdirectory(PowerElectronicsComponents) + diff --git a/ComponentLib/Generator/GeneratorPV.hpp b/ComponentLib/Generator/GeneratorPV.hpp index 5e3be204..c34df84a 100644 --- a/ComponentLib/Generator/GeneratorPV.hpp +++ b/ComponentLib/Generator/GeneratorPV.hpp @@ -124,14 +124,18 @@ namespace ModelLib return P_; } + /// @brief Reactive power excess on PV bus + /// @return reference to negative PV generator reactive power virtual ScalarT& Q() { - return bus_->Q(); + return (bus_->Q()); } + /// @brief Reactive power excess on PV bus + /// @return const reference to negative PV generator reactive power virtual const ScalarT& Q() const { - return bus_->Q(); + return (bus_->Q()); } private: diff --git a/ComponentLib/PowerElectronicsComponents/CMakeLists.txt b/ComponentLib/PowerElectronicsComponents/CMakeLists.txt new file mode 100644 index 00000000..2cfe97bc --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/CMakeLists.txt @@ -0,0 +1,14 @@ + + +add_subdirectory(Capacitor) +add_subdirectory(Resistor) +add_subdirectory(VoltageSource) +add_subdirectory(Inductor) +add_subdirectory(LinearTransformer) +add_subdirectory(InductionMotor) +add_subdirectory(SynchronousMachine) +add_subdirectory(DistributedGenerator) +add_subdirectory(TransmissionLine) +add_subdirectory(MicrogridLoad) +add_subdirectory(MicrogridLine) +add_subdirectory(MicrogridBusDQ) \ No newline at end of file diff --git a/ComponentLib/PowerElectronicsComponents/Capacitor/CMakeLists.txt b/ComponentLib/PowerElectronicsComponents/Capacitor/CMakeLists.txt new file mode 100644 index 00000000..62efc179 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/Capacitor/CMakeLists.txt @@ -0,0 +1,8 @@ + + + +gridkit_add_library(powerelec_capacitor + SOURCES + Capacitor.cpp + OUTPUT_NAME + gridkit_powerelec_capacitor) \ No newline at end of file diff --git a/ComponentLib/PowerElectronicsComponents/Capacitor/Capacitor.cpp b/ComponentLib/PowerElectronicsComponents/Capacitor/Capacitor.cpp new file mode 100644 index 00000000..779b8a62 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/Capacitor/Capacitor.cpp @@ -0,0 +1,145 @@ + + + +#include +#include +#include +#include "Capacitor.hpp" + +namespace ModelLib { + +/*! + * @brief Constructor for Capacitor + * + * @todo this needs to be tested on some circuit + * + * Calls default ModelEvaluatorImpl constructor. + */ + +template +Capacitor::Capacitor(IdxT id, ScalarT C) + : C_(C) +{ + size_ = 3; + n_intern_ = 1; + n_extern_ = 2; + extern_indices_ = {0,1}; + idc_ = id; +} + +template +Capacitor::~Capacitor() +{ +} + +/*! + * @brief allocate method creates memory for vectors + */ +template +int Capacitor::allocate() +{ + y_.resize(size_); + yp_.resize(size_); + f_.resize(size_); + + return 0; +} + +/** + * Initialization of the grid model + */ +template +int Capacitor::initialize() +{ + return 0; +} + +/* + * \brief Identify differential variables + */ +template +int Capacitor::tagDifferentiable() +{ + return 0; +} + +/** + * @brief Evaluate the resisdual of the Capcitor + * + */ +template +int Capacitor::evaluateResidual() +{ + //input + f_[0] = C_ * yp_[2]; + //output + f_[1] = -C_ * yp_[2]; + + //internal + f_[2] = -C_ * yp_[2] + y_[0] - y_[1] - y_[2]; + return 0; +} + +/** + * @brief Compute the Jacobian dF/dy - a dF/dy' + * + * @tparam ScalarT + * @tparam IdxT + * @return int + */ +template +int Capacitor::evaluateJacobian() +{ + J_.zeroMatrix(); + //Create dF/dy + std::vector rcord{2,2,2}; + std::vector ccord{0,1,2}; + std::vector vals{1.0, -1.0, -1.0}; + J_.setValues(rcord, ccord, vals); + + //Create dF/dy' + std::vector rcordder{0,1,2}; + std::vector ccordder{2,2,2}; + std::vector valsder{C_, -C_, -C_}; + COO_Matrix Jacder = COO_Matrix(rcordder, ccordder, valsder,3,3); + + //Perform dF/dy + \alpha dF/dy' + J_.axpy(alpha_, Jacder); + + return 0; +} + +template +int Capacitor::evaluateIntegrand() +{ + return 0; +} + +template +int Capacitor::initializeAdjoint() +{ + return 0; +} + +template +int Capacitor::evaluateAdjointResidual() +{ + return 0; +} + +template +int Capacitor::evaluateAdjointIntegrand() +{ + return 0; +} + + + + +// Available template instantiations +template class Capacitor; +template class Capacitor; + + +} //namespace ModelLib + diff --git a/ComponentLib/PowerElectronicsComponents/Capacitor/Capacitor.hpp b/ComponentLib/PowerElectronicsComponents/Capacitor/Capacitor.hpp new file mode 100644 index 00000000..5b5fb713 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/Capacitor/Capacitor.hpp @@ -0,0 +1,68 @@ + + +#ifndef _CAP_HPP_ +#define _CAP_HPP_ + +#include +#include +#include + + +namespace ModelLib +{ + template class BaseBus; +} + + +namespace ModelLib +{ + /*! + * @brief Declaration of a Capacitor class. + * + */ + template + class Capacitor : public CircuitComponent + { + using CircuitComponent::size_; + using CircuitComponent::nnz_; + using CircuitComponent::time_; + using CircuitComponent::alpha_; + using CircuitComponent::y_; + using CircuitComponent::yp_; + using CircuitComponent::tag_; + using CircuitComponent::f_; + using CircuitComponent::g_; + using CircuitComponent::yB_; + using CircuitComponent::ypB_; + using CircuitComponent::fB_; + using CircuitComponent::gB_; + using CircuitComponent::J_; + using CircuitComponent::param_; + using CircuitComponent::idc_; + + using CircuitComponent::extern_indices_; + using CircuitComponent::n_extern_; + using CircuitComponent::n_intern_; + + public: + Capacitor(IdxT id, ScalarT C); + virtual ~Capacitor(); + + int allocate(); + int initialize(); + int tagDifferentiable(); + int evaluateResidual(); + int evaluateJacobian(); + int evaluateIntegrand(); + + int initializeAdjoint(); + int evaluateAdjointResidual(); + //int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); + + private: + ScalarT C_; + }; +} + +#endif diff --git a/ComponentLib/PowerElectronicsComponents/CircuitComponent.hpp b/ComponentLib/PowerElectronicsComponents/CircuitComponent.hpp new file mode 100644 index 00000000..8186c727 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/CircuitComponent.hpp @@ -0,0 +1,86 @@ + + +#ifndef _CIRCCOMP_HPP_ +#define _CIRCCOMP_HPP_ + +#include +#include +#include +#include + + +namespace ModelLib +{ + /*! + * @brief Declaration of a CircuitComponent class. + * + */ + template + class CircuitComponent : public ModelEvaluatorImpl + { + + public: + + + void updateTime(ScalarT t, ScalarT a) + { + this->time_ = t; + this->alpha_ = a; + } + + bool hasJacobian() { return true;} + + size_t getExternSize() + { + return n_extern_; + } + + size_t getInternalSize() + { + return this->n_intern_; + } + + std::set getExternIndices() + { + return this->extern_indices_; + } + + /** + * @brief Create the mappings from local to global indexes + * + * @param local_index + * @param global_index + * @return int + */ + int setExternalConnectionNodes(IdxT local_index, IdxT global_index) + { + connection_nodes_[local_index] = global_index; + return 0; + } + + /** + * @brief Given the location of value in the local vector map to global index + * + * f(local_index) = global_index + * + * @param local_index index of local value in vector + * @return IdxT Index of the same value in the global vector + */ + IdxT getNodeConnection(IdxT local_index) + { + return connection_nodes_.at(local_index); + } + + protected: + size_t n_extern_; + size_t n_intern_; + std::set extern_indices_; + ///@todo may want to replace the mapping of connection_nodes to Node objects instead of IdxT. Allows for container free setup + std::map connection_nodes_; + + }; + + +} + +#endif diff --git a/ComponentLib/PowerElectronicsComponents/DistributedGenerator/CMakeLists.txt b/ComponentLib/PowerElectronicsComponents/DistributedGenerator/CMakeLists.txt new file mode 100644 index 00000000..22652f53 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/DistributedGenerator/CMakeLists.txt @@ -0,0 +1,8 @@ + + + +gridkit_add_library(powerelec_disgen + SOURCES + DistributedGenerator.cpp + OUTPUT_NAME + gridkit_powerelec_disgen) \ No newline at end of file diff --git a/ComponentLib/PowerElectronicsComponents/DistributedGenerator/DistributedGenerator.cpp b/ComponentLib/PowerElectronicsComponents/DistributedGenerator/DistributedGenerator.cpp new file mode 100644 index 00000000..21b91d75 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/DistributedGenerator/DistributedGenerator.cpp @@ -0,0 +1,363 @@ + + + +#include +#include +#include +#include "DistributedGenerator.hpp" + +namespace ModelLib { + + +/*! + * @brief Constructor for a Discrete Generator + * @todo Maybe have parameters be templated in. Variables cannot be changed and are unlikely to. Allows for compile time optimizations + * + * Calls default ModelEvaluatorImpl constructor. + */ + +template +DistributedGenerator::DistributedGenerator(IdxT id, DistributedGeneratorParameters parm, bool reference_frame) + : wb_(parm.wb_), + wc_(parm.wc_), + mp_(parm.mp_), + Vn_(parm.Vn_), + nq_(parm.nq_), + F_(parm.F_), + Kiv_(parm.Kiv_), + Kpv_(parm.Kpv_), + Kic_(parm.Kic_), + Kpc_(parm.Kpc_), + Cf_(parm.Cf_), + rLf_(parm.rLf_), + Lf_(parm.Lf_), + rLc_(parm.rLc_), + Lc_(parm.Lc_), + refframe_(reference_frame) +{ + // 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; + extern_indices_ = {0,1,2}; + idc_ = id; +} + +template +DistributedGenerator::~DistributedGenerator() +{ +} + +/*! + * @brief allocate method computes sparsity pattern of the Jacobian. + */ +template +int DistributedGenerator::allocate() +{ + y_.resize(size_); + yp_.resize(size_); + f_.resize(size_); + + return 0; +} + +/** + * Initialization of the grid model + */ +template +int DistributedGenerator::initialize() +{ + return 0; +} + +/* + * \brief Identify differential variables + */ +template +int DistributedGenerator::tagDifferentiable() +{ + return 0; +} + +/** + * @brief Contributes to the resisdual of the Distributed Generator + * + */ +template +int DistributedGenerator::evaluateResidual() +{ + // ### Externals Componenets ### + + ScalarT omega = wb_ - mp_ * y_[4]; + //ref common ref motor angle + /// @todo fix boolian conditional, unclear result + if (refframe_) + { + f_[0] = omega - y_[0]; + } + else + { + f_[0] = 0.0; + } + + + //output + //current transformed to common frame + f_[1] = cos(y_[3]) * y_[14] - sin(y_[3]) * y_[15]; + f_[2] = sin(y_[3]) * y_[14] + cos(y_[3]) * y_[15]; + + //Take incoming voltages to current rotator reference frame + ScalarT vbd_in = cos(y_[3]) * y_[1] + sin(y_[3]) * y_[2]; + ScalarT vbq_in = -sin(y_[3]) * y_[1] + cos(y_[3]) * y_[2]; + + // ### Internal Componenets ## + // Rotor difference angle + f_[3] = -yp_[3] + omega - y_[0]; + + // P and Q equations + f_[4] = -yp_[4] + wc_ * (y_[12] * y_[14] + y_[13] * y_[15] - y_[4]); + f_[5] = -yp_[5] + wc_ * (-y_[12] * y_[15] + y_[13] * y_[14] - y_[5]); + + //Voltage control + ScalarT vod_star = Vn_ - nq_ * y_[5]; + ScalarT voq_star = 0.0; + + f_[6] = -yp_[6] + vod_star - y_[12]; + f_[7] = -yp_[7] + voq_star - y_[13]; + + + ScalarT ild_star = F_ * y_[14] - wb_ * Cf_ * y_[13] + Kpv_ * (vod_star - y_[12]) + Kiv_ * y_[6]; + ScalarT ilq_star = F_ * y_[15] + wb_ * Cf_ * y_[12] + Kpv_ * (voq_star - y_[13]) + Kiv_ * y_[7]; + + //Current control + f_[8] = -yp_[8] + ild_star - y_[10]; + f_[9] = -yp_[9] + ilq_star - y_[11]; + + ScalarT vid_star = -wb_ * Lf_ * y_[11] + Kpc_ * (ild_star - y_[10]) + Kic_ * y_[8]; + ScalarT viq_star = wb_ * Lf_ * y_[10] + Kpc_ * (ilq_star - y_[11]) + Kic_ * y_[9]; + + //Output LC Filter + f_[10] = -yp_[10] - (rLf_ / Lf_) * y_[10] + omega * y_[11] + (vid_star - y_[12]) / Lf_; + f_[11] = -yp_[11] - (rLf_ / Lf_) * y_[11] - omega * y_[10] + (viq_star - y_[13]) / Lf_; + + f_[12] = -yp_[12] + omega * y_[13] + (y_[10] - y_[14]) / Cf_; + f_[13] = -yp_[13] - omega * y_[12] + (y_[11] - y_[15]) / Cf_; + + //Output Connector + f_[14] = -yp_[14] - (rLc_ / Lc_) * y_[14] + omega * y_[15] + (y_[12] - vbd_in) / Lc_; + f_[15] = -yp_[15] - (rLc_ / Lc_) * y_[15] - omega * y_[14] + (y_[13] - vbq_in) / Lc_; + return 0; +} + +/** + * @brief Compute the jacobian of the DistributedGenerator for iteration. dF/dy - \alpha dF/dy' + * + * The matrix dF/dy should be + * +[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] +[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] +[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] +[-1, 0, 0, 0, -mp, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] +[ 0, 0, 0, 0, -wc, 0, 0, 0, 0, 0, 0, 0, wc*x15, wc*x16, wc*x13, wc*x14] +[ 0, 0, 0, 0, 0, -wc, 0, 0, 0, 0, 0, 0, -wc*x16, wc*x15, wc*x14, -wc*x13] +[ 0, 0, 0, 0, 0, -nq, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0] +[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0] +[ 0, 0, 0, 0, 0, -Kpv*nq, Kiv, 0, 0, 0, -1, 0, -Kpv, -Cf*wb, F, 0] +[ 0, 0, 0, 0, 0, 0, 0, Kiv, 0, 0, 0, -1, Cf*wb, -Kpv, 0, F] +[ 0, 0, 0, 0, -mp*x12, -(Kpc*Kpv*nq)/Lf, (Kiv*Kpc)/Lf, 0, Kic/Lf, 0, - Kpc/Lf - rLf/Lf, -mp*x5, -(Kpc*Kpv + 1)/Lf, -(Cf*Kpc*wb)/Lf, (F*Kpc)/Lf, 0] +[ 0, 0, 0, 0, mp*x11, 0, 0, (Kiv*Kpc)/Lf, 0, Kic/Lf, mp*x5, - Kpc/Lf - rLf/Lf, (Cf*Kpc*wb)/Lf, -(Kpc*Kpv + 1)/Lf, 0, (F*Kpc)/Lf] +[ 0, 0, 0, 0, -mp*x14, 0, 0, 0, 0, 0, 1/Cf, 0, 0, wb - mp*x5, -1/Cf, 0] +[ 0, 0, 0, 0, mp*x13, 0, 0, 0, 0, 0, 0, 1/Cf, mp*x5 - wb, 0, 0, -1/Cf] +[ 0, -cos(x4)/Lc, -sin(x4)/Lc, -(x3*cos(x4) - x2*sin(x4))/Lc, -mp*x16, 0, 0, 0, 0, 0, 0, 0, 1/Lc, 0, -rLc/Lc, wb - mp*x5] +[ 0, sin(x4)/Lc, -cos(x4)/Lc, (x2*cos(x4) + x3*sin(x4))/Lc, mp*x15, 0, 0, 0, 0, 0, 0, 0, 0, 1/Lc, mp*x5 - wb, -rLc/Lc] + * 'Generated from MATLAB symbolic' + * + * @tparam ScalarT + * @tparam IdxT + * @return int + */ +template +int DistributedGenerator::evaluateJacobian() +{ + J_.zeroMatrix(); + //Create dF/dy' + std::vector rcordder(13); + std::vector valsder(13, -1.0); + for (int i = 0; i < 13; i++) + { + rcordder[i] = i + 3; + } + COO_Matrix Jacder = COO_Matrix(rcordder, rcordder, valsder,16,16); + + std::vector ctemp{}; + std::vector rtemp{}; + std::vector valtemp{}; + + + + //Create dF/dy + //r = 1 + + ctemp = {3, 14, 15}; + rtemp.clear(); + for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(1); + valtemp = { - sin(y_[3]) * y_[14] - cos(y_[3]) * y_[15], cos(y_[3]),-sin(y_[3])}; + J_.setValues(rtemp, ctemp, valtemp); + + //r = 2 + + ctemp = {3, 14, 15}; + rtemp.clear(); + for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(2); + valtemp = { cos(y_[3]) * y_[14] - sin(y_[3]) * y_[15], sin(y_[3]),cos(y_[3])}; + J_.setValues(rtemp, ctemp, valtemp); + + //r = 3 + + ctemp = {0, 4}; + rtemp.clear(); + for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(3); + valtemp = {-1.0, -mp_}; + J_.setValues(rtemp, ctemp, valtemp); + + //r = 0 + if (refframe_) + { + ctemp = {0, 4}; + rtemp.clear(); + for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(0); + valtemp = {-1.0, -mp_}; + J_.setValues(rtemp, ctemp, valtemp); + } + + + //r = 4 + ctemp = {4, 12, 13, 14, 15}; + rtemp.clear(); + for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(4); + valtemp = {-wc_, wc_*y_[14], wc_*y_[15], wc_*y_[12], wc_*y_[13]}; + J_.setValues(rtemp, ctemp, valtemp); + + //r = 5 + ctemp = {5, 12, 13, 14, 15}; + rtemp.clear(); + for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(5); + valtemp = {-wc_, -wc_*y_[15], wc_*y_[14], wc_*y_[13], -wc_*y_[12]}; + J_.setValues(rtemp, ctemp, valtemp); + + //r = 6 + ctemp = {5, 12}; + rtemp.clear(); + for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(6); + valtemp = {-nq_, -1.0}; + J_.setValues(rtemp, ctemp, valtemp); + + //r = 7 + ctemp = {13}; + rtemp.clear(); + for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(7); + valtemp = {-1.0}; + J_.setValues(rtemp, ctemp, valtemp); + + //r = 8 + ctemp = {5,6,10,12,13,14}; + rtemp.clear(); + for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(8); + valtemp = {-Kpv_*nq_, Kiv_, -1.0, -Kpv_, -Cf_*wb_, F_}; + J_.setValues(rtemp, ctemp, valtemp); + + //r = 9 + ctemp = {7, 11, 12, 13, 15}; + rtemp.clear(); + for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(9); + valtemp = {Kiv_, -1.0, Cf_*wb_,-Kpv_,F_}; + J_.setValues(rtemp, ctemp, valtemp); + + //r = 10 + ctemp = {4, 5, 6, 8, 10, 11, 12, 13, 14}; + rtemp.clear(); + for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(10); + valtemp = {-mp_ * y_[11], -(Kpc_ * Kpv_ * nq_) / Lf_, (Kpc_ * Kiv_) / Lf_, Kic_ / Lf_, -(Kpc_ + rLf_) / Lf_, -mp_ * y_[4], -(Kpc_ * Kpv_ + 1.0) / Lf_, -(Cf_ * Kpc_ * wb_) / Lf_, (F_ * Kpc_) / Lf_}; + J_.setValues(rtemp, ctemp, valtemp); + + //r = 11 + ctemp = {4, 7, 9, 10, 11, 12, 13, 15}; + rtemp.clear(); + for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(11); + valtemp = {mp_ * y_[10], (Kiv_ * Kpc_) / Lf_, Kic_ / Lf_, mp_ * y_[4], -(Kpc_ + rLf_) / Lf_, (Cf_ * Kpc_ * wb_) / Lf_, -(Kpc_ * Kpv_ + 1.0) / Lf_, (F_ * Kpc_) / Lf_}; + J_.setValues(rtemp, ctemp, valtemp); + + //r = 12 + ctemp = {4, 10, 13, 14}; + rtemp.clear(); + for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(12); + valtemp = {-mp_ * y_[13], 1.0 / Cf_, wb_ - mp_ * y_[4], -1.0 / Cf_}; + J_.setValues(rtemp, ctemp, valtemp); + + + //r = 13 + ctemp = {4, 11, 12, 15}; + rtemp.clear(); + for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(13); + valtemp = {mp_ * y_[12], 1.0 / Cf_, -wb_ + mp_ * y_[4], -1.0 / Cf_}; + J_.setValues(rtemp, ctemp, valtemp); + + + //r = 14 + ctemp = {1, 2, 3, 4, 12, 14, 15}; + rtemp.clear(); + for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(14); + valtemp = {(1.0/Lc_) * -cos(y_[3]) , (1.0/Lc_) * -sin(y_[3]) , (1.0/Lc_) * (sin(y_[3]) * y_[1] - cos(y_[3]) * y_[2]), -mp_ * y_[15], 1.0 / Lc_, -rLc_ / Lc_, wb_ - mp_ * y_[4]}; + J_.setValues(rtemp, ctemp, valtemp); + + + //r = 15 + ctemp = {1, 2, 3, 4, 13, 14, 15}; + rtemp.clear(); + for (size_t i = 0; i < ctemp.size(); i++) rtemp.push_back(15); + valtemp = {(1.0/Lc_) * sin(y_[3]) , (1.0/Lc_) * -cos(y_[3]), (1.0/Lc_) * (cos(y_[3]) * y_[1] + sin(y_[3]) * y_[2]), mp_ * y_[14], 1.0 / Lc_, -wb_ + mp_ * y_[4], -rLc_ / Lc_}; + J_.setValues(rtemp, ctemp, valtemp); + + + //Perform dF/dy + \alpha dF/dy' + + J_.axpy(alpha_, Jacder); + + return 0; +} + +template +int DistributedGenerator::evaluateIntegrand() +{ + return 0; +} + +template +int DistributedGenerator::initializeAdjoint() +{ + return 0; +} + +template +int DistributedGenerator::evaluateAdjointResidual() +{ + return 0; +} + +template +int DistributedGenerator::evaluateAdjointIntegrand() +{ + return 0; +} + + + + +// Available template instantiations +template class DistributedGenerator; +template class DistributedGenerator; + + +} //namespace ModelLib + diff --git a/ComponentLib/PowerElectronicsComponents/DistributedGenerator/DistributedGenerator.hpp b/ComponentLib/PowerElectronicsComponents/DistributedGenerator/DistributedGenerator.hpp new file mode 100644 index 00000000..9946f9a2 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/DistributedGenerator/DistributedGenerator.hpp @@ -0,0 +1,103 @@ + + +#ifndef _CAP_HPP_ +#define _CAP_HPP_ + +#include +#include +#include + + +namespace ModelLib +{ + template class BaseBus; + + template + struct DistributedGeneratorParameters + { + ScalarT wb_; + ScalarT wc_; + ScalarT mp_; + ScalarT Vn_; + ScalarT nq_; + ScalarT F_; + ScalarT Kiv_; + ScalarT Kpv_; + ScalarT Kic_; + ScalarT Kpc_; + ScalarT Cf_; + ScalarT rLf_; + ScalarT Lf_; + ScalarT rLc_; + ScalarT Lc_; + }; +} + + +namespace ModelLib +{ + /*! + * @brief Declaration of a DistributedGenerator class. + * + */ + template + class DistributedGenerator : public CircuitComponent + { + using CircuitComponent::size_; + using CircuitComponent::nnz_; + using CircuitComponent::time_; + using CircuitComponent::alpha_; + using CircuitComponent::y_; + using CircuitComponent::yp_; + using CircuitComponent::tag_; + using CircuitComponent::f_; + using CircuitComponent::g_; + using CircuitComponent::yB_; + using CircuitComponent::ypB_; + using CircuitComponent::fB_; + using CircuitComponent::gB_; + using CircuitComponent::J_; + using CircuitComponent::param_; + using CircuitComponent::idc_; + + using CircuitComponent::extern_indices_; + using CircuitComponent::n_extern_; + using CircuitComponent::n_intern_; + + + public: + DistributedGenerator(IdxT id, DistributedGeneratorParameters parm, bool reference_frame); + virtual ~DistributedGenerator(); + + int allocate(); + int initialize(); + int tagDifferentiable(); + int evaluateResidual(); + int evaluateJacobian(); + int evaluateIntegrand(); + int initializeAdjoint(); + int evaluateAdjointResidual(); + //int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); + + private: + ScalarT wb_; + ScalarT wc_; + ScalarT mp_; + ScalarT Vn_; + ScalarT nq_; + ScalarT F_; + ScalarT Kiv_; + ScalarT Kpv_; + ScalarT Kic_; + ScalarT Kpc_; + ScalarT Cf_; + ScalarT rLf_; + ScalarT Lf_; + ScalarT rLc_; + ScalarT Lc_; + bool refframe_; + }; +} + +#endif diff --git a/ComponentLib/PowerElectronicsComponents/InductionMotor/CMakeLists.txt b/ComponentLib/PowerElectronicsComponents/InductionMotor/CMakeLists.txt new file mode 100644 index 00000000..95ef4243 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/InductionMotor/CMakeLists.txt @@ -0,0 +1,8 @@ + + + +gridkit_add_library(powerelec_inductionmotor + SOURCES + InductionMotor.cpp + OUTPUT_NAME + gridkit_powerelec_inductionmotor) \ No newline at end of file diff --git a/ComponentLib/PowerElectronicsComponents/InductionMotor/InductionMotor.cpp b/ComponentLib/PowerElectronicsComponents/InductionMotor/InductionMotor.cpp new file mode 100644 index 00000000..fa0da7fb --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/InductionMotor/InductionMotor.cpp @@ -0,0 +1,145 @@ + + + +#include +#include +#include +#include "InductionMotor.hpp" + + +namespace ModelLib { + + + +/*! + * @brief Constructor for a constant InductionMotor model + * + * Calls default ModelEvaluatorImpl constructor. + * @todo create a test case utilizing the component. + * @todo create a unit test to check correctness of component + */ + +template +InductionMotor::InductionMotor(IdxT id, ScalarT Lls, ScalarT Rs, ScalarT Llr, ScalarT Rr, ScalarT Lms, ScalarT RJ, ScalarT P) + : Lls_(Lls), + Rs_(Rs), + Llr_(Llr), + Rr_(Rr), + Lms_(Lms), + RJ_(RJ), + P_(P) +{ + size_ = 10; + n_intern_ = 5; + n_extern_ = 5; + extern_indices_ = {0,1,2,3,4}; + idc_ = id; +} + +template +InductionMotor::~InductionMotor() +{ +} + +/*! + * @brief allocate method computes sparsity pattern of the Jacobian. + */ +template +int InductionMotor::allocate() +{ + y_.resize(size_); + yp_.resize(size_); + f_.resize(size_); + return 0; +} + +/** + * Initialization of the grid model + */ +template +int InductionMotor::initialize() +{ + return 0; +} + +/* + * \brief Identify differential variables + */ +template +int InductionMotor::tagDifferentiable() +{ + return 0; +} + +/** + * @brief Contributes to the resisdual + * + */ +template +int InductionMotor::evaluateResidual() +{ + + f_[0] = y_[5] + y_[7]; + f_[1] = (-1.0/2.0) * y_[5] - (sqrt(3.0)/2.0)*y_[6] + y_[7]; + f_[2] = (-1.0/2.0) * y_[5] + (sqrt(3.0)/2.0)*y_[6] + y_[7]; + f_[3] = RJ_ * yp_[3] - (3.0/4.0)*P_*Lms_ * (y_[5]*y_[9] - y_[6]*y_[8]); + f_[4] = yp_[4] - y_[3]; + f_[5] = (1.0/3.0)*(2.0* y_[0] - y_[1] - y_[2]) - Rs_*y_[5] - (Lls_ + Lms_) * yp_[5] - Lms_ * yp_[6]; + f_[6] = (1.0/sqrt(3.0))*(-y_[1] + y_[2]) - Rs_*y_[6] - (Lls_ + Lms_) * yp_[6] - Lms_ * yp_[5]; + f_[7] = (y_[0] + y_[1] + y_[2])/3.0 - Rs_*y_[7] - Lls_ * yp_[7]; + f_[8] = Rr_*y_[8] + (Llr_ + Lms_)*yp_[8] + Lms_ * yp_[5] - (P_/2)*y_[3]*((Llr_+Lms_)*y_[9] + Lms_*y_[6]); + f_[9] = Rr_*y_[9] + (Llr_ + Lms_)*yp_[9] + Lms_ * yp_[6] + (P_/2)*y_[3]*((Llr_+Lms_)*y_[8] + Lms_*y_[5]); + return 0; +} + +/** + * @brief Compute component Jacobian + * + * @todo need to implement + * + * @tparam ScalarT + * @tparam IdxT + * @return int + */ +template +int InductionMotor::evaluateJacobian() +{ + + return 0; +} + +template +int InductionMotor::evaluateIntegrand() +{ + return 0; +} + +template +int InductionMotor::initializeAdjoint() +{ + return 0; +} + +template +int InductionMotor::evaluateAdjointResidual() +{ + return 0; +} + +template +int InductionMotor::evaluateAdjointIntegrand() +{ + return 0; +} + + + + + +// Available template instantiations +template class InductionMotor; +template class InductionMotor; + + +} //namespace ModelLib + diff --git a/ComponentLib/PowerElectronicsComponents/InductionMotor/InductionMotor.hpp b/ComponentLib/PowerElectronicsComponents/InductionMotor/InductionMotor.hpp new file mode 100644 index 00000000..cb49f761 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/InductionMotor/InductionMotor.hpp @@ -0,0 +1,74 @@ + + +#ifndef _IMOTOR_HPP_ +#define _IMOTOR_HPP_ + +#include +#include +#include + + +namespace ModelLib +{ + template class BaseBus; +} + + +namespace ModelLib +{ + /*! + * @brief Declaration of a InductionMotor class. + * + */ + template + class InductionMotor : public CircuitComponent + { + using CircuitComponent::size_; + using CircuitComponent::nnz_; + using CircuitComponent::time_; + using CircuitComponent::alpha_; + using CircuitComponent::y_; + using CircuitComponent::yp_; + using CircuitComponent::tag_; + using CircuitComponent::f_; + using CircuitComponent::g_; + using CircuitComponent::yB_; + using CircuitComponent::ypB_; + using CircuitComponent::fB_; + using CircuitComponent::gB_; + using CircuitComponent::J_; + using CircuitComponent::param_; + using CircuitComponent::idc_; + + using CircuitComponent::extern_indices_; + using CircuitComponent::n_extern_; + using CircuitComponent::n_intern_; + + public: + InductionMotor(IdxT id, ScalarT Lls, ScalarT Rs, ScalarT Llr, ScalarT Rr, ScalarT Lms, ScalarT RJ, ScalarT P); + virtual ~InductionMotor(); + + int allocate(); + int initialize(); + int tagDifferentiable(); + int evaluateResidual(); + int evaluateJacobian(); + int evaluateIntegrand(); + + int initializeAdjoint(); + int evaluateAdjointResidual(); + //int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); + + private: + ScalarT Lls_; + ScalarT Rs_; + ScalarT Llr_; + ScalarT Rr_; + ScalarT Lms_; + ScalarT RJ_; + ScalarT P_; + }; +} + +#endif diff --git a/ComponentLib/PowerElectronicsComponents/Inductor/CMakeLists.txt b/ComponentLib/PowerElectronicsComponents/Inductor/CMakeLists.txt new file mode 100644 index 00000000..4c620030 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/Inductor/CMakeLists.txt @@ -0,0 +1,8 @@ + + + +gridkit_add_library(powerelec_inductor + SOURCES + Inductor.cpp + OUTPUT_NAME + gridkit_powerelec_inductor) \ No newline at end of file diff --git a/ComponentLib/PowerElectronicsComponents/Inductor/Inductor.cpp b/ComponentLib/PowerElectronicsComponents/Inductor/Inductor.cpp new file mode 100644 index 00000000..2bb4dd1d --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/Inductor/Inductor.cpp @@ -0,0 +1,144 @@ + + + +#include +#include +#include +#include "Inductor.hpp" + +namespace ModelLib { + +/*! + * @brief Constructor for a inductor + * + * Calls default ModelEvaluatorImpl constructor. + */ + +template +Inductor::Inductor(IdxT id, ScalarT L) + : L_(L) +{ + size_ = 3; + n_intern_ = 1; + n_extern_ = 2; + extern_indices_ = {0,1}; + idc_ = id; +} + +template +Inductor::~Inductor() +{ +} + +/*! + * @brief allocate method computes sparsity pattern of the Jacobian. + */ +template +int Inductor::allocate() +{ + + y_.resize(size_); + yp_.resize(size_); + f_.resize(size_); + + return 0; +} + +/** + * Initialization of the grid model + */ +template +int Inductor::initialize() +{ + return 0; +} + +/* + * \brief Identify differential variables + */ +template +int Inductor::tagDifferentiable() +{ + return 0; +} + +/** + * @brief Compute the resisdual of the component + * + */ +template +int Inductor::evaluateResidual() +{ + //input + f_[0] = -y_[2]; + //output + f_[1] = y_[2]; + //internal + f_[2] = -L_ * yp_[2] + y_[1] - y_[0] ; + return 0; +} + +/** + * @brief Evaluate the jacobian of the component + * + * @tparam ScalarT + * @tparam IdxT + * @return int + */ +template +int Inductor::evaluateJacobian() +{ + J_.zeroMatrix(); + + //Create dF/dy + std::vector rcord{0,1,2,2}; + std::vector ccord{2,2,0,1}; + std::vector vals{-1.0, 1.0, -1.0, 1.0}; + J_.setValues(rcord, ccord, vals); + + //Create dF/dy' + std::vector rcordder{2}; + std::vector ccordder{2}; + std::vector valsder{-L_}; + COO_Matrix Jacder = COO_Matrix(rcordder, ccordder, valsder,3,3); + + //Perform dF/dy + \alpha dF/dy' + J_.axpy(alpha_, Jacder); + + return 0; +} + +template +int Inductor::evaluateIntegrand() +{ + return 0; +} + +template +int Inductor::initializeAdjoint() +{ + return 0; +} + +template +int Inductor::evaluateAdjointResidual() +{ + return 0; +} + +template +int Inductor::evaluateAdjointIntegrand() +{ + return 0; +} + + + + +// Available template instantiations +template class Inductor; +template class Inductor; + + +} //namespace ModelLib + diff --git a/ComponentLib/PowerElectronicsComponents/Inductor/Inductor.hpp b/ComponentLib/PowerElectronicsComponents/Inductor/Inductor.hpp new file mode 100644 index 00000000..2b0f8ec6 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/Inductor/Inductor.hpp @@ -0,0 +1,72 @@ + + +#ifndef _IND_HPP_ +#define _IND_HPP_ + +#include +#include +#include + + +namespace ModelLib +{ + template class BaseBus; +} + + +namespace ModelLib +{ + /*! + * @brief Declaration of a Inductor class. + * + */ + template + class Inductor : public CircuitComponent + { + using CircuitComponent::size_; + using CircuitComponent::nnz_; + using CircuitComponent::time_; + using CircuitComponent::alpha_; + using CircuitComponent::y_; + using CircuitComponent::yp_; + using CircuitComponent::tag_; + using CircuitComponent::f_; + using CircuitComponent::g_; + using CircuitComponent::yB_; + using CircuitComponent::ypB_; + using CircuitComponent::fB_; + using CircuitComponent::gB_; + using CircuitComponent::J_; + using CircuitComponent::param_; + using CircuitComponent::idc_; + + using CircuitComponent::extern_indices_; + using CircuitComponent::n_extern_; + using CircuitComponent::n_intern_; + + + + + public: + Inductor(IdxT id, ScalarT L); + virtual ~Inductor(); + + int allocate(); + int initialize(); + int tagDifferentiable(); + int evaluateResidual(); + int evaluateJacobian(); + int evaluateIntegrand(); + + int initializeAdjoint(); + int evaluateAdjointResidual(); + //int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); + + + private: + ScalarT L_; + }; +} + +#endif diff --git a/ComponentLib/PowerElectronicsComponents/LinearTransformer/CMakeLists.txt b/ComponentLib/PowerElectronicsComponents/LinearTransformer/CMakeLists.txt new file mode 100644 index 00000000..eaf95c04 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/LinearTransformer/CMakeLists.txt @@ -0,0 +1,8 @@ + + + +gridkit_add_library(powerelec_lineartrasnformer + SOURCES + LinearTransformer.cpp + OUTPUT_NAME + gridkit_powerelec_lineartrasnformer) \ No newline at end of file diff --git a/ComponentLib/PowerElectronicsComponents/LinearTransformer/LinearTransformer.cpp b/ComponentLib/PowerElectronicsComponents/LinearTransformer/LinearTransformer.cpp new file mode 100644 index 00000000..b53fbd27 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/LinearTransformer/LinearTransformer.cpp @@ -0,0 +1,123 @@ + + + +#include +#include +#include +#include "LinearTransformer.hpp" + +namespace ModelLib { + +/*! + * @brief Constructor for a LinearTransformer model + * + * Calls default ModelEvaluatorImpl constructor. + * @todo Not tested in any model yet. Should be + * @todo Has not been tested for correctness + */ + +template +LinearTransformer::LinearTransformer(IdxT id, ScalarT L1, ScalarT L2, ScalarT R1, ScalarT R2, ScalarT M) + : L1_(L1), + L2_(L2), + R1_(R1), + R2_(R2), + M_(M) +{ + size_ = 4; + n_intern_ = 2; + n_extern_ = 2; + extern_indices_ = {0,1}; + idc_ = id; +} + +template +LinearTransformer::~LinearTransformer() +{ +} + +/*! + * @brief allocate method computes sparsity pattern of the Jacobian. + */ +template +int LinearTransformer::allocate() +{ + y_.resize(size_); + yp_.resize(size_); + f_.resize(size_); + return 0; +} + +/** + * Initialization of the grid model + */ +template +int LinearTransformer::initialize() +{ + return 0; +} + +/* + * \brief Identify differential variables + */ +template +int LinearTransformer::tagDifferentiable() +{ + return 0; +} + +/** + * @brief Computes the component resisdual + * + */ +template +int LinearTransformer::evaluateResidual() +{ + f_[0] = y_[2]; + f_[1] = y_[3]; + f_[2] = y_[0] - R1_ * y_[2] - L1_ * yp_[2] - M_ * yp_[3]; + f_[2] = y_[1] - R2_ * y_[3] - M_ * yp_[2] - L2_ * yp_[3]; + return 0; +} + +template +int LinearTransformer::evaluateJacobian() +{ + return 0; +} + +template +int LinearTransformer::evaluateIntegrand() +{ + return 0; +} + +template +int LinearTransformer::initializeAdjoint() +{ + return 0; +} + +template +int LinearTransformer::evaluateAdjointResidual() +{ + return 0; +} + +template +int LinearTransformer::evaluateAdjointIntegrand() +{ + return 0; +} + + + + + +// Available template instantiations +template class LinearTransformer; +template class LinearTransformer; + + +} //namespace ModelLib + diff --git a/ComponentLib/PowerElectronicsComponents/LinearTransformer/LinearTransformer.hpp b/ComponentLib/PowerElectronicsComponents/LinearTransformer/LinearTransformer.hpp new file mode 100644 index 00000000..6d9a920c --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/LinearTransformer/LinearTransformer.hpp @@ -0,0 +1,72 @@ + + +#ifndef _LTRANS_HPP_ +#define _LTRANS_HPP_ + +#include +#include +#include + + +namespace ModelLib +{ + template class BaseBus; +} + + +namespace ModelLib +{ + /*! + * @brief Declaration of a LinearTransformer class. + * + */ + template + class LinearTransformer : public CircuitComponent + { + using CircuitComponent::size_; + using CircuitComponent::nnz_; + using CircuitComponent::time_; + using CircuitComponent::alpha_; + using CircuitComponent::y_; + using CircuitComponent::yp_; + using CircuitComponent::tag_; + using CircuitComponent::f_; + using CircuitComponent::g_; + using CircuitComponent::yB_; + using CircuitComponent::ypB_; + using CircuitComponent::fB_; + using CircuitComponent::gB_; + using CircuitComponent::J_; + using CircuitComponent::param_; + using CircuitComponent::idc_; + + using CircuitComponent::extern_indices_; + using CircuitComponent::n_extern_; + using CircuitComponent::n_intern_; + + public: + LinearTransformer(IdxT id, ScalarT L1, ScalarT L2, ScalarT R1, ScalarT R2, ScalarT M); + virtual ~LinearTransformer(); + + int allocate(); + int initialize(); + int tagDifferentiable(); + int evaluateResidual(); + int evaluateJacobian(); + int evaluateIntegrand(); + + int initializeAdjoint(); + int evaluateAdjointResidual(); + //int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); + + private: + ScalarT L1_; + ScalarT L2_; + ScalarT R1_; + ScalarT R2_; + ScalarT M_; + }; +} + +#endif diff --git a/ComponentLib/PowerElectronicsComponents/MicrogridBusDQ/CMakeLists.txt b/ComponentLib/PowerElectronicsComponents/MicrogridBusDQ/CMakeLists.txt new file mode 100644 index 00000000..1bfb8dc3 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/MicrogridBusDQ/CMakeLists.txt @@ -0,0 +1,8 @@ + + + +gridkit_add_library(powerelec_mircobusdq + SOURCES + MicrogridBusDQ.cpp + OUTPUT_NAME + gridkit_powerelec_mircobusdq) \ No newline at end of file diff --git a/ComponentLib/PowerElectronicsComponents/MicrogridBusDQ/MicrogridBusDQ.cpp b/ComponentLib/PowerElectronicsComponents/MicrogridBusDQ/MicrogridBusDQ.cpp new file mode 100644 index 00000000..0355e44d --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/MicrogridBusDQ/MicrogridBusDQ.cpp @@ -0,0 +1,139 @@ + +#include +#include +#include +#include "MicrogridBusDQ.hpp" + +namespace ModelLib { + +/*! + * @brief Constructor for a constant MicrogridBusDQ model + * + * Calls default ModelEvaluatorImpl constructor. + * + * In DQ space + * Each microgrid line has a virtual resistance RN + * Model is from paper: " + "Modeling, Analysis and Testing of Autonomous Operation of an Inverter-Based Microgrid" Nagaraju Pogaku, Milan Prodanovic, and Timothy C. Green" + * Section E + */ + +template +MicrogridBusDQ::MicrogridBusDQ(IdxT id, ScalarT RN) + : RN_(RN) +{ + // externals [vbus_d, vbus_q] + size_ = 2; + n_intern_ = 0; + n_extern_ = 2; + extern_indices_ = {0,1}; + idc_ = id; +} + +template +MicrogridBusDQ::~MicrogridBusDQ() +{ +} + +/*! + * @brief allocate method computes sparsity pattern of the Jacobian. + */ +template +int MicrogridBusDQ::allocate() +{ + y_.resize(size_); + yp_.resize(size_); + f_.resize(size_); + + return 0; +} + +/** + * Initialization of the grid model + */ +template +int MicrogridBusDQ::initialize() +{ + return 0; +} + +/* + * \brief Identify differential variables + */ +template +int MicrogridBusDQ::tagDifferentiable() +{ + return 0; +} + +/** + * @brief Evaluate residual of microgrid line + * This model has "Virtual resistors". The voltage of the bus divided by its virtual resistance. + * The components are external to allow for outside components to add inductances to the terms. + * + * refernce to equations in class header + * + */ +template +int MicrogridBusDQ::evaluateResidual() +{ + //bus voltage + f_[0] = -y_[0] / RN_; + f_[1] = -y_[1] / RN_; + + return 0; +} + +/** + * @brief Generate Jacobian + * + * @tparam ScalarT + * @tparam IdxT + * @return int + */ +template +int MicrogridBusDQ::evaluateJacobian() +{ + J_.zeroMatrix(); + + //Create dF/dy + std::vector rtemp{0,1}; + std::vector ctemp{0,1}; + std::vector vals{-1.0 / RN_,-1.0 / RN_}; + J_.setValues(rtemp, ctemp, vals); + + return 0; +} + +template +int MicrogridBusDQ::evaluateIntegrand() +{ + return 0; +} + +template +int MicrogridBusDQ::initializeAdjoint() +{ + return 0; +} + +template +int MicrogridBusDQ::evaluateAdjointResidual() +{ + return 0; +} + +template +int MicrogridBusDQ::evaluateAdjointIntegrand() +{ + return 0; +} + + +// Available template instantiations +template class MicrogridBusDQ; +template class MicrogridBusDQ; + + +} //namespace ModelLib + diff --git a/ComponentLib/PowerElectronicsComponents/MicrogridBusDQ/MicrogridBusDQ.hpp b/ComponentLib/PowerElectronicsComponents/MicrogridBusDQ/MicrogridBusDQ.hpp new file mode 100644 index 00000000..cc5280fa --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/MicrogridBusDQ/MicrogridBusDQ.hpp @@ -0,0 +1,70 @@ + + +#ifndef _VIRBUSDQ_HPP_ +#define _VIRBUSDQ_HPP_ + +#include +#include +#include + + +namespace ModelLib +{ + template class BaseBus; +} + + +namespace ModelLib +{ + /*! + * @brief Declaration of a MicrogridBusDQ class. + * + */ + template + class MicrogridBusDQ : public CircuitComponent + { + using CircuitComponent::size_; + using CircuitComponent::nnz_; + using CircuitComponent::time_; + using CircuitComponent::alpha_; + using CircuitComponent::y_; + using CircuitComponent::yp_; + using CircuitComponent::tag_; + using CircuitComponent::f_; + using CircuitComponent::g_; + using CircuitComponent::yB_; + using CircuitComponent::ypB_; + using CircuitComponent::fB_; + using CircuitComponent::gB_; + using CircuitComponent::J_; + using CircuitComponent::param_; + using CircuitComponent::idc_; + + using CircuitComponent::extern_indices_; + using CircuitComponent::n_extern_; + using CircuitComponent::n_intern_; + + + public: + MicrogridBusDQ(IdxT id, ScalarT RN); + virtual ~MicrogridBusDQ(); + + int allocate(); + int initialize(); + int tagDifferentiable(); + int evaluateResidual(); + int evaluateJacobian(); + int evaluateIntegrand(); + + int initializeAdjoint(); + int evaluateAdjointResidual(); + //int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); + + + private: + ScalarT RN_; + }; +} + +#endif diff --git a/ComponentLib/PowerElectronicsComponents/MicrogridLine/CMakeLists.txt b/ComponentLib/PowerElectronicsComponents/MicrogridLine/CMakeLists.txt new file mode 100644 index 00000000..d59cebb9 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/MicrogridLine/CMakeLists.txt @@ -0,0 +1,8 @@ + + + +gridkit_add_library(powerelec_mircoline + SOURCES + MicrogridLine.cpp + OUTPUT_NAME + gridkit_powerelec_mircoline) \ No newline at end of file diff --git a/ComponentLib/PowerElectronicsComponents/MicrogridLine/MicrogridLine.cpp b/ComponentLib/PowerElectronicsComponents/MicrogridLine/MicrogridLine.cpp new file mode 100644 index 00000000..46c97fa2 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/MicrogridLine/MicrogridLine.cpp @@ -0,0 +1,177 @@ + +#include +#include +#include +#include "MicrogridLine.hpp" + +namespace ModelLib { + +/*! + * @brief Constructor for a constant MicrogridLine model + * + * Calls default ModelEvaluatorImpl constructor. + * + * + * Model is from paper: " + "Modeling, Analysis and Testing of Autonomous Operation of an Inverter-Based Microgrid" Nagaraju Pogaku, Milan Prodanovic, and Timothy C. Green" + * Section C + * + * @todo Consider having \omegaref as a global constant, not a node variable. + */ + +template +MicrogridLine::MicrogridLine(IdxT id, ScalarT R,ScalarT L) + : R_(R), + L_(L) +{ + // internals [id, iq] + // externals [\omegaref, vbd_in, vbq_in, vbd_out, vbq_out] + size_ = 7; + n_intern_ = 2; + n_extern_ = 5; + extern_indices_ = {0,1,2,3,4}; + idc_ = id; +} + +template +MicrogridLine::~MicrogridLine() +{ +} + +/*! + * @brief allocate method computes sparsity pattern of the Jacobian. + */ +template +int MicrogridLine::allocate() +{ + y_.resize(size_); + yp_.resize(size_); + f_.resize(size_); + + return 0; +} + +/** + * Initialization of the grid model + */ +template +int MicrogridLine::initialize() +{ + return 0; +} + +/* + * \brief Identify differential variables + */ +template +int MicrogridLine::tagDifferentiable() +{ + return 0; +} + +/** + * @brief Evaluate residual of microgrid line + * + */ +template +int MicrogridLine::evaluateResidual() +{ + //ref motor + f_[0] = 0.0; + + //input + f_[1] = -y_[5] ; + f_[2] = -y_[6] ; + + //output + f_[3] = y_[5] ; + f_[4] = y_[6] ; + + //Internal variables + f_[5] = -yp_[5] - (R_ / L_) * y_[5] + y_[0]*y_[6] + (y_[1] - y_[3])/L_; + f_[6] = -yp_[6] - (R_ / L_) * y_[6] - y_[0]*y_[5] + (y_[2] - y_[4])/L_; + + + return 0; +} + +/** + * @brief Generate Jacobian for Microgrid Line + * + * @tparam ScalarT + * @tparam IdxT + * @return int + */ +template +int MicrogridLine::evaluateJacobian() +{ + J_.zeroMatrix(); + + //Create dF/dy + std::vector rtemp{1,2,3,4}; + std::vector ctemp{5,6,5,6}; + std::vector vals{-1.0,-1.0,1.0,1.0}; + J_.setValues(rtemp, ctemp, vals); + + + std::vector ccord{0, 1, 3, 5, 6}; + + std::vector rcord(ccord.size(),5); + vals = {y_[6], (1.0 / L_) , -(1.0 / L_), - (R_ / L_) , y_[0]}; + J_.setValues(rcord, ccord, vals); + + + std::vector ccor2{0, 2, 4, 5, 6}; + std::fill(rcord.begin(), rcord.end(), 6); + vals = {-y_[5], (1.0 / L_) , -(1.0 / L_), -y_[0], - (R_ / L_)}; + J_.setValues(rcord, ccor2, vals); + + + //Create -dF/dy' + std::vector rcordder{5,6}; + std::vector ccordder{5,6}; + std::vector valsder{-1.0, -1.0}; + COO_Matrix Jacder = COO_Matrix(rcordder, ccordder, valsder,7,7); + + //Perform dF/dy + \alpha dF/dy' + J_.axpy(alpha_, Jacder); + + + return 0; +} + +template +int MicrogridLine::evaluateIntegrand() +{ + return 0; +} + +template +int MicrogridLine::initializeAdjoint() +{ + return 0; +} + +template +int MicrogridLine::evaluateAdjointResidual() +{ + return 0; +} + +template +int MicrogridLine::evaluateAdjointIntegrand() +{ + return 0; +} + + + + + +// Available template instantiations +template class MicrogridLine; +template class MicrogridLine; + + +} //namespace ModelLib + diff --git a/ComponentLib/PowerElectronicsComponents/MicrogridLine/MicrogridLine.hpp b/ComponentLib/PowerElectronicsComponents/MicrogridLine/MicrogridLine.hpp new file mode 100644 index 00000000..362cb4b6 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/MicrogridLine/MicrogridLine.hpp @@ -0,0 +1,68 @@ + + +#ifndef _TRANLINE_HPP_ +#define _TRANLINE_HPP_ + +#include +#include +#include + +namespace ModelLib +{ + template + class BaseBus; +} + +namespace ModelLib +{ + /*! + * @brief Declaration of a MicrogridLine class. + * + */ + template + class MicrogridLine : public CircuitComponent + { + using CircuitComponent::size_; + using CircuitComponent::nnz_; + using CircuitComponent::time_; + using CircuitComponent::alpha_; + using CircuitComponent::y_; + using CircuitComponent::yp_; + using CircuitComponent::tag_; + using CircuitComponent::f_; + using CircuitComponent::g_; + using CircuitComponent::yB_; + using CircuitComponent::ypB_; + using CircuitComponent::fB_; + using CircuitComponent::gB_; + using CircuitComponent::J_; + using CircuitComponent::param_; + using CircuitComponent::idc_; + + using CircuitComponent::extern_indices_; + using CircuitComponent::n_extern_; + using CircuitComponent::n_intern_; + + public: + MicrogridLine(IdxT id, ScalarT R, ScalarT L); + virtual ~MicrogridLine(); + + int allocate(); + int initialize(); + int tagDifferentiable(); + int evaluateResidual(); + int evaluateJacobian(); + int evaluateIntegrand(); + + int initializeAdjoint(); + int evaluateAdjointResidual(); + // int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); + + private: + ScalarT R_; + ScalarT L_; + }; +} + +#endif diff --git a/ComponentLib/PowerElectronicsComponents/MicrogridLoad/CMakeLists.txt b/ComponentLib/PowerElectronicsComponents/MicrogridLoad/CMakeLists.txt new file mode 100644 index 00000000..313d4768 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/MicrogridLoad/CMakeLists.txt @@ -0,0 +1,8 @@ + + + +gridkit_add_library(powerelec_microload + SOURCES + MicrogridLoad.cpp + OUTPUT_NAME + gridkit_powerelec_microload) \ No newline at end of file diff --git a/ComponentLib/PowerElectronicsComponents/MicrogridLoad/MicrogridLoad.cpp b/ComponentLib/PowerElectronicsComponents/MicrogridLoad/MicrogridLoad.cpp new file mode 100644 index 00000000..ddfa9ab9 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/MicrogridLoad/MicrogridLoad.cpp @@ -0,0 +1,172 @@ + +#include +#include +#include +#include "MicrogridLoad.hpp" + +namespace ModelLib { + +/*! + * @brief Constructor for a constant MicrogridLoad model + * + * Calls default ModelEvaluatorImpl constructor. + * + * + * Model is from paper: " + "Modeling, Analysis and Testing of Autonomous Operation of an Inverter-Based Microgrid" Nagaraju Pogaku, Milan Prodanovic, and Timothy C. Green" + * Section D + */ + +template +MicrogridLoad::MicrogridLoad(IdxT id, ScalarT R,ScalarT L) + : R_(R), + L_(L) +{ + // internals [id, iq] + // externals [\omegaref, vbd_out, vbq_out] + size_ = 5; + n_intern_ = 2; + n_extern_ = 3; + extern_indices_ = {0,1,2}; + idc_ = id; + +} + +template +MicrogridLoad::~MicrogridLoad() +{ +} + +/*! + * @brief allocate method computes sparsity pattern of the Jacobian. + */ +template +int MicrogridLoad::allocate() +{ + y_.resize(size_); + yp_.resize(size_); + f_.resize(size_); + + return 0; +} + +/** + * Initialization of the grid model + */ +template +int MicrogridLoad::initialize() +{ + return 0; +} + +/* + * \brief Identify differential variables + */ +template +int MicrogridLoad::tagDifferentiable() +{ + return 0; +} + +/** + * @brief Eval Micro Load + */ +template +int MicrogridLoad::evaluateResidual() +{ + //ref motor + f_[0] = 0.0; + + //only input for loads + + //input + f_[1] = -y_[3] ; + f_[2] = -y_[4] ; + + //Internal variables + f_[3] = -yp_[3] - (R_ / L_) * y_[3] + y_[0]*y_[4] + y_[1] / L_; + f_[4] = -yp_[4] - (R_ / L_) * y_[4] - y_[0]*y_[3] + y_[2] / L_; + + + return 0; +} + +/** + * @brief Generate Jacobian for Micro Load + * + * @tparam ScalarT + * @tparam IdxT + * @return int + */ +template +int MicrogridLoad::evaluateJacobian() +{ + J_.zeroMatrix(); + + //Create dF/dy + std::vector rtemp{1,2}; + std::vector ctemp{3,4}; + std::vector vals{-1.0,-1.0}; + J_.setValues(rtemp, ctemp, vals); + + + std::vector ccord{0, 1, 3, 4}; + + std::vector rcord(ccord.size(),3); + vals = {y_[4], (1.0 / L_) , - (R_ / L_) , y_[0]}; + J_.setValues(rcord, ccord, vals); + + + std::vector ccor2{0, 2, 3, 4}; + std::fill(rcord.begin(), rcord.end(), 4); + vals = {-y_[3], (1.0 / L_) , -y_[0], - (R_ / L_)}; + J_.setValues(rcord, ccor2, vals); + + + //Create -dF/dy' + std::vector rcordder{3,4}; + std::vector ccordder{3,4}; + std::vector valsder{-1.0, -1.0}; + COO_Matrix Jacder = COO_Matrix(rcordder, ccordder, valsder,5,5); + + //Perform dF/dy + \alpha dF/dy' + J_.axpy(alpha_, Jacder); + + return 0; +} + +template +int MicrogridLoad::evaluateIntegrand() +{ + return 0; +} + +template +int MicrogridLoad::initializeAdjoint() +{ + return 0; +} + +template +int MicrogridLoad::evaluateAdjointResidual() +{ + return 0; +} + +template +int MicrogridLoad::evaluateAdjointIntegrand() +{ + return 0; +} + + + + + +// Available template instantiations +template class MicrogridLoad; +template class MicrogridLoad; + + +} //namespace ModelLib + diff --git a/ComponentLib/PowerElectronicsComponents/MicrogridLoad/MicrogridLoad.hpp b/ComponentLib/PowerElectronicsComponents/MicrogridLoad/MicrogridLoad.hpp new file mode 100644 index 00000000..69ad822a --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/MicrogridLoad/MicrogridLoad.hpp @@ -0,0 +1,71 @@ + + +#ifndef _TRANLOAD_HPP_ +#define _TRANLOAD_HPP_ + +#include +#include +#include + + +namespace ModelLib +{ + template class BaseBus; +} + + +namespace ModelLib +{ + /*! + * @brief Declaration of a passive MicrogridLoad class. + * + */ + template + class MicrogridLoad : public CircuitComponent + { + using CircuitComponent::size_; + using CircuitComponent::nnz_; + using CircuitComponent::time_; + using CircuitComponent::alpha_; + using CircuitComponent::y_; + using CircuitComponent::yp_; + using CircuitComponent::tag_; + using CircuitComponent::f_; + using CircuitComponent::g_; + using CircuitComponent::yB_; + using CircuitComponent::ypB_; + using CircuitComponent::fB_; + using CircuitComponent::gB_; + using CircuitComponent::J_; + using CircuitComponent::param_; + using CircuitComponent::idc_; + + using CircuitComponent::extern_indices_; + using CircuitComponent::n_extern_; + using CircuitComponent::n_intern_; + + + public: + MicrogridLoad(IdxT id, ScalarT R, ScalarT L); + virtual ~MicrogridLoad(); + + int allocate(); + int initialize(); + int tagDifferentiable(); + int evaluateResidual(); + int evaluateJacobian(); + int evaluateIntegrand(); + + int initializeAdjoint(); + int evaluateAdjointResidual(); + //int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); + + + private: + ScalarT R_; + ScalarT L_; + }; +} + +#endif diff --git a/ComponentLib/PowerElectronicsComponents/Resistor/CMakeLists.txt b/ComponentLib/PowerElectronicsComponents/Resistor/CMakeLists.txt new file mode 100644 index 00000000..9386bda8 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/Resistor/CMakeLists.txt @@ -0,0 +1,8 @@ + + + +gridkit_add_library(powerelec_resistor + SOURCES + Resistor.cpp + OUTPUT_NAME + gridkit_powerelec_resistor) \ No newline at end of file diff --git a/ComponentLib/PowerElectronicsComponents/Resistor/Resistor.cpp b/ComponentLib/PowerElectronicsComponents/Resistor/Resistor.cpp new file mode 100644 index 00000000..caf0335f --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/Resistor/Resistor.cpp @@ -0,0 +1,126 @@ + + + +#include +#include +#include +#include "Resistor.hpp" + +namespace ModelLib { + +/*! + * @brief Constructor for a resistor model + * + * Calls default ModelEvaluatorImpl constructor. + */ + +template +Resistor::Resistor(IdxT id, ScalarT R) + : R_(R) +{ + size_ = 2; + n_intern_ = 0; + n_extern_ = 2; + extern_indices_ = {0,1}; + idc_ = id; +} + +template +Resistor::~Resistor() +{ +} + +/*! + * @brief allocate method computes sparsity pattern of the Jacobian. + */ +template +int Resistor::allocate() +{ + y_.resize(size_); + yp_.resize(size_); + f_.resize(size_); + + return 0; +} + +/** + * Initialization of the grid model + */ +template +int Resistor::initialize() +{ + return 0; +} + +/* + * \brief Identify differential variables + */ +template +int Resistor::tagDifferentiable() +{ + return 0; +} + +/** + * @brief Computes the resistors resisdual + * + */ +template +int Resistor::evaluateResidual() +{ + //input + f_[0] = (y_[0] - y_[1])/R_ ; + //ouput + f_[1] = (y_[1] - y_[0])/R_ ; + return 0; +} + +template +int Resistor::evaluateJacobian() +{ + + //Create dF/dy + //does compiler make constant??? + std::vector rcord{0,0,1,1}; + std::vector ccord{0,1,0,1}; + std::vector vals{1.0 / R_, -1.0 / R_, -1.0 / R_, 1.0 / R_}; + J_.setValues(rcord, ccord, vals); + + return 0; +} + +template +int Resistor::evaluateIntegrand() +{ + return 0; +} + +template +int Resistor::initializeAdjoint() +{ + return 0; +} + +template +int Resistor::evaluateAdjointResidual() +{ + return 0; +} + +template +int Resistor::evaluateAdjointIntegrand() +{ + return 0; +} + + + + + +// Available template instantiations +template class Resistor; +template class Resistor; + + +} //namespace ModelLib + diff --git a/ComponentLib/PowerElectronicsComponents/Resistor/Resistor.hpp b/ComponentLib/PowerElectronicsComponents/Resistor/Resistor.hpp new file mode 100644 index 00000000..5ec5587b --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/Resistor/Resistor.hpp @@ -0,0 +1,71 @@ + + +#ifndef _RES_HPP_ +#define _RES_HPP_ + +#include +#include +#include + + +namespace ModelLib +{ + template class BaseBus; +} + + +namespace ModelLib +{ + /*! + * @brief Declaration of a Resistor class. + * + */ + template + class Resistor : public CircuitComponent + { + using CircuitComponent::size_; + using CircuitComponent::nnz_; + using CircuitComponent::time_; + using CircuitComponent::alpha_; + using CircuitComponent::y_; + using CircuitComponent::yp_; + using CircuitComponent::tag_; + using CircuitComponent::f_; + using CircuitComponent::g_; + using CircuitComponent::yB_; + using CircuitComponent::ypB_; + using CircuitComponent::fB_; + using CircuitComponent::gB_; + using CircuitComponent::J_; + using CircuitComponent::param_; + using CircuitComponent::idc_; + + + using CircuitComponent::extern_indices_; + using CircuitComponent::n_extern_; + using CircuitComponent::n_intern_; + + + public: + Resistor(IdxT id, ScalarT R); + virtual ~Resistor(); + + int allocate(); + int initialize(); + int tagDifferentiable(); + int evaluateResidual(); + int evaluateJacobian(); + int evaluateIntegrand(); + + int initializeAdjoint(); + int evaluateAdjointResidual(); + //int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); + + + private: + ScalarT R_; + }; +} + +#endif diff --git a/ComponentLib/PowerElectronicsComponents/SynchronousMachine/CMakeLists.txt b/ComponentLib/PowerElectronicsComponents/SynchronousMachine/CMakeLists.txt new file mode 100644 index 00000000..bfaeace4 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/SynchronousMachine/CMakeLists.txt @@ -0,0 +1,8 @@ + + + +gridkit_add_library(powerelec_synmachine + SOURCES + SynchronousMachine.cpp + OUTPUT_NAME + gridkit_powerelec_synmachine) \ No newline at end of file diff --git a/ComponentLib/PowerElectronicsComponents/SynchronousMachine/SynchronousMachine.cpp b/ComponentLib/PowerElectronicsComponents/SynchronousMachine/SynchronousMachine.cpp new file mode 100644 index 00000000..e3f65bed --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/SynchronousMachine/SynchronousMachine.cpp @@ -0,0 +1,153 @@ + + + +#include +#include +#include +#include "SynchronousMachine.hpp" + + +namespace ModelLib { + + + +/*! + * @brief Constructor for a constant SynchronousMachine model + * + * Calls default ModelEvaluatorImpl constructor. + * @todo This models equations are not finish + * @todo needs to be tested for correctness + */ + +template +SynchronousMachine::SynchronousMachine(IdxT id, ScalarT Lls, std::tuple Llkq, ScalarT Llfd, ScalarT Llkd, ScalarT Lmq, ScalarT Lmd, ScalarT Rs, std::tuple Rkq, ScalarT Rfd, ScalarT Rkd, ScalarT RJ, ScalarT P, ScalarT mub) + : Lls_(Lls), + Llkq_(Llkq), + Llfd_(Llfd), + Llkd_(Llkd), + Lmq_(Lmq), + Lmd_(Lmd), + Rs_(Rs), + Rkq_(Rkq), + Rfd_(Rfd), + Rkd_(Rkd), + RJ_(RJ), + P_(P), + mub_(mub) +{ + size_ = 13; + n_intern_ = 6; + n_extern_ = 7; + extern_indices_ = {0,1,2,3,4}; + idc_ = id; +} + +template +SynchronousMachine::~SynchronousMachine() +{ +} + +/*! + * @brief allocate method computes sparsity pattern of the Jacobian. + */ +template +int SynchronousMachine::allocate() +{ + y_.resize(size_); + yp_.resize(size_); + f_.resize(size_); + return 0; +} + +/** + * Initialization of the grid model + */ +template +int SynchronousMachine::initialize() +{ + return 0; +} + +/* + * \brief Identify differential variables + */ +template +int SynchronousMachine::tagDifferentiable() +{ + return 0; +} + +/** + * @brief Compute the resisdual of the component. + * + * @todo not finished + */ +template +int SynchronousMachine::evaluateResidual() +{ + ScalarT rkq1 = std::get<0>(Rkq_); + ScalarT rkq2 = std::get<1>(Rkq_); + ScalarT llkq1 = std::get<0>(Llkq_); + ScalarT llkq2 = std::get<1>(Llkq_); + + ScalarT cos1 = cos((P_/2.0)*y_[5]); + ScalarT sin1 = sin((P_/2.0)*y_[5]); + ScalarT cos23m = cos((P_/2.0)*y_[5] - (2.0/3.0)*M_PI); + ScalarT sin23m = sin((P_/2.0)*y_[5] - (2.0/3.0)*M_PI); + ScalarT cos23p = cos((P_/2.0)*y_[5] + (2.0/3.0)*M_PI); + ScalarT sin23p = sin((P_/2.0)*y_[5] + (2.0/3.0)*M_PI); + + f_[0] = y_[6]*cos1 + y_[7]*sin1 + y_[8]; + f_[1] = y_[6]*cos23m + y_[7]*sin23m + y_[8]; + f_[2] = y_[6]*cos23p + y_[7]*sin23p + y_[8]; + f_[3] = RJ_ * yp_[4] - (3.0/4.0)*P_*(Lmd_ *y_[6]* (y_[7] + y_[11] + y_[12]) - Lmq_*y_[7]*(y_[6] + y_[9] + y_[0])); + f_[4] = yp_[5] - y_[4]; + f_[5] = (-2.0/3.0)*(y_[0]*cos1 +y_[1]*cos23m + y_[2]*cos23p) + Rs_*y_[6] + (Lls_ + Lmq_)*yp_[6] + Lmq_*yp_[9] + Lmq_*yp_[10] + y_[4]*(P_/2.0)*((Lls_ + Lmd_)*y_[7] + Lmd_*y_[11] + Lmd_*y_[12]); + f_[6] = (-2.0/3.0)*(y_[0]*sin1 -y_[1]*sin23m - y_[2]*sin23p) + Rs_*y_[7] + (Lls_ + Lmd_)*yp_[7] + Lmd_*yp_[11] + Lmd_*yp_[12] - y_[4]*(P_/2.0)*((Lls_ + Lmq_)*y_[6] + Lmq_*y_[9] + Lmq_*y_[10]); + f_[7] = (-1.0/3.0)*(y_[0] + y_[1] + y_[2]) + Rs_*y_[8] + Lls_*yp_[8]; + f_[8] = rkq1*y_[9] + (llkq1 + Lmq_)*yp_[9] + Lmq_*yp_[6] + Lmq_*yp_[10]; + f_[9] = rkq1*y_[9] + (llkq1 + Lmq_)*yp_[9] + Lmq_*yp_[6] + Lmq_*yp_[10]; + return 0; +} + +template +int SynchronousMachine::evaluateJacobian() +{ + return 0; +} + +template +int SynchronousMachine::evaluateIntegrand() +{ + return 0; +} + +template +int SynchronousMachine::initializeAdjoint() +{ + return 0; +} + +template +int SynchronousMachine::evaluateAdjointResidual() +{ + return 0; +} + +template +int SynchronousMachine::evaluateAdjointIntegrand() +{ + return 0; +} + + + + + +// Available template instantiations +template class SynchronousMachine; +template class SynchronousMachine; + + +} //namespace ModelLib + diff --git a/ComponentLib/PowerElectronicsComponents/SynchronousMachine/SynchronousMachine.hpp b/ComponentLib/PowerElectronicsComponents/SynchronousMachine/SynchronousMachine.hpp new file mode 100644 index 00000000..4681cb6f --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/SynchronousMachine/SynchronousMachine.hpp @@ -0,0 +1,81 @@ + + +#ifndef _SYNMACH_HPP_ +#define _SYNMACH_HPP_ + +#include +#include +#include +#include + + +namespace ModelLib +{ + template class BaseBus; +} + + +namespace ModelLib +{ + /*! + * @brief Declaration of a SynchronousMachine class. + * + */ + template + class SynchronousMachine : public CircuitComponent + { + using CircuitComponent::size_; + using CircuitComponent::nnz_; + using CircuitComponent::time_; + using CircuitComponent::alpha_; + using CircuitComponent::y_; + using CircuitComponent::yp_; + using CircuitComponent::tag_; + using CircuitComponent::f_; + using CircuitComponent::g_; + using CircuitComponent::yB_; + using CircuitComponent::ypB_; + using CircuitComponent::fB_; + using CircuitComponent::gB_; + using CircuitComponent::J_; + using CircuitComponent::param_; + using CircuitComponent::idc_; + + using CircuitComponent::extern_indices_; + using CircuitComponent::n_extern_; + using CircuitComponent::n_intern_; + + public: + SynchronousMachine(IdxT id, ScalarT Lls, std::tuple Llkq, ScalarT Llfd, ScalarT Llkd, ScalarT Lmq, ScalarT Lmd, ScalarT Rs, std::tuple Rkq, ScalarT Rfd, ScalarT Rkd, ScalarT RJ, ScalarT P, ScalarT mub); + virtual ~SynchronousMachine(); + + int allocate(); + int initialize(); + int tagDifferentiable(); + int evaluateResidual(); + int evaluateJacobian(); + int evaluateIntegrand(); + + int initializeAdjoint(); + int evaluateAdjointResidual(); + //int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); + + private: + ScalarT Lls_; + std::tuple Llkq_; + ScalarT Llfd_; + ScalarT Llkd_; + ScalarT Lmq_; + ScalarT Lmd_; + ScalarT Rs_; + std::tuple Rkq_; + ScalarT Rfd_; + ScalarT Rkd_; + ScalarT RJ_; + ScalarT P_; + ScalarT mub_; + }; +} + +#endif diff --git a/ComponentLib/PowerElectronicsComponents/TransmissionLine/CMakeLists.txt b/ComponentLib/PowerElectronicsComponents/TransmissionLine/CMakeLists.txt new file mode 100644 index 00000000..6e2cea4f --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/TransmissionLine/CMakeLists.txt @@ -0,0 +1,8 @@ + + + +gridkit_add_library(powerelec_tranline + SOURCES + TransmissionLine.cpp + OUTPUT_NAME + gridkit_powerelec_tranline) \ No newline at end of file diff --git a/ComponentLib/PowerElectronicsComponents/TransmissionLine/TransmissionLine.cpp b/ComponentLib/PowerElectronicsComponents/TransmissionLine/TransmissionLine.cpp new file mode 100644 index 00000000..010f38e1 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/TransmissionLine/TransmissionLine.cpp @@ -0,0 +1,205 @@ + +#include +#include +#include +#include "TransmissionLine.hpp" + +namespace ModelLib { + +/*! + * @brief Constructor for a TransmissionLine model + * + * Calls default ModelEvaluatorImpl constructor. + * + * This is the Medium distance form with the use of the admittance matrix. + * Since the line is of medium length then there is no real part for shunt admittance + * @todo needs to used in a model + * @todo test for correctness + */ + +template +TransmissionLine::TransmissionLine(IdxT id, ScalarT R,ScalarT X, ScalarT B) + : R_(R), + X_(X), + B_(B) +{ + // internals [Iret1, Iimt1, Iret2, Iimt2] + // externals [Vre11, Vim11, Vre12, Vim12, Vre21, Vim21, Vre22, Vim22] + size_ = 12; + n_intern_ = 4; + n_extern_ = 8; + extern_indices_ = {0,1,2,3,4,5,6,7}; + idc_ = id; + + ScalarT magImpendence = 1 / (R_*R_ + X_*X_); + YReMat_ = magImpendence * R_; + YImMatOff_ = magImpendence * X_; + YImMatDi_ = B_ / (2.0) - YImMatOff_; +} + +template +TransmissionLine::~TransmissionLine() +{ +} + +/*! + * @brief allocate method computes sparsity pattern of the Jacobian. + */ +template +int TransmissionLine::allocate() +{ + y_.resize(size_); + yp_.resize(size_); + f_.resize(size_); + + return 0; +} + +/** + * Initialization of the grid model + */ +template +int TransmissionLine::initialize() +{ + return 0; +} + +/* + * \brief Identify differential variables + */ +template +int TransmissionLine::tagDifferentiable() +{ + return 0; +} + +/** + * @brief Evaluate residual of transmission line + * + * The complex admittance matrix is: + * [[ Y/2 + 1/Z, -1/Z]; + * [ -1/Z, Y/2 + 1/Z]] = + * [[R/|Z|, -R/|Z|]; + * [-R/|Z|, R/|Z|]] + + * i [[B/2 - X/|Z|, X/|Z|]; + * [X/|Z|, B/2 - X/|Z|]] + * = Dre + i Dim + * + * Then + * Ire = Dre Vre - Dim Vim + * Iim = Dre Vim + Dim Vre + * + * To express this for Modified Nodal Analysis the Voltages of the admittance matrix are put into voltage drops + */ +template +int TransmissionLine::evaluateResidual() +{ + //input + f_[0] = y_[8] ; + f_[1] = y_[9] ; + + f_[2] = y_[10] ; + f_[3] = y_[11] ; + //ouput + f_[4] = -y_[8] ; + f_[5] = -y_[9] ; + + f_[6] = -y_[10] ; + f_[7] = -y_[11] ; + + //Voltage drop accross terminals + ScalarT V1re = y_[0] - y_[4]; + ScalarT V1im = y_[1] - y_[5]; + ScalarT V2re = y_[2] - y_[6]; + ScalarT V2im = y_[3] - y_[7]; + + //Internal variables + //row 1 + f_[8] = YReMat_ * (V1re - V2re) - (YImMatDi_ * V1im + YImMatOff_ * V2im) - y_[8] ; + f_[9] = YReMat_ * (V1im - V2im) + (YImMatDi_ * V1re + YImMatOff_ * V2re) - y_[9] ; + + //row2 + f_[10] = YReMat_ * (V2re - V1re) - (YImMatOff_ * V1im + YImMatDi_ * V2im) - y_[10]; + f_[11] = YReMat_ * (V2im - V1im) + (YImMatOff_ * V1re + YImMatDi_ * V2re) - y_[11]; + + + return 0; +} + +/** + * @brief Generate Jacobian for Transmission Line + * + * @tparam ScalarT + * @tparam IdxT + * @return int + */ +template +int TransmissionLine::evaluateJacobian() +{ + + //Create dF/dy + std::vector rtemp{0,1,2,3,4,5,6,7,8,9,10,11}; + std::vector ctemp{8,9,10,11,8,9,10,11,8,9,10,11}; + std::vector vals{1.0,1.0,1.0,1.0,-1.0,-1.0,-1.0,-1.0,-1.0,-1.0,-1.0,-1.0}; + J_.setValues(rtemp, ctemp, vals); + + + std::vector ccord{0,1,2,3,4,5,6,7}; + + std::vector rcord(ccord.size(),8); + vals = {YReMat_, -YImMatDi_ ,-YReMat_, -YImMatOff_,-YReMat_, YImMatDi_ ,YReMat_, YImMatOff_}; + J_.setValues(rtemp, ctemp, vals); + + + std::fill(rcord.begin(), rcord.end(), 9); + vals = {YImMatDi_ ,YReMat_, YImMatOff_, -YReMat_,-YImMatDi_ ,-YReMat_, -YImMatOff_, YReMat_}; + J_.setValues(rtemp, ctemp, vals); + + + std::fill(rcord.begin(), rcord.end(), 10); + vals = {-YReMat_, -YImMatDi_ ,YReMat_, -YImMatOff_,YReMat_, YImMatDi_ ,-YReMat_, YImMatOff_}; + J_.setValues(rtemp, ctemp, vals); + + + std::fill(rcord.begin(), rcord.end(), 11); + vals = {YImMatDi_ ,-YReMat_, YImMatOff_, YReMat_,-YImMatDi_ ,YReMat_, -YImMatOff_, -YReMat_}; + J_.setValues(rtemp, ctemp, vals); + + return 0; +} + +template +int TransmissionLine::evaluateIntegrand() +{ + return 0; +} + +template +int TransmissionLine::initializeAdjoint() +{ + return 0; +} + +template +int TransmissionLine::evaluateAdjointResidual() +{ + return 0; +} + +template +int TransmissionLine::evaluateAdjointIntegrand() +{ + return 0; +} + + + + + +// Available template instantiations +template class TransmissionLine; +template class TransmissionLine; + + +} //namespace ModelLib + diff --git a/ComponentLib/PowerElectronicsComponents/TransmissionLine/TransmissionLine.hpp b/ComponentLib/PowerElectronicsComponents/TransmissionLine/TransmissionLine.hpp new file mode 100644 index 00000000..63cd1e88 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/TransmissionLine/TransmissionLine.hpp @@ -0,0 +1,78 @@ + + +#ifndef _TRANLINE_HPP_ +#define _TRANLINE_HPP_ + +#include +#include +#include + + +namespace ModelLib +{ + template class BaseBus; +} + + +namespace ModelLib +{ + /*! + * @brief Declaration of a TransmissionLine class. + * + * Model from Adam Birchfield paper (medium distances < 2km). + * See also textbooks "Power System Analysis" by Grainger and "Power System Dynamics and Stability" by Sauer & Pai + * + * @note Not used in the Microgrid model. + */ + template + class TransmissionLine : public CircuitComponent + { + using CircuitComponent::size_; + using CircuitComponent::nnz_; + using CircuitComponent::time_; + using CircuitComponent::alpha_; + using CircuitComponent::y_; + using CircuitComponent::yp_; + using CircuitComponent::tag_; + using CircuitComponent::f_; + using CircuitComponent::g_; + using CircuitComponent::yB_; + using CircuitComponent::ypB_; + using CircuitComponent::fB_; + using CircuitComponent::gB_; + using CircuitComponent::J_; + using CircuitComponent::param_; + using CircuitComponent::idc_; + + using CircuitComponent::extern_indices_; + using CircuitComponent::n_extern_; + using CircuitComponent::n_intern_; + + public: + TransmissionLine(IdxT id, ScalarT R, ScalarT X, ScalarT B); + virtual ~TransmissionLine(); + + int allocate(); + int initialize(); + int tagDifferentiable(); + int evaluateResidual(); + int evaluateJacobian(); + int evaluateIntegrand(); + + int initializeAdjoint(); + int evaluateAdjointResidual(); + //int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); + + + private: + ScalarT R_; + ScalarT X_; + ScalarT B_; + ScalarT YReMat_; + ScalarT YImMatDi_; + ScalarT YImMatOff_; + }; +} + +#endif diff --git a/ComponentLib/PowerElectronicsComponents/VoltageSource/CMakeLists.txt b/ComponentLib/PowerElectronicsComponents/VoltageSource/CMakeLists.txt new file mode 100644 index 00000000..7196f4d4 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/VoltageSource/CMakeLists.txt @@ -0,0 +1,8 @@ + + + +gridkit_add_library(powerelec_voltagesource + SOURCES + VoltageSource.cpp + OUTPUT_NAME + gridkit_powerelec_voltagesource) \ No newline at end of file diff --git a/ComponentLib/PowerElectronicsComponents/VoltageSource/VoltageSource.cpp b/ComponentLib/PowerElectronicsComponents/VoltageSource/VoltageSource.cpp new file mode 100644 index 00000000..2e404151 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/VoltageSource/VoltageSource.cpp @@ -0,0 +1,122 @@ + + + +#include +#include +#include +#include "VoltageSource.hpp" + +namespace ModelLib { + +/*! + * @brief Constructor for a constant VoltageSource model + * + * Calls default ModelEvaluatorImpl constructor. + */ + +template +VoltageSource::VoltageSource(IdxT id, ScalarT V) + : V_(V) +{ + size_ = 3; + n_intern_ = 1; + n_extern_ = 2; + extern_indices_ = {0,1}; + idc_ = id; +} + +template +VoltageSource::~VoltageSource() +{ +} + +/*! + * @brief allocate method computes sparsity pattern of the Jacobian. + */ +template +int VoltageSource::allocate() +{ + y_.resize(size_); + yp_.resize(size_); + f_.resize(size_); + + return 0; +} + +/** + * Initialization of the grid model + */ +template +int VoltageSource::initialize() +{ + return 0; +} + +/* + * \brief Identify differential variables + */ +template +int VoltageSource::tagDifferentiable() +{ + return 0; +} + +/** + * @brief Evaluate resisdual of component + */ +template +int VoltageSource::evaluateResidual() +{ + //input + f_[0] = -y_[2]; + //ouput + f_[1] = y_[2]; + //internal + f_[2] = y_[1] - y_[0] - V_; + return 0; +} + +template +int VoltageSource::evaluateJacobian() +{ + //Create dF/dy + std::vector rcord{0,1,2,2}; + std::vector ccord{2,2,0,1}; + std::vector vals{-1.0, 1.0, -1.0, 1.0}; + J_.setValues(rcord, ccord, vals); + + return 0; +} + +template +int VoltageSource::evaluateIntegrand() +{ + return 0; +} + +template +int VoltageSource::initializeAdjoint() +{ + return 0; +} + +template +int VoltageSource::evaluateAdjointResidual() +{ + return 0; +} + +template +int VoltageSource::evaluateAdjointIntegrand() +{ + return 0; +} + + +// Available template instantiations +template class VoltageSource; +template class VoltageSource; + + +} //namespace ModelLib + diff --git a/ComponentLib/PowerElectronicsComponents/VoltageSource/VoltageSource.hpp b/ComponentLib/PowerElectronicsComponents/VoltageSource/VoltageSource.hpp new file mode 100644 index 00000000..89c13d3f --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/VoltageSource/VoltageSource.hpp @@ -0,0 +1,68 @@ + + +#ifndef _VOSO_HPP_ +#define _VOSO_HPP_ + +#include +#include +#include + + +namespace ModelLib +{ + template class BaseBus; +} + + +namespace ModelLib +{ + /*! + * @brief Declaration of a VoltageSource class. + * + */ + template + class VoltageSource : public CircuitComponent + { + using CircuitComponent::size_; + using CircuitComponent::nnz_; + using CircuitComponent::time_; + using CircuitComponent::alpha_; + using CircuitComponent::y_; + using CircuitComponent::yp_; + using CircuitComponent::tag_; + using CircuitComponent::f_; + using CircuitComponent::g_; + using CircuitComponent::yB_; + using CircuitComponent::ypB_; + using CircuitComponent::fB_; + using CircuitComponent::gB_; + using CircuitComponent::J_; + using CircuitComponent::param_; + using CircuitComponent::idc_; + + using CircuitComponent::extern_indices_; + using CircuitComponent::n_extern_; + using CircuitComponent::n_intern_; + + public: + VoltageSource(IdxT id, ScalarT V); + virtual ~VoltageSource(); + + int allocate(); + int initialize(); + int tagDifferentiable(); + int evaluateResidual(); + int evaluateJacobian(); + int evaluateIntegrand(); + + int initializeAdjoint(); + int evaluateAdjointResidual(); + //int evaluateAdjointJacobian(); + int evaluateAdjointIntegrand(); + + private: + ScalarT V_; + }; +} + +#endif diff --git a/Examples/CMakeLists.txt b/Examples/CMakeLists.txt index feec6623..0965f568 100644 --- a/Examples/CMakeLists.txt +++ b/Examples/CMakeLists.txt @@ -56,6 +56,10 @@ # add_subdirectory(MatPowerTesting) +add_subdirectory(RLCircuit) +add_subdirectory(Microgrid) +add_subdirectory(SparseTest) +add_subdirectory(DistributedGeneratorTest) if(TARGET SUNDIALS::kinsol) add_subdirectory(Grid3Bus) diff --git a/Examples/DistributedGeneratorTest/CMakeLists.txt b/Examples/DistributedGeneratorTest/CMakeLists.txt new file mode 100644 index 00000000..41b124cb --- /dev/null +++ b/Examples/DistributedGeneratorTest/CMakeLists.txt @@ -0,0 +1,12 @@ + + + + +add_executable(dgtest DGTest.cpp) +target_link_libraries(dgtest GRIDKIT::powerelec_disgen + GRIDKIT::powerelec_mircoline + GRIDKIT::powerelec_microload + GRIDKIT::solvers_dyn) + +add_test(NAME DistributedGeneratorTest COMMAND $) +install(TARGETS dgtest RUNTIME DESTINATION bin) diff --git a/Examples/DistributedGeneratorTest/DGTest.cpp b/Examples/DistributedGeneratorTest/DGTest.cpp new file mode 100644 index 00000000..d1f58c5c --- /dev/null +++ b/Examples/DistributedGeneratorTest/DGTest.cpp @@ -0,0 +1,85 @@ + + +#include +#include +#include +#include +#include +#include + +#include + + +/** + * @brief Testing for the Distributed Generators outputs + * + * @param argc + * @param argv + * @return int + */ +int main(int argc, char const *argv[]) +{ + + ModelLib::DistributedGeneratorParameters parms; + //Parameters from MATLAB Microgrid code for first DG + parms.wb_ = 2.0*M_PI*50.0; + parms.wc_ = 31.41; + parms.mp_ = 9.4e-5; + parms.Vn_ = 380; + parms.nq_ = 1.3e-3; + parms.F_ = 0.75; + parms.Kiv_ = 420.0; + parms.Kpv_ = 0.1; + parms.Kic_ = 20.0 * 1.0e3; + parms.Kpc_ = 15.0; + parms.Cf_ = 50.0e-6; + parms.rLf_ = 0.1; + parms.Lf_ = 1.35e-3; + parms.rLc_ = 0.03; + parms.Lc_ = 0.35e-3; + + ModelLib::DistributedGenerator *dg = new ModelLib::DistributedGenerator(0, parms, true); + + std::vector t1(16,0.0); + std::vector t2{0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5}; + + dg->allocate(); + + dg->y() = t2; + dg->yp() = t1; + + dg->evaluateResidual(); + + std::cout << "Output: {"; + for (double i : dg->getResidual()) + { + printf("%e ,", i); + } + std::cout << "}\n"; + + //Generated from matlab code with same parameters and inputs + std::vector true_vec{3.141592277589793e+02, + 8.941907747838389e-01, + 1.846733023014284e+00, + 3.141592277589793e+02, + 1.014543000000000e+02, + -1.507680000000000e+01, + 3.787993500000000e+02, + -1.300000000000000e+00, + 2.899095146477517e+02, + 2.939138495559215e+02, + 1.507210571826699e+07, + 1.659799832843673e+07, + -7.591593003913325e+03, + -8.376991073310774e+03, + 3.337988298081817e+03, + 2.684419146397466e+03}; + + std::cout << "Test the Relative Error\n"; + for (size_t i = 0; i < true_vec.size(); i++) + { + printf("%e ,\n", (true_vec[i] - dg->getResidual()[i]) / true_vec[i]); + } + + return 0; +} diff --git a/Examples/GenInfiniteBus/GenInfiniteBus.cpp b/Examples/GenInfiniteBus/GenInfiniteBus.cpp index d4305ae0..f366a5a8 100644 --- a/Examples/GenInfiniteBus/GenInfiniteBus.cpp +++ b/Examples/GenInfiniteBus/GenInfiniteBus.cpp @@ -93,6 +93,7 @@ int main() // allocate model components model->allocate(); + // Create numerical integrator and configure it for the generator model Ida* idas = new Ida(model); diff --git a/Examples/Grid3Bus/3bus.mat b/Examples/Grid3Bus/3bus.mat new file mode 100644 index 00000000..4663eb33 --- /dev/null +++ b/Examples/Grid3Bus/3bus.mat @@ -0,0 +1,47 @@ +( +function mpc = case5 +% Created by Reid Gomillion + +% MATPOWER + +%% MATPOWER Case Format : Version 2 +mpc.version = '2'; + +%%----- Power Flow Data -----%% +%% system MVA base +mpc.baseMVA = 100; + +%% bus data +% bus_i type Pd Qd Gs Bs area Vm Va baseKV zone Vmax Vmin +mpc.bus = [ + 1 3 2.0 0.0 0 0 0 1 0.0 0 0 0 0.0; + 2 1 2.5 -0.8 0 0 0 1 0.0 0 0 0 0.0; + 3 2 0 0 0 0 0 1.1 0.0 0 0 0 0.0; +]; + +%% generator data +% bus Pg Qg Qmax Qmin Vg mBase status Pmax Pmin Pc1 Pc2 Qc1min Qc1max Qc2min Qc2max ramp_agc ramp_10 ramp_30 ramp_q apf +mpc.gen = [ + 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; + 3 2.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; +]; + +%% branch data +% fbus tbus r x b rateA rateB rateC ratio angle status angmin angmax +mpc.branch = [ + 1 2 0 0.1 0 0 0 0 0 0 0 0 0; + 1 3 0 0.0666666 0 0 0 0 0 0 0 0 0; + 2 3 0 0.0833333 0 0 0 0 0 0 0 0 0; +]; + +%%----- OPF Data -----%% +%% generator cost data +% 1 startup shutdown n x1 y1 ... xn yn +% 2 startup shutdown n c(n-1) ... c0 +mpc.gencost = [ + 2 0 0 3 0 14 0; + 2 0 0 3 0 15 0; + 2 0 0 3 0 30 0; +]; + +) \ No newline at end of file diff --git a/Examples/Grid3Bus/Grid3BusSys.cpp b/Examples/Grid3Bus/Grid3BusSys.cpp index a379f401..1248f6aa 100644 --- a/Examples/Grid3Bus/Grid3BusSys.cpp +++ b/Examples/Grid3Bus/Grid3BusSys.cpp @@ -101,24 +101,24 @@ mpc.baseMVA = 100; %% bus data % bus_i type Pd Qd Gs Bs area Vm Va baseKV zone Vmax Vmin mpc.bus = [ - 1 3 2.0 0.0 0 0 0 1 0.0 0 0 0 0.0; - 2 1 2.5 -0.8 0 0 0 1 0.0 0 0 0 0.0; - 3 2 0 0 0 0 0 1.1 0.0 0 0 0 0.0; + 1 3 2.0 0.0 0 0 0 1.0 0.0 0 0 0 0.0; + 2 1 2.5 -0.8 0 0 0 1.0 0.0 0 0 0 0.0; + 3 2 0 0 0 0 0 1.1 0.0 0 0 0 0.0; ]; %% generator data % bus Pg Qg Qmax Qmin Vg mBase status Pmax Pmin Pc1 Pc2 Qc1min Qc1max Qc2min Qc2max ramp_agc ramp_10 ramp_30 ramp_q apf mpc.gen = [ - 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; - 3 2.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; + 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; + 3 2.0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; ]; %% branch data % fbus tbus r x b rateA rateB rateC ratio angle status angmin angmax mpc.branch = [ - 1 2 0 0.1 0 0 0 0 0 0 0 0 0; - 1 3 0 0.0666666 0 0 0 0 0 0 0 0 0; - 2 3 0 0.0833333 0 0 0 0 0 0 0 0 0; + 1 2 0 0.1 0 0 0 0 0 0 0 0 0; + 1 3 0 0.0666666 0 0 0 0 0 0 0 0 0; + 2 3 0 0.0833333 0 0 0 0 0 0 0 0 0; ]; %%----- OPF Data -----%% @@ -126,9 +126,9 @@ mpc.branch = [ % 1 startup shutdown n x1 y1 ... xn yn % 2 startup shutdown n c(n-1) ... c0 mpc.gencost = [ - 2 0 0 3 0 14 0; - 2 0 0 3 0 15 0; - 2 0 0 3 0 30 0; + 2 0 0 3 0 14 0; + 2 0 0 3 0 15 0; + 2 0 0 3 0 30 0; ]; )"; @@ -140,12 +140,15 @@ using namespace AnalysisManager; using namespace GridKit::Testing; using namespace GridKit::PowerSystemData; +constexpr double theta2_ref = -4.87979; // [deg] +constexpr double V2_ref = 1.08281; // [p.u.] +constexpr double theta3_ref = 1.46241; // [deg] /** * Testing the monlithic case via the class MiniGrid * @return returns 0 if pass o.w. fails */ -int monolithic_case() +int monolithicCase() { std::cout << "\nSolving power flow for a 3-bus monolithic model ...\n\n"; // Create a 3-bus model @@ -153,6 +156,7 @@ int monolithic_case() // allocate model model->allocate(); + model->initialize(); std::cout << "Model size: " << model->size() << "\n\n"; // Create numerical solver and attach the model to it. @@ -170,17 +174,17 @@ int monolithic_case() double V2 = model->V2(); double th3 = model->th3() * 180.0/M_PI; std::cout << "Solution:\n"; - std::cout << " theta2 = " << th2 << " deg, expected = " << " -4.87979 deg\n"; - std::cout << " V2 = " << V2 << " p.u., expected = " << " 1.08281 p.u.\n"; - std::cout << " theta3 = " << th3 << " deg, expected = " << " 1.46241 deg\n\n"; + std::cout << " theta2 = " << th2 << " deg, expected = " << theta2_ref << " deg\n"; + std::cout << " V2 = " << V2 << " p.u., expected = " << V2_ref << " p.u.\n"; + std::cout << " theta3 = " << th3 << " deg, expected = " << theta3_ref << " deg\n\n"; // Print solver performance statistics kinsol->printFinalStats(); int retval1 = 0; - retval1 += !isEqual(th2, -4.87979, 1e-4); - retval1 += !isEqual(V2, 1.08281, 1e-4); - retval1 += !isEqual(th3, 1.46241, 1e-4); + retval1 += !isEqual(th2, theta2_ref, 1e-4); + retval1 += !isEqual(V2, V2_ref , 1e-4); + retval1 += !isEqual(th3, theta3_ref, 1e-4); if(retval1 == 0) std::cout << "\nSuccess!\n\n\n"; @@ -197,7 +201,7 @@ int monolithic_case() * Run the Testing case for parser setup * @return returns 0 if pass o.w. fail */ -int parser_case() +int parserCase() { std::cout << "Solving same problem, but assembled from components via a parser ...\n\n"; @@ -213,6 +217,7 @@ int parser_case() // allocate model sysmodel->allocate(); + sysmodel->initialize(); std::cout << "Model size: " << sysmodel->size() << "\n\n"; // Create numerical solver and attach the model to it. @@ -232,17 +237,17 @@ int parser_case() std::cout << "Solution:\n"; - std::cout << " theta2 = " << th2 << " deg, expected = " << " -4.87979 deg\n"; - std::cout << " V2 = " << V2 << " p.u., expected = " << " 1.08281 p.u.\n"; - std::cout << " theta3 = " << th3 << " deg, expected = " << " 1.46241 deg\n\n"; + std::cout << " theta2 = " << th2 << " deg, expected = " << theta2_ref << " deg\n"; + std::cout << " V2 = " << V2 << " p.u., expected = " << V2_ref << " p.u.\n"; + std::cout << " theta3 = " << th3 << " deg, expected = " << theta3_ref << " deg\n\n"; // Print solver performance statistics kinsol->printFinalStats(); int retval2 = 0; - retval2 += !isEqual(th2, -4.87979, 1e-4); - retval2 += !isEqual(V2, 1.08281, 1e-4); - retval2 += !isEqual(th3, 1.46241, 1e-4); + retval2 += !isEqual(th2, theta2_ref, 1e-4); + retval2 += !isEqual(V2, V2_ref , 1e-4); + retval2 += !isEqual(th3, theta3_ref, 1e-4); if(retval2 == 0) std::cout << "\nSuccess!\n\n\n"; @@ -262,7 +267,7 @@ int parser_case() * Hardwired Test Case * @return 0 if pass otherwise fails */ -int hardwired_case() +int hardwiredCase() { std::cout << "Solving same problem, but assembled from components manually ...\n\n"; @@ -336,6 +341,7 @@ int hardwired_case() // allocate model sysmodel->allocate(); + sysmodel->initialize(); std::cout << "Model size: " << sysmodel->size() << "\n\n"; // Create numerical solver and attach the model to it. @@ -351,21 +357,21 @@ int hardwired_case() // Print solution double th2 = bus2->theta() * 180.0/M_PI; double V2 = bus2->V(); - double th3 = bus3->theta() * 180.0/M_PI; + double th3 = bus3->theta() * 180.0/M_PI; std::cout << "Solution:\n"; - std::cout << " theta2 = " << th2 << " deg, expected = " << " -4.87979 deg\n"; - std::cout << " V2 = " << V2 << " p.u., expected = " << " 1.08281 p.u.\n"; - std::cout << " theta3 = " << th3 << " deg, expected = " << " 1.46241 deg\n\n"; + std::cout << " theta2 = " << th2 << " deg, expected = " << theta2_ref << " deg\n"; + std::cout << " V2 = " << V2 << " p.u., expected = " << V2_ref << " p.u.\n"; + std::cout << " theta3 = " << th3 << " deg, expected = " << theta3_ref << " deg\n\n"; // Print solver performance statistics kinsol->printFinalStats(); int retval2 = 0; - retval2 += !isEqual(th2, -4.87979, 1e-4); - retval2 += !isEqual(V2, 1.08281, 1e-4); - retval2 += !isEqual(th3, 1.46241, 1e-4); + retval2 += !isEqual(th2, theta2_ref, 1e-4); + retval2 += !isEqual(V2, V2_ref , 1e-4); + retval2 += !isEqual(th3, theta3_ref, 1e-4); if(retval2 == 0) std::cout << "\nSuccess!\n\n\n"; @@ -385,10 +391,21 @@ int main() //return the results of each case int resolve = 0; std::cout << std::string(32,'-') << std::endl; - resolve += monolithic_case(); + resolve |= monolithicCase(); std::cout << std::string(32,'-') << std::endl; - resolve += hardwired_case(); + resolve |= parserCase(); std::cout << std::string(32,'-') << std::endl; - resolve += parser_case(); + resolve |= hardwiredCase(); + + if (resolve) + { + std::cout << "Failure!\n"; + } + else + { + std::cout << "Success!\n"; + } + + return resolve; } diff --git a/Examples/Microgrid/CMakeLists.txt b/Examples/Microgrid/CMakeLists.txt new file mode 100644 index 00000000..77c59d75 --- /dev/null +++ b/Examples/Microgrid/CMakeLists.txt @@ -0,0 +1,13 @@ + + + + +add_executable(microgrid Microgrid.cpp) +target_link_libraries(microgrid GRIDKIT::powerelec_disgen + GRIDKIT::powerelec_mircoline + GRIDKIT::powerelec_microload + GRIDKIT::solvers_dyn + GRIDKIT::powerelec_mircobusdq) + +add_test(NAME Microgrid COMMAND $) +install(TARGETS microgrid RUNTIME DESTINATION bin) diff --git a/Examples/Microgrid/Microgrid.cpp b/Examples/Microgrid/Microgrid.cpp new file mode 100644 index 00000000..30639437 --- /dev/null +++ b/Examples/Microgrid/Microgrid.cpp @@ -0,0 +1,447 @@ + + + +#include +#include +#include +#include +#include + + +#include +#include +#include +#include + +#include +#include +#include + + +int main(int argc, char const *argv[]) +{ + ///@todo Needs to be modified. Some components are small relative to others thus there error is high (or could be matlab vector issue) + double abstol = 1.0e-8; + double reltol = 1.0e-8; + size_t max_step_amount = 3000; + bool usejac = true; + + //Create model + ModelLib::PowerElectronicsModel* sysmodel = new ModelLib::PowerElectronicsModel(reltol, abstol, usejac, max_step_amount); + + //Modeled after the problem in the paper + double RN = 1.0e4; + + //DG Params + ModelLib::DistributedGeneratorParameters parms1; + parms1.wb_ = 2.0*M_PI*50.0; + parms1.wc_ = 31.41; + parms1.mp_ = 9.4e-5; + parms1.Vn_ = 380.0; + parms1.nq_ = 1.3e-3; + parms1.F_ = 0.75; + parms1.Kiv_ = 420.0; + parms1.Kpv_ = 0.1; + parms1.Kic_ = 2.0e4; + parms1.Kpc_ = 15.0; + parms1.Cf_ = 5.0e-5; + parms1.rLf_ = 0.1; + parms1.Lf_ = 1.35e-3; + parms1.rLc_ = 0.03; + parms1.Lc_ = 0.35e-3; + + ModelLib::DistributedGeneratorParameters parms2; + //Parameters from MATLAB Microgrid code for first DG + parms2.wb_ = 2.0*M_PI*50.0; + parms2.wc_ = 31.41; + parms2.mp_ = 12.5e-5; + parms2.Vn_ = 380.0; + parms2.nq_ = 1.5e-3; + parms2.F_ = 0.75; + parms2.Kiv_ = 390.0; + parms2.Kpv_ = 0.05; + parms2.Kic_ = 16.0e3; + parms2.Kpc_ = 10.5; + parms2.Cf_ = 50.0e-6; + parms2.rLf_ = 0.1; + parms2.Lf_ = 1.35e-3; + parms2.rLc_ = 0.03; + parms2.Lc_ = 0.35e-3; + + //Line params + 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); + + double rline3 = 0.23; + double Lline3 = 0.1 / (2.0 * M_PI * 50.0); + + //load parms + 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); + + + //indexing sets + size_t Nsize = 2; + // DGs + - refframe Lines + Loads + size_t vec_size_internals = 13*(2*Nsize) - 1 + (2 + 4*(Nsize - 1)) + 2*Nsize; + // \omegaref + BusDQ + size_t vec_size_externals = 1 + 2*(2*Nsize); + size_t dqbus1 = vec_size_internals + 1; + size_t dqbus2 = vec_size_internals + 3; + size_t dqbus3 = vec_size_internals + 5; + size_t dqbus4 = vec_size_internals + 7; + + size_t vec_size_total = vec_size_internals + vec_size_externals; + + + size_t indexv = 0; + + //dg 1 + ModelLib::DistributedGenerator *dg1 = new ModelLib::DistributedGenerator(0, parms1, true); + //ref motor + dg1->setExternalConnectionNodes(0,vec_size_internals); + //outputs + dg1->setExternalConnectionNodes(1,dqbus1); + dg1->setExternalConnectionNodes(2,dqbus1 + 1); + //"grounding" of the difference + dg1->setExternalConnectionNodes(3,-1); + //internal connections + for (size_t i = 0; i < 12; i++) + { + + dg1->setExternalConnectionNodes(4 + i,indexv + i); + } + indexv += 12; + sysmodel->addComponent(dg1); + + //dg 2 + ModelLib::DistributedGenerator *dg2 = new ModelLib::DistributedGenerator(1, parms1, false); + //ref motor + dg2->setExternalConnectionNodes(0,vec_size_internals); + //outputs + dg2->setExternalConnectionNodes(1,dqbus2); + dg2->setExternalConnectionNodes(2,dqbus2 + 1); + //internal connections + for (size_t i = 0; i < 13; i++) + { + + dg2->setExternalConnectionNodes(3 + i,indexv + i); + } + indexv += 13; + sysmodel->addComponent(dg2); + + + //dg 3 + ModelLib::DistributedGenerator *dg3 = new ModelLib::DistributedGenerator(2, parms2, false); + //ref motor + dg3->setExternalConnectionNodes(0,vec_size_internals); + //outputs + dg3->setExternalConnectionNodes(1,dqbus3); + dg3->setExternalConnectionNodes(2,dqbus3 + 1); + //internal connections + for (size_t i = 0; i < 13; i++) + { + + dg3->setExternalConnectionNodes(3 + i,indexv + i); + } + indexv += 13; + sysmodel->addComponent(dg3); + + + //dg 4 + ModelLib::DistributedGenerator *dg4 = new ModelLib::DistributedGenerator(3, parms2, false); + //ref motor + dg4->setExternalConnectionNodes(0,vec_size_internals); + //outputs + dg4->setExternalConnectionNodes(1,dqbus4); + dg4->setExternalConnectionNodes(2,dqbus4 + 1); + + //internal connections + for (size_t i = 0; i < 13; i++) + { + + dg4->setExternalConnectionNodes(3 + i,indexv + i); + } + indexv += 13; + sysmodel->addComponent(dg4); + + // Lines + + //line 1 + ModelLib::MicrogridLine *l1 = new ModelLib::MicrogridLine(4, rline1, Lline1); + //ref motor + l1->setExternalConnectionNodes(0,vec_size_internals); + //input connections + l1->setExternalConnectionNodes(1,dqbus1); + l1->setExternalConnectionNodes(2,dqbus1 + 1); + //output connections + l1->setExternalConnectionNodes(3,dqbus2); + l1->setExternalConnectionNodes(4,dqbus2 + 1); + //internal connections + for (size_t i = 0; i < 2; i++) + { + + l1->setExternalConnectionNodes(5 + i,indexv + i); + } + indexv += 2; + sysmodel->addComponent(l1); + + + //line 2 + ModelLib::MicrogridLine *l2 = new ModelLib::MicrogridLine(5, rline2, Lline2); + //ref motor + l2->setExternalConnectionNodes(0,vec_size_internals); + //input connections + l2->setExternalConnectionNodes(1,dqbus2); + l2->setExternalConnectionNodes(2,dqbus2 + 1); + //output connections + l2->setExternalConnectionNodes(3,dqbus3); + l2->setExternalConnectionNodes(4,dqbus3 + 1); + //internal connections + for (size_t i = 0; i < 2; i++) + { + + l2->setExternalConnectionNodes(5 + i,indexv + i); + } + indexv += 2; + sysmodel->addComponent(l2); + + //line 3 + ModelLib::MicrogridLine *l3 = new ModelLib::MicrogridLine(6, rline3, Lline3); + //ref motor + l3->setExternalConnectionNodes(0,vec_size_internals); + //input connections + l3->setExternalConnectionNodes(1,dqbus3); + l3->setExternalConnectionNodes(2,dqbus3 + 1); + //output connections + l3->setExternalConnectionNodes(3,dqbus4); + l3->setExternalConnectionNodes(4,dqbus4 + 1); + //internal connections + for (size_t i = 0; i < 2; i++) + { + + l3->setExternalConnectionNodes(5 + i,indexv + i); + } + indexv += 2; + sysmodel->addComponent(l3); + + // loads + + //load 1 + ModelLib::MicrogridLoad *load1 = new ModelLib::MicrogridLoad(7, rload1, Lload1); + //ref motor + load1->setExternalConnectionNodes(0,vec_size_internals); + //input connections + load1->setExternalConnectionNodes(1,dqbus1); + load1->setExternalConnectionNodes(2,dqbus1 + 1); + //internal connections + for (size_t i = 0; i < 2; i++) + { + + load1->setExternalConnectionNodes(3 + i,indexv + i); + } + indexv += 2; + sysmodel->addComponent(load1); + + //load 2 + ModelLib::MicrogridLoad *load2 = new ModelLib::MicrogridLoad(8, rload2, Lload2); + //ref motor + load2->setExternalConnectionNodes(0,vec_size_internals); + //input connections + load2->setExternalConnectionNodes(1,dqbus3); + load2->setExternalConnectionNodes(2,dqbus3 + 1); + //internal connections + for (size_t i = 0; i < 2; i++) + { + + load2->setExternalConnectionNodes(3 + i,indexv + i); + } + indexv += 2; + sysmodel->addComponent(load2); + + //Virtual PQ Buses + ModelLib::MicrogridBusDQ *bus1 = new ModelLib::MicrogridBusDQ(9, RN); + + bus1->setExternalConnectionNodes(0,dqbus1); + bus1->setExternalConnectionNodes(1,dqbus1 + 1); + sysmodel->addComponent(bus1); + + ModelLib::MicrogridBusDQ *bus2 = new ModelLib::MicrogridBusDQ(10, RN); + + bus2->setExternalConnectionNodes(0,dqbus2); + bus2->setExternalConnectionNodes(1,dqbus2 + 1); + sysmodel->addComponent(bus2); + + ModelLib::MicrogridBusDQ *bus3 = new ModelLib::MicrogridBusDQ(11, RN); + + bus3->setExternalConnectionNodes(0,dqbus3); + bus3->setExternalConnectionNodes(1,dqbus3 + 1); + sysmodel->addComponent(bus3); + + ModelLib::MicrogridBusDQ *bus4 = new ModelLib::MicrogridBusDQ(12, RN); + + bus4->setExternalConnectionNodes(0,dqbus4); + bus4->setExternalConnectionNodes(1,dqbus4 + 1); + sysmodel->addComponent(bus4); + + sysmodel->allocate(vec_size_total); + + std::cout << sysmodel->y().size() << std::endl; + std::cout << vec_size_internals << ", " << vec_size_externals << "\n"; + + //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; + } + + // Create Intial derivatives specifics generated in MATLAB + //DGs 1 + sysmodel->yp()[2] = parms1.Vn_; + sysmodel->yp()[4] = parms1.Kpv_ * parms1.Vn_; + sysmodel->yp()[6] = (parms1.Kpc_ * parms1.Kpv_ * parms1.Vn_) / parms1.Lf_; + sysmodel->yp()[12 + 3] = parms1.Vn_; + sysmodel->yp()[12 + 5] = parms1.Kpv_ * parms1.Vn_; + sysmodel->yp()[12 + 7] = (parms1.Kpc_ * parms1.Kpv_ * parms1.Vn_) / parms1.Lf_; + for (size_t i = 2; i < 4; i++) + { + sysmodel->yp()[13*i - 1 + 3] = parms2.Vn_; + sysmodel->yp()[13*i - 1 + 5] = parms2.Kpv_ * parms2.Vn_; + sysmodel->yp()[13*i - 1 + 7] = (parms2.Kpc_ * parms2.Kpv_ * parms2.Vn_) / parms2.Lf_; + } + + //since the intial P_com = 0 + sysmodel->y()[vec_size_internals] = parms1.wb_; + + + + sysmodel->initialize(); + sysmodel->evaluateResidual(); + + std::vector& fres = sysmodel->getResidual(); + std::cout << "Verify Intial Resisdual is Zero: {\n"; + for (size_t i = 0; i < fres.size(); i++) + { + printf("%lu : %e \n", i, fres[i]); + } + std::cout << "}\n"; + + sysmodel->updateTime(0.0, 1.0e-8); + sysmodel->evaluateJacobian(); + std::cout << "Intial Jacobian with alpha:\n"; + // std::cout << sysmodel->getJacobian().frobnorm() << "\n"; + + + //Create numerical integrator and configure it for the generator model + AnalysisManager::Sundials::Ida* idas = new AnalysisManager::Sundials::Ida(sysmodel); + + double t_init = 0.0; + double t_final = 1.0; + + // setup simulation + idas->configureSimulation(); + idas->getDefaultInitialCondition(); + idas->initializeSimulation(t_init); + + idas->runSimulation(t_final); + + std::vector& yfinial = sysmodel->y(); + + std::cout << "Final Vector y\n"; + for (size_t i = 0; i < yfinial.size(); i++) + { + std::cout << yfinial[i] << "\n"; + } + + //Generate from MATLAB code ODE form with tolerances of 1e-12 + std::vectortrue_vec{ + 2.297543153595780e+04, + 1.275311524125022e+04, + 3.763060183116022e-02, + -2.098153459325261e-02, + 1.848285659119097e-02, + -1.563291404944864e-04, + 6.321941907011718e+01, + -2.942264300846256e+01, + 3.634209302905854e+02, + -2.668928293656362e-06, + 6.321941919221522e+01, + -3.509200178595996e+01, + -7.555954467454730e-03, + 2.297580486511343e+04, + 8.742028429066131e+03, + 3.710079564796484e-02, + -1.421122598056797e-02, + 1.874079517807597e-02, + -9.891304812687215e-05, + 6.232933298360234e+01, + -1.796494061423331e+01, + 3.686353885026506e+02, + 3.465673854181523e-05, + 6.232933406188410e+01, + -2.371564475187742e+01, + -8.273939686941580e-02, + 1.727775042678524e+04, + 1.649365247247288e+04, + 3.116555157570849e-02, + -2.985990066758010e-02, + 2.250012115906506e-02, + -2.643873146501096e-04, + 4.861823510250247e+01, + -4.088592755441309e+01, + 3.552597163751238e+02, + -1.496407194199739e-04, + 4.861823504694532e+01, + -4.642797132602495e+01, + -8.445727984408551e-02, + 1.727723725566433e+04, + 9.182386962936238e+03, + 3.024959333190777e-02, + -1.617250828202081e-02, + 2.318056864131751e-02, + -1.295918667730514e-04, + 4.718938244522050e+01, + -1.935782085675469e+01, + 3.662262287803608e+02, + 1.076423957830039e-04, + 4.718938116520511e+01, + -2.507094256286497e+01, + -1.881248349415025e+01, + 2.114714832305742e+01, + 4.329946674909793e+01, + -3.037887936225145e+00, + -4.487023117352992e+01, + 2.895883729832657e+01, + 8.199613345691378e+01, + -5.623856502948122e+01, + 1.327498499660322e+02, + -8.228065162347022e+01, + 3.119995747945993e+02, + 3.576922945168803e+02, + -5.850795361581618e+00, + 3.641193316268954e+02, + -8.846325267612976e+00, + 3.472146752739036e+02, + -3.272400970143252e+01, + 3.604108939430972e+02, + -3.492842627398574e+01 + }; + + std::cout << "Test the Relative Error\n"; + for (size_t i = 0; i < true_vec.size(); i++) + { + printf("%lu : %e ,\n", i, abs(true_vec[i] - yfinial[i]) / abs(true_vec[i])); + } + + delete idas; + delete sysmodel; + + return 0; +} diff --git a/Examples/RLCircuit/CMakeLists.txt b/Examples/RLCircuit/CMakeLists.txt new file mode 100644 index 00000000..e7e05260 --- /dev/null +++ b/Examples/RLCircuit/CMakeLists.txt @@ -0,0 +1,13 @@ + + + + +add_executable(rlcircuit RLCircuit.cpp) +target_link_libraries(rlcircuit GRIDKIT::powerelec_capacitor + GRIDKIT::powerelec_inductor + GRIDKIT::powerelec_resistor + GRIDKIT::powerelec_voltagesource + GRIDKIT::solvers_dyn) + +add_test(NAME RLCircuit COMMAND $) +install(TARGETS rlcircuit RUNTIME DESTINATION bin) diff --git a/Examples/RLCircuit/RLCircuit.cpp b/Examples/RLCircuit/RLCircuit.cpp new file mode 100644 index 00000000..770307d8 --- /dev/null +++ b/Examples/RLCircuit/RLCircuit.cpp @@ -0,0 +1,146 @@ + + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + + +int main(int argc, char const *argv[]) +{ + double abstol = 1.0e-8; + double reltol = 1.0e-8; + bool usejac = true; + + //TODO:setup as named parameters + //Create circuit model + ModelLib::PowerElectronicsModel* sysmodel = new ModelLib::PowerElectronicsModel(reltol, abstol, usejac); + + size_t idoff = 0; + + //RL circuit parameters + double rinit = 1.0; + double linit = 1.0; + double vinit = 1.0; + + + //inductor + ModelLib::Inductor* induct = new ModelLib::Inductor(idoff,linit); + //Form index to node uid realations + // input + induct->setExternalConnectionNodes(0,1); + //output + induct->setExternalConnectionNodes(1,-1); + //internal + induct->setExternalConnectionNodes(2,2); + //add component + sysmodel->addComponent(induct); + + + //resistor + idoff++; + ModelLib::Resistor* resis = new ModelLib::Resistor(idoff, rinit); + //Form index to node uid realations + //input + resis->setExternalConnectionNodes(0,0); + //output + resis->setExternalConnectionNodes(1,1); + //add + sysmodel->addComponent(resis); + + //voltage source + idoff++; + ModelLib::VoltageSource* vsource = new ModelLib::VoltageSource(idoff, vinit); + //Form index to node uid realations + //input + vsource->setExternalConnectionNodes(0,-1); + //output + vsource->setExternalConnectionNodes(1,0); + //internal + vsource->setExternalConnectionNodes(2,3); + + + sysmodel->addComponent(vsource); + + sysmodel->allocate(4); + + std::cout << sysmodel->y().size() << std::endl; + + //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->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->initialize(); + sysmodel->evaluateResidual(); + + std::cout << "Verify Intial Resisdual is Zero: {"; + for (double i : sysmodel->getResidual()) + { + std::cout << i << ", "; + } + std::cout << "}\n"; + + + sysmodel->updateTime(0.0, 1.0); + sysmodel->evaluateJacobian(); + std::cout << "Intial Jacobian with alpha = 1:\n"; + sysmodel->getJacobian().printMatrix(); + + + // Create numerical integrator and configure it for the generator model + AnalysisManager::Sundials::Ida* idas = new AnalysisManager::Sundials::Ida(sysmodel); + + double t_init = 0.0; + double t_final = 1.0; + + // setup simulation + idas->configureSimulation(); + idas->getDefaultInitialCondition(); + idas->initializeSimulation(t_init); + + idas->runSimulation(t_final); + + std::vector& yfinial = sysmodel->y(); + + std::cout << "Final Vector y\n"; + for (size_t i = 0; i < yfinial.size(); i++) + { + std::cout << yfinial[i] << "\n"; + } + + 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]; + + std::cout << "Element-wise Relative error at t=" << t_final << "\n"; + for (size_t i = 0; i < yfinial.size(); i++) + { + std::cout << abs((yfinial[i] - yexact[i]) / yexact[i]) << "\n"; + } + + return 0; +} diff --git a/Examples/SparseTest/CMakeLists.txt b/Examples/SparseTest/CMakeLists.txt new file mode 100644 index 00000000..ef5c8ae0 --- /dev/null +++ b/Examples/SparseTest/CMakeLists.txt @@ -0,0 +1,7 @@ + + +add_executable(spmattest SparseTest.cpp) +target_link_libraries(spmattest GRIDKIT::SparseMatrix) + +add_test(NAME SparseMatrixTest COMMAND $) +install(TARGETS spmattest RUNTIME DESTINATION bin) diff --git a/Examples/SparseTest/SparseTest.cpp b/Examples/SparseTest/SparseTest.cpp new file mode 100644 index 00000000..c563dbd4 --- /dev/null +++ b/Examples/SparseTest/SparseTest.cpp @@ -0,0 +1,102 @@ + + +#include +#include +#include +#include +#include +#include +#include + +#include + +int main(int argc, char const *argv[]) +{ + std::vector val{0.1, 0.2, 0.3, 0.4}; + std::vector x{2,1,3,1}; + std::vector y{1,3,2,2}; + size_t n = 4; + size_t m = 4; + + COO_Matrix A = COO_Matrix(x,y,val,m,n); + + std::vector valn(4); + std::vector xn(4); + std::vector yn(4); + + std::tie(xn, yn, valn) = A.getEntries(); + + for (size_t i = 0; i < valn.size(); i++) + { + std::cout << valn[i] << "\n"; + } + + std::cout << "A:\n"; + A.printMatrix(); + + std::vector val2{0.5, 0.6, 0.7, 0.8, 1.0}; + std::vector x2{0,2,0,2,1}; + std::vector y2{3,3,2,2,3}; + COO_Matrix B = COO_Matrix(x2,y2,val2,m,n); + + std::cout << "B:\n"; + B.printMatrix(); + + A.axpy(2.0, B); + + std::cout << "A + 2B:\n"; + A.printMatrix(); + + std::vector r; + std::vector c; + std::vector v; + std::tie(r,c,v) = A.getDataToCSR(); + + for (size_t i = 0; i < r.size() - 1; i++) + { + std::cout << r[i] << std::endl; + size_t rdiff = r[i+1] - r[i]; + for (size_t j = 0; j < rdiff; j++) + { + std::cout << c[j + r[i]] << ", " << v[j + r[i]] << std::endl; + } + } + std::cout << r[r.size()-1] << std::endl; + + //Basic Verification test + std::vector rtest = {0, 2, 4, 7, 8}; + std::vector ctest = {2,3,2,3,1,2,3,2}; + std::vector valtest = {1.4, 1.0, 0.4, 2.2, 0.1, 1.6, 1.2, 0.3}; + + assert(rtest.size() == r.size()); + assert(ctest.size() == c.size()); + assert(valtest.size() == v.size()); + + int failval = 0; + for (size_t i = 0; i < rtest.size(); i++) + { + if (r[i] != rtest[i]) + { + failval--; + } + } + for (size_t i = 0; i < ctest.size(); i++) + { + double vdiff = v[i] - valtest[i]; + if (c[i] != ctest[i] || -1e-14 > vdiff || vdiff > 1e-14) + { + failval--; + } + } + + if (failval == 0) + { + std::cout << "Success!" << std::endl; + } + else + { + std::cout << "Failed!" << std::endl; + } + + return failval; +} diff --git a/ModelEvaluator.hpp b/ModelEvaluator.hpp index 5eda481d..ad56fbd4 100644 --- a/ModelEvaluator.hpp +++ b/ModelEvaluator.hpp @@ -62,6 +62,7 @@ #include #include +#include namespace ModelLib { @@ -88,15 +89,25 @@ namespace ModelLib virtual int initializeAdjoint() = 0; virtual int evaluateAdjointResidual() = 0; - //virtual int evaluateAdjointJacobian() = 0; + // virtual int evaluateAdjointJacobian() = 0; virtual int evaluateAdjointIntegrand() = 0; virtual IdxT size() = 0; virtual IdxT nnz() = 0; + + /** + * @brief Is the Jacobian defined. Used in IDA to determine wether DQ is used or not + * + * @return true + * @return false + */ + virtual bool hasJacobian() = 0; + virtual IdxT size_quad() = 0; virtual IdxT size_opt() = 0; virtual void updateTime(real_type t, real_type a) = 0; virtual void setTolerances(real_type& rtol, real_type& atol) const = 0; + virtual void setMaxSteps(IdxT& msa) const = 0; virtual std::vector& y() = 0; virtual const std::vector& y() const = 0; @@ -125,6 +136,10 @@ namespace ModelLib virtual std::vector& getResidual() = 0; virtual const std::vector& getResidual() const = 0; + + virtual COO_Matrix& getJacobian() = 0; + virtual const COO_Matrix& getJacobian() const = 0; + virtual std::vector& getIntegrand() = 0; virtual const std::vector& getIntegrand() const = 0; diff --git a/ModelEvaluatorImpl.hpp b/ModelEvaluatorImpl.hpp index a4ef43a0..0f8969f0 100644 --- a/ModelEvaluatorImpl.hpp +++ b/ModelEvaluatorImpl.hpp @@ -94,10 +94,12 @@ namespace ModelLib ypB_(size_), fB_(size_), gB_(size_opt_), + J_(COO_Matrix()), param_(size_opt_), param_up_(size_opt_), param_lo_(size_opt_) - {} + { + } virtual IdxT size() { @@ -109,6 +111,11 @@ namespace ModelLib return nnz_; } + virtual bool hasJacobian() + { + return false; + } + virtual IdxT size_quad() { return size_quad_; @@ -132,6 +139,11 @@ namespace ModelLib atol = atol_; } + virtual void setMaxSteps(IdxT& msa) const + { + msa = max_steps_; + } + std::vector& y() { return y_; @@ -222,6 +234,16 @@ namespace ModelLib return f_; } + COO_Matrix& getJacobian() + { + return J_; + } + + const COO_Matrix& getJacobian() const + { + return J_; + } + std::vector& getIntegrand() { return g_; @@ -252,6 +274,12 @@ namespace ModelLib return gB_; } + //@todo Fix ID naming + IdxT getIDcomponent() + { + return idc_; + } + protected: @@ -271,6 +299,8 @@ namespace ModelLib std::vector fB_; std::vector gB_; + COO_Matrix J_; + std::vector param_; std::vector param_up_; std::vector param_lo_; @@ -281,6 +311,10 @@ namespace ModelLib real_type rtol_; real_type atol_; + IdxT max_steps_; + + IdxT idc_; + }; diff --git a/PowerElectronicsModel.hpp b/PowerElectronicsModel.hpp new file mode 100644 index 00000000..ef109b92 --- /dev/null +++ b/PowerElectronicsModel.hpp @@ -0,0 +1,332 @@ + + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include + +namespace ModelLib +{ + + template + class PowerElectronicsModel : public ModelEvaluatorImpl + { + typedef CircuitComponent component_type; + + using ModelEvaluatorImpl::size_; + // using ModelEvaluatorImpl::size_quad_; + // using ModelEvaluatorImpl::size_opt_; + using ModelEvaluatorImpl::nnz_; + using ModelEvaluatorImpl::time_; + using ModelEvaluatorImpl::alpha_; + using ModelEvaluatorImpl::y_; + using ModelEvaluatorImpl::yp_; + // using ModelEvaluatorImpl::yB_; + // using ModelEvaluatorImpl::ypB_; + // using ModelEvaluatorImpl::tag_; + using ModelEvaluatorImpl::f_; + // using ModelEvaluatorImpl::fB_; + // using ModelEvaluatorImpl::g_; + // using ModelEvaluatorImpl::gB_; + using ModelEvaluatorImpl::J_; + using ModelEvaluatorImpl::rtol_; + using ModelEvaluatorImpl::atol_; + // using ModelEvaluatorImpl::param_; + // using ModelEvaluatorImpl::param_up_; + // using ModelEvaluatorImpl::param_lo_; + + public: + /** + * @brief Default constructor for the system model + */ + PowerElectronicsModel() : ModelEvaluatorImpl(0, 0, 0) + { + // Set system model tolerances as default + rtol_ = 1e-4; + atol_ = 1e-4; + this->max_steps_ = 2000; + // By default not use jacobian + usejac_ = false; + } + + /** + * @brief Constructor for the system model + */ + PowerElectronicsModel(double rt = 1e-4, double at = 1e-4, bool ju = false, IdxT msa = 2000) : ModelEvaluatorImpl(0, 0, 0) + { + // Set system model tolerances from input + rtol_ = rt; + atol_ = at; + this->max_steps_ = msa; + // Can choose if to use jacobain + usejac_ = ju; + } + + /** + * @brief Destructor for the system model + */ + virtual ~PowerElectronicsModel() + { + for (auto comp : this->components_) + delete comp; + } + + /** + * @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 + */ + int allocate() + { + + return 1; + } + + /** + * @brief Will check if each component has jacobian avalible. If one doesn't then jacobain is false. + * + * + * @return true + * @return false + */ + bool hasJacobian() + { + if (!this->usejac_) + return false; + + for (const auto &component : components_) + { + if (!component->hasJacobian()) + { + return false; + } + } + return true; + } + + /** + * @brief Allocate the vector data with size amount + * @todo Add capability to go through component model connection to get the size of the actual vector + * + * @param s + * @return int + */ + int allocate(IdxT s) + { + + // Allocate all components + size_ = s; + for (const auto &component : components_) + { + component->allocate(); + } + + // Allocate global vectors + y_.resize(size_); + yp_.resize(size_); + f_.resize(size_); + + return 0; + } + + /** + * @brief Set intial y and y' of each component + * + * @return int + */ + int initialize() + { + + // Initialize components + for (const auto &component : components_) + { + component->initialize(); + } + this->distributeVectors(); + + return 0; + } + + /** + * @brief Distribute y and y' to each component based of node connection graph + * + * @return int + */ + int distributeVectors() + { + for (const auto &component : components_) + { + for (IdxT j = 0; j < component->size(); ++j) + { + if (component->getNodeConnection(j) != neg1_) + { + component->y()[j] = y_[component->getNodeConnection(j)]; + component->yp()[j] = yp_[component->getNodeConnection(j)]; + } + else + { + component->y()[j] = 0.0; + component->yp()[j] = 0.0; + } + } + } + return 0; + } + + int tagDifferentiable() + { + return 0; + } + + /** + * @brief Evaluate Residuals at each component then collect them + * + * @return int + */ + int evaluateResidual() + { + for (IdxT i = 0; i < this->f_.size(); i++) + { + f_[i] = 0.0; + } + + this->distributeVectors(); + + // Update system residual vector + + for (const auto &component : components_) + { + // TODO:check return type + component->evaluateResidual(); + for (IdxT j = 0; j < component->size(); ++j) + { + //@todo should do a different grounding check + if (component->getNodeConnection(j) != neg1_) + { + f_[component->getNodeConnection(j)] += component->getResidual()[j]; + } + } + } + + return 0; + } + + /** + * @brief Creates the Sparse COO Jacobian representing \alpha dF/dy' + dF/dy + * + * @return int + */ + int evaluateJacobian() + { + J_.zeroMatrix(); + distributeVectors(); + + // Evaluate component jacs + for (const auto &component : components_) + { + component->evaluateJacobian(); + + // get references to local jacobain + std::tuple&, std::vector&, std::vector&> tpm = component->getJacobian().getEntries(); + const auto& [r, c, v] = tpm; + + // Create copies of data to handle groundings + std::vector rgr; + std::vector cgr; + std::vector vgr; + for (IdxT i = 0; i < static_cast(r.size()); i++) + { + if (component->getNodeConnection(r[i]) != neg1_ && component->getNodeConnection(c[i]) != neg1_) + { + rgr.push_back(component->getNodeConnection(r[i])); + cgr.push_back(component->getNodeConnection(c[i])); + vgr.push_back(v[i]); + } + } + + // AXPY to Global Jacobian + J_.axpy(1.0, rgr, cgr, vgr); + } + + return 0; + } + + /** + * @brief Evaluate integrands for the system quadratures. + */ + int evaluateIntegrand() + { + + return 0; + } + + /** + * @brief Initialize system adjoint. + * + * Updates variables and optimization parameters, then initializes + * adjoints locally and copies them to the system adjoint vector. + */ + int initializeAdjoint() + { + return 0; + } + + /** + * @brief Compute adjoint residual for the system model. + * + * + */ + int evaluateAdjointResidual() + { + return 0; + } + + /** + * @brief Evaluate adjoint integrand for the system model. + * + * + */ + int evaluateAdjointIntegrand() + { + return 0; + } + + /** + * @brief Distribute time and time scaling for each component + * + * @param t + * @param a + */ + void updateTime(ScalarT t, ScalarT a) + { + for (const auto &component : components_) + { + component->updateTime(t, a); + } + time_ = t; + alpha_ = a; + } + + void addComponent(component_type *component) + { + components_.push_back(component); + } + + private: + + static constexpr IdxT neg1_ = static_cast(-1); + + std::vector components_; + bool usejac_; + + }; // class PowerElectronicsModel + +} // namespace ModelLib diff --git a/Solver/Dynamic/CMakeLists.txt b/Solver/Dynamic/CMakeLists.txt index 75982996..ab38756c 100644 --- a/Solver/Dynamic/CMakeLists.txt +++ b/Solver/Dynamic/CMakeLists.txt @@ -66,6 +66,7 @@ gridkit_add_library(solvers_dyn LINK_LIBRARIES PUBLIC SUNDIALS::nvecserial PUBLIC SUNDIALS::idas + PUBLIC SUNDIALS::sunlinsolklu OUTPUT_NAME gridkit_solvers_dyn) diff --git a/Solver/Dynamic/Ida.cpp b/Solver/Dynamic/Ida.cpp index c71a8fd0..637d1106 100644 --- a/Solver/Dynamic/Ida.cpp +++ b/Solver/Dynamic/Ida.cpp @@ -61,10 +61,13 @@ #include #include -// #include #include /* access to IDADls interface */ #include +//Sundials Sparse KLU +#include +#include + #include "ModelEvaluator.hpp" #include "Ida.hpp" @@ -87,9 +90,35 @@ namespace Sundials tag_ = NULL; } + /** + * @brief Destroy the Ida< Scalar T, Idx T>:: Ida object + * + * @note if sysmodel is freed before this will fail. May want something agnostic to this + * + * @tparam ScalarT + * @tparam IdxT + */ template Ida::~Ida() { + N_VDestroy(yy_); + N_VDestroy(yp_); + N_VDestroy(yy0_); + N_VDestroy(yp0_); + if (model_->hasJacobian()) + { + SUNLinSolFree_KLU(linearSolver_); + SUNMatDestroy_Sparse(JacobianMat_); + } + else + { + SUNLinSolFree_Dense(linearSolver_); + SUNMatDestroy_Dense(JacobianMat_); + } + ///@todo this free is needed but on geninfbus this seg faults + // IDAFree(&solver_); + SUNContext_Free(&context_); + } template @@ -103,6 +132,9 @@ namespace Sundials yp_ = N_VClone(yy_); checkAllocation((void*) yp_, "N_VClone"); + //get intial conditions + this->getDefaultInitialCondition(); + // Create vectors to store restart initial condition yy0_ = N_VClone(yy_); checkAllocation((void*) yy0_, "N_VClone"); @@ -127,6 +159,13 @@ namespace Sundials model_->setTolerances(rtol, atol); ///< \todo Function name should be "getTolerances"! retval = IDASStolerances(solver_, rtol, atol); checkOutput(retval, "IDASStolerances"); + + IdxT msa; + model_->setMaxSteps(msa); + + /// \todo Need to set max number of steps based on user input! + retval = IDASetMaxNumSteps(solver_, msa); + checkOutput(retval, "IDASetMaxNumSteps"); // Tag differential variables std::vector& tag = model_->tag(); @@ -144,14 +183,7 @@ namespace Sundials } // Set up linear solver - JacobianMat_ = SUNDenseMatrix(model_->size(), model_->size(), context_); - checkAllocation((void*) JacobianMat_, "SUNDenseMatrix"); - - linearSolver_ = SUNLinSol_Dense(yy_, JacobianMat_, context_); - checkAllocation((void*) linearSolver_, "SUNLinSol_Dense"); - - retval = IDASetLinearSolver(solver_, linearSolver_, JacobianMat_); - checkOutput(retval, "IDASetLinearSolver"); + this->configureLinearSolver(); return retval; } @@ -160,16 +192,32 @@ namespace Sundials int Ida::configureLinearSolver() { int retval = 0; + if (model_->hasJacobian()) + { + JacobianMat_ = SUNSparseMatrix(model_->size(), model_->size(), model_->size() * model_->size(), CSR_MAT, context_); + checkAllocation((void*) JacobianMat_, "SUNSparseMatrix"); - // Set up linear solver - JacobianMat_ = SUNDenseMatrix(model_->size(), model_->size(), context_); - checkAllocation((void*) JacobianMat_, "SUNDenseMatrix"); + linearSolver_ = SUNLinSol_KLU(yy_, JacobianMat_, context_); + checkAllocation((void*) linearSolver_, "SUNLinSol_KLU"); + + retval = IDASetLinearSolver(solver_, linearSolver_, JacobianMat_); + checkOutput(retval, "IDASetLinearSolver"); - linearSolver_ = SUNLinSol_Dense(yy_, JacobianMat_, context_); - checkAllocation((void*) linearSolver_, "SUNLinSol_Dense"); + retval = IDASetJacFn(solver_, this->Jac); + checkOutput(retval, "IDASetJacFn"); + } + else + { + JacobianMat_ = SUNDenseMatrix(model_->size(), model_->size(), context_); + checkAllocation((void*) JacobianMat_, "SUNDenseMatrix"); + + linearSolver_ = SUNLinSol_Dense(yy_, JacobianMat_, context_); + checkAllocation((void*) linearSolver_, "SUNLinSol_Dense"); - retval = IDASetLinearSolver(solver_, linearSolver_, JacobianMat_); - checkOutput(retval, "IDASetLinearSolver"); + retval = IDASetLinearSolver(solver_, linearSolver_, JacobianMat_); + checkOutput(retval, "IDASetLinearSolver"); + + } return retval; } @@ -525,6 +573,47 @@ namespace Sundials return 0; } + template + int Ida::Jac(realtype t, realtype cj, N_Vector yy, N_Vector yp, N_Vector resvec, SUNMatrix J, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3) + { + + ModelLib::ModelEvaluator* model = static_cast*>(user_data); + + + model->updateTime(t, cj); + copyVec(yy, model->y()); + copyVec(yp, model->yp()); + + model->evaluateJacobian(); + COO_Matrix Jac = model->getJacobian(); + + //Get reference to the jacobian entries + std::tuple&, std::vector&, std::vector&> tpm = Jac.getEntries(); + const auto [r, c, val] = tpm; + + //get the CSR row pointers from COO matrix + std::vector csrrowdata = Jac.getCSRRowData(); + + SUNMatZero(J); + + //Set row pointers + sunindextype *rowptrs = SUNSparseMatrix_IndexPointers(J); + for (unsigned int i = 0; i < csrrowdata.size() ; i++) + { + rowptrs[i] = csrrowdata[i]; + } + + sunindextype *colvals = SUNSparseMatrix_IndexValues(J); + realtype *data = SUNSparseMatrix_Data(J); + //Copy data from model jac to sundials + for (unsigned int i = 0; i < c.size(); i++ ) + { + colvals[i] = c[i]; + data[i] = val[i]; + } + + return 0; + } template int Ida::Integrand(realtype tt, N_Vector yy, N_Vector yp, N_Vector rhsQ, void *user_data) diff --git a/Solver/Dynamic/Ida.hpp b/Solver/Dynamic/Ida.hpp index 3fa3368c..a4cea956 100644 --- a/Solver/Dynamic/Ida.hpp +++ b/Solver/Dynamic/Ida.hpp @@ -65,7 +65,7 @@ #include #include #include /* access to sparse SUNMatrix */ -// #include /* access to KLU linear solver */ +#include /* access to KLU linear solver */ #include /* access to dense linear solver */ #include "ModelEvaluator.hpp" @@ -165,6 +165,11 @@ namespace AnalysisManager static int Residual(realtype t, N_Vector yy, N_Vector yp, N_Vector rr, void *user_data); + + static int Jac(realtype t, realtype cj, + N_Vector yy, N_Vector yp, N_Vector resvec, + SUNMatrix J, void *user_data, + N_Vector tmp1, N_Vector tmp2, N_Vector tmp3); static int Integrand(realtype t, N_Vector yy, N_Vector yp, diff --git a/SparseMatrix/CMakeLists.txt b/SparseMatrix/CMakeLists.txt new file mode 100644 index 00000000..b51d669e --- /dev/null +++ b/SparseMatrix/CMakeLists.txt @@ -0,0 +1,7 @@ + +add_library(SparseMatrix INTERFACE) +include_directories(SparseMatrix INTERFACE ${CMAKE_CURRENT_LIST_DIR}) + +add_library(GRIDKIT::SparseMatrix ALIAS SparseMatrix) + + diff --git a/SparseMatrix/COO_Matrix.hpp b/SparseMatrix/COO_Matrix.hpp new file mode 100644 index 00000000..9036085a --- /dev/null +++ b/SparseMatrix/COO_Matrix.hpp @@ -0,0 +1,766 @@ + + +#include +#include +#include +#include +#include +#include +#include +#include + +/** + * @brief Quick class to provide sparse matrices of COO type. Simplifies data movement + * + * @todo add functionality to keep track of multiple sorted_ list. Faster adding of new entries and will have a threshold to sort completely. + * + * m x n sparse matrix + */ +template +class COO_Matrix +{ +private: + std::vector values_; + std::vector row_indexes_; + std::vector column_indexes_; + Intdx rows_size_; + Intdx columns_size_; + bool sorted_; +public: + //Constructors + COO_Matrix(std::vector r, std::vector c, std::vector v, Intdx m, Intdx n); + COO_Matrix(Intdx m, Intdx n); + COO_Matrix(); + ~COO_Matrix(); + + + //Operations + + // --- Functions which call sort --- + std::tuple, std::vector> getRowCopy(Intdx r); + std::tuple&, std::vector&, std::vector&> getEntries(); + std::tuple, std::vector, std::vector> getEntrieCopies(); + std::tuple, std::vector, std::vector> getEntrieCopiesSubMatrix(std::vector submap); + + std::tuple, std::vector, std::vector> getDataToCSR(); + std::vector getCSRRowData(); + + // BLAS. Will sort before running + void setValues(std::vector r, std::vector c, std::vector v); + void axpy(ScalarT alpha, COO_Matrix& a); + void axpy(ScalarT alpha, std::vector r, std::vector c, std::vector v); + void scal(ScalarT alpha); + ScalarT frobnorm(); + + // --- Permutation Operations --- + //No sorting is actually done. Only done when nesscary + void permutation(std::vector row_perm, std::vector col_perm); + void permutationSizeMap(std::vector row_perm, std::vector col_perm, Intdx m, Intdx n); + + void zeroMatrix(); + + void identityMatrix(Intdx n); + + //Resort values_ + void sortSparse(); + bool isSorted(); + Intdx nnz(); + + std::tuple getDimensions(); + + void printMatrix(); + + + static void sortSparseCOO(std::vector &rows, std::vector &columns, std::vector &vals); + +private: + Intdx indexStartRow(const std::vector &rows, Intdx r); + Intdx sparseCordBinarySearch(const std::vector &rows, const std::vector &columns, Intdx ri, Intdx ci); + bool checkIncreaseSize(Intdx r, Intdx c); + +}; + +/** + * @brief Get copy of row values_ + * + * @tparam ScalarT + * @tparam Intdx + * @param r + * @return std::tuple, std::vector> + */ +template +inline std::tuple, std::vector> COO_Matrix::getRowCopy(Intdx r) +{ + if (!this->sorted_) + { + this->sortSparse(); + } + Intdx rowindex = this->indexStartRow(r); + + + if (rowindex == -1) + { + return {std::vector(),std::vector()}; + } + + Intdx rsize = rowindex; + do + { + rsize++; + } while (rsize < this->values_.size() && this->row_indexes_[rsize] == r); + + return {{this->column_indexes_.begin() + rowindex, this->column_indexes_.begin() + rsize},{this->values_.begin() + rowindex, this->values_.begin() + rsize}}; +} + +/** + * @brief Get all Entries pointers. Will sort before returnings + * + * @tparam ScalarT + * @tparam Intdx + * @return std::tuple, std::vector, std::vector> + */ +template +inline std::tuple&, std::vector&, std::vector&> COO_Matrix::getEntries() +{ + if (!this->sorted_) + { + this->sortSparse(); + } + return {this->row_indexes_, this->column_indexes_, this->values_}; +} + +/** + * @brief Get copies of the data. Sorted_ before returning + * + * @tparam ScalarT + * @tparam Intdx + * @return std::tuple, std::vector, std::vector> + */ +template +inline std::tuple, std::vector, std::vector> COO_Matrix::getEntrieCopies() +{ + if (!this->sorted_) + { + this->sortSparse(); + } + return {this->row_indexes_, this->column_indexes_, this->values_}; +} + +/** + * @brief Returns the data into CSR Format + * + * @tparam ScalarT + * @tparam Intdx + * @return std::tuple, std::vector, std::vector> + */ +template +inline std::tuple, std::vector, std::vector> COO_Matrix::getDataToCSR() +{ + if (!this->isSorted()) this->sortSparse(); + std::vector rowsizevec(this->rows_size_ + 1, 0); + Intdx counter = 0; + for (Intdx i = 0; i < static_cast(rowsizevec.size() - 1); i++) + { + rowsizevec[i + 1] = rowsizevec[i]; + while (counter < static_cast(this->row_indexes_.size()) && i == this->row_indexes_[counter]) + { + rowsizevec[i+1]++; + counter++; + } + } + return {rowsizevec, this->column_indexes_, this->values_}; +} + +/** + * @brief Only creates the row data + * + * @todo swap this with having the matrix store the data and updates. This can then be passed by reference + * + * + * @tparam ScalarT + * @tparam Intdx + * @return std::vector + */ +template +inline std::vector COO_Matrix::getCSRRowData() +{ + if (!this->isSorted()) this->sortSparse(); + std::vector rowsizevec(this->rows_size_ + 1, 0); + Intdx counter = 0; + for (Intdx i = 0; i < static_cast(rowsizevec.size() - 1); i++) + { + rowsizevec[i + 1] = rowsizevec[i]; + while (counter < static_cast(this->row_indexes_.size()) && i == this->row_indexes_[counter]) + { + rowsizevec[i+1]++; + counter++; + } + } + return rowsizevec; +} + +/** + * @brief Given set of vector data it will set the values_ into the matrix + * + * @tparam ScalarT + * @tparam Intdx + * @param r + * @param c + * @param v + */ +template +inline void COO_Matrix::setValues(std::vector r, std::vector c, std::vector v) +{ + //sort input + this->sortSparseCOO(r, c, v); + + + //Duplicated with axpy. Could replace with function depdent on lambda expression + Intdx aiter = 0; + //iterate for all current values_ in matrix + for (Intdx i = 0; i < static_cast(this->row_indexes_.size()); i++) + { + //pushback values_ when they are not in current matrix + while(aiter < static_cast(r.size()) && (r[aiter] < this->row_indexes_[i] || (r[aiter] == this->row_indexes_[i] && c[aiter] < this->column_indexes_[i]))) + { + this->row_indexes_.push_back(r[aiter]); + this->column_indexes_.push_back(c[aiter]); + this->values_.push_back(v[aiter]); + this->checkIncreaseSize(r[aiter], c[aiter]); + aiter++; + } + if (aiter >= static_cast(r.size())) + { + break; + } + + + if (r[aiter] == this->row_indexes_[i] && c[aiter] == this->column_indexes_[i]) + { + this->values_[i] = v[aiter]; + aiter++; + } + } + //push back rest that was not found sorted + for (Intdx i = aiter; i < static_cast(r.size()); i++) + { + this->row_indexes_.push_back(r[i]); + this->column_indexes_.push_back(c[i]); + this->values_.push_back(v[i]); + + this->checkIncreaseSize(r[i], c[i]); + } + + this->sorted_ = false; + +} + +/** + * @brief BLAS axpy operation on another COO matrix. Will sort both matrices before acting + * + * @tparam ScalarT + * @tparam Intdx + * @param alpha + * @param a + */ +template +inline void COO_Matrix::axpy(ScalarT alpha, COO_Matrix& a) +{ + if (alpha == 0) + { + return; + } + + if (!this->sorted_) + { + this->sortSparse(); + } + if (!a.isSorted()) + { + a.sortSparse(); + } + Intdx m = 0; + Intdx n = 0; + std::tuple&, std::vector&, std::vector&> tpm = a.getEntries(); + const auto& [r, c, val] = tpm; + std::tie(m,n) = a.getDimensions(); + + //Increase size as nesscary + this->rows_size_ = this->rows_size_ > m ? this->rows_size_ : m; + this->columns_size_ = this->columns_size_ > n ? this->columns_size_ : n; + + Intdx aiter = 0; + //iterate for all current values_ in matrix + for (Intdx i = 0; i < static_cast(this->row_indexes_.size()); i++) + { + //pushback values_ when they are not in current matrix + while(aiter < static_cast(r.size()) && (r[aiter] < this->row_indexes_[i] || (r[aiter] == this->row_indexes_[i] && c[aiter] < this->column_indexes_[i]))) + { + this->row_indexes_.push_back(r[aiter]); + this->column_indexes_.push_back(c[aiter]); + this->values_.push_back(alpha * val[aiter]); + + this->checkIncreaseSize(r[aiter], c[aiter]); + aiter++; + } + if (aiter >= static_cast(r.size())) + { + break; + } + + + if (r[aiter] == this->row_indexes_[i] && c[aiter] == this->column_indexes_[i]) + { + this->values_[i] += alpha * val[aiter]; + aiter++; + } + } + //push back rest that was not found sorted_ + for (Intdx i = aiter; i < static_cast(r.size()); i++) + { + this->row_indexes_.push_back(r[i]); + this->column_indexes_.push_back(c[i]); + this->values_.push_back(alpha * val[i]); + + this->checkIncreaseSize(r[i], c[i]); + } + + this->sorted_ = false; +} + +/** + * @brief axpy on 3list. + * + * @tparam ScalarT + * @tparam Intdx + * @param alpha + * @param r + * @param c + * @param v + */ +template +inline void COO_Matrix::axpy(ScalarT alpha, std::vector r, std::vector c, std::vector v) +{ + if (alpha == 0) return; + + if (!this->sorted_) + { + this->sortSparse(); + } + + //sort input + this->sortSparseCOO(r, c, v); + + Intdx aiter = 0; + //iterate for all current values_ in matrix + for (Intdx i = 0; i < static_cast(this->row_indexes_.size()); i++) + { + //pushback values_ when they are not in current matrix + while(aiter < static_cast(r.size()) && (r[aiter] < this->row_indexes_[i] || (r[aiter] == this->row_indexes_[i] && c[aiter] < this->column_indexes_[i]))) + { + this->row_indexes_.push_back(r[aiter]); + this->column_indexes_.push_back(c[aiter]); + this->values_.push_back(alpha * v[aiter]); + + this->checkIncreaseSize(r[aiter], c[aiter]); + aiter++; + } + if (aiter >= static_cast(r.size())) + { + break; + } + + + if (r[aiter] == this->row_indexes_[i] && c[aiter] == this->column_indexes_[i]) + { + this->values_[i] += alpha * v[aiter]; + aiter++; + } + } + //push back rest that was not found sorted_ + for (Intdx i = aiter; i < static_cast(r.size()); i++) + { + this->row_indexes_.push_back(r[i]); + this->column_indexes_.push_back(c[i]); + this->values_.push_back(alpha * v[i]); + + this->checkIncreaseSize(r[i], c[i]); + } + + this->sorted_ = false; +} + +/** + * @brief Scale all values_ + * + * @tparam ScalarT + * @tparam Intdx + * @param alpha + */ +template +inline void COO_Matrix::scal(ScalarT alpha) +{ + for (auto i = this->values_.begin(); i < this->values_.end(); i++) + { + *i *= alpha; + } +} + +/** + * @brief Frobenius Norm of the Matrix + * + * @tparam ScalarT + * @tparam Intdx + * @return ScalarT + */ +template +inline ScalarT COO_Matrix::frobnorm() +{ + ScalarT totsum = 0.0; + for (auto i = this->values_.begin(); i < this->values_.end(); i++) + { + totsum += abs(*i)^2; + } + return totsum; +} + +/** + * @brief Permutate the matrix to a different one. Only changes the coordinates + * + * @tparam ScalarT + * @tparam Intdx + * @param row_perm + * @param col_perm + */ +template +inline void COO_Matrix::permutation(std::vector row_perm, std::vector col_perm) +{ + assert(row_perm.size() = this->rows_size_); + assert(col_perm.size() = this->columns_size_); + + for (int i = 0; i < this->values_.size(); i++) + { + this->row_indexes_[i] = row_perm[this->row_indexes_[i]]; + this->column_indexes_[i] = col_perm[this->column_indexes_[i]]; + } + this->sorted_ = false; + //cycle sorting maybe useful since permutations are already known +} + +/** + * @brief Permutates the matrix and can change its size efficently + * if size is shrinking and value is to be removed the negative one + * + * @tparam ScalarT + * @tparam Intdx + * @param row_perm size of m + * @param col_perm size of n + * @param m + * @param n + */ +template +inline void COO_Matrix::permutationSizeMap(std::vector row_perm, std::vector col_perm, Intdx m, Intdx n) +{ + assert(row_perm.size() == this->rows_size_); + assert(col_perm.size() == this->columns_size_); + + this->rows_size_ = m; + this->columns_size_ = n; + + for (int i = 0; i < this->values_.size(); i++) + { + if (row_perm[this->row_indexes_[i]] == -1 || col_perm[this->column_indexes_[i]] == -1) + { + this->values_[i] = 0; + } + else + { + this->row_indexes_[i] = row_perm[this->row_indexes_[i]]; + this->column_indexes_[i] = col_perm[this->column_indexes_[i]]; + } + } + this->sorted_ = false; +} + +/** + * @brief Turn matrix into the zero matrix. Does not actual delete memory + * + * @tparam ScalarT + * @tparam Intdx + */ +template +inline void COO_Matrix::zeroMatrix() +{ + //resize doesn't effect capacity if smaller + this->column_indexes_.resize(0); + this->row_indexes_.resize(0); + this->values_.resize(0); + this->sorted_ = true; +} + +template +inline void COO_Matrix::identityMatrix(Intdx n) +{ + //Reset Matrix + this->zeroMatrix(); + for (Intdx i = 0; i < n; i++) + { + this->column_indexes_[i] = i; + this->row_indexes_[i] = i; + this->values_[i] = 1.0; + } + this->sorted_ = true; +} + +/** + * @brief Restructure the sparse matrix for faster accesses and modifications + * + * @tparam ScalarT + * @tparam Intdx + */ +template +inline void COO_Matrix::sortSparse() +{ + this->sortSparseCOO(this->row_indexes_, this->column_indexes_, this->values_); + this->sorted_ = true; +} + +template +inline bool COO_Matrix::isSorted() +{ + return this->sorted_; +} + +template +inline Intdx COO_Matrix::nnz() +{ + return static_cast(this->values_.size); +} + +template +inline std::tuple COO_Matrix::getDimensions() +{ + return std::tuple(this->rows_size_, this->columns_size_); +} + +/** + * @brief Print matrix that is sorted_ + * + * @tparam ScalarT + * @tparam Intdx + */ +template +inline void COO_Matrix::printMatrix() +{ + if (this->sorted_ == false) + { + this->sortSparse(); + } + + std::cout << "Sparse COO Matrix\n"; + std::cout << "(x , y, value)\n"; + for (size_t i = 0; i < this->values_.size(); i++) + { + std::cout << "(" << this->row_indexes_[i] + << ", " << this->column_indexes_[i] + << ", " << this->values_[i] << ")\n"; + } + std::cout << std::flush; +} + +/** + * @brief Find the lowest row cordinate from set of provided cordinates + * + * Assumes rows and columns are sorted_ + * @tparam ScalarT + * @tparam Intdx + * @param r + * @return Intdx + */ +template +inline Intdx COO_Matrix::indexStartRow(const std::vector &rows, Intdx r) +{ + //Specialized Binary Search for Lowest Row + Intdx i1 = 0; + Intdx i2 = rows->size()-1; + Intdx m_smallest = -1; + Intdx m = -1; + while (i1 <= i2) + { + m = (i2 + i1) / 2; + //rows + if (rows[m] < r) + { + i1 = m + 1; + } + else if (r < rows[m]) + { + i2 = m - 1; + } + else + { + if (i1 == i2) + { + return m_smallest; + } + + //Keep track of smallest cordinate + m_smallest = m; + i2 = m - 1; + } + } + return m_smallest; +} + +/** + * @brief Basic binary search + * + * @tparam ScalarT + * @tparam Intdx + * @param rows + * @param columns + * @param ri + * @param ci + * @return Intdx + */ +template +inline Intdx COO_Matrix::sparseCordBinarySearch(const std::vector &rows, const std::vector &columns, Intdx ri, Intdx ci) +{ + assert(rows.size() == columns.size()); + //basic binary search + Intdx i1 = 0; + Intdx i2 = rows.size()-1; + Intdx m = 0; + while (i1 <= i2) + { + m = (i2 + i1) / 2; + //rows + if (rows[m] < ri) + { + i1 = m + 1; + } + else if (ri < rows[m]) + { + i2 = m - 1; + } + else + { + if (columns[m] < ci) + { + i1 = m + 1; + } + else if (ci < columns[m]) + { + i2 = m - 1; + } + break; + } + } + + return m; +} + +template +inline bool COO_Matrix::checkIncreaseSize(Intdx r, Intdx c) +{ + bool changed = false; + if (r + 1 > this->rows_size_) + { + this->rows_size_ = r + 1; + changed = true; + } + if (c + 1 > this->columns_size_) + { + this->columns_size_ = c + 1; + changed = true; + } + + return changed; +} + +/** + * @brief Sort a disoreded set of values_. Assume nothing on order. + * + * @todo simple setup. Should add stable sorting since list are pre-sorted_ + * + * @tparam ScalarT + * @tparam Intdx + * @param rows + * @param columns + * @param values_ + */ +template +inline void COO_Matrix::sortSparseCOO(std::vector &rows, std::vector &columns, std::vector &vals) +{ + + //index based sort code + // https://stackoverflow.com/questions/25921706/creating-a-vector-of-indices-of-a-sorted_-vector + //cannot call sort since two arrays are used instead + std::vector ordervec(rows.size()); + std::size_t n(0); + std::generate(std::begin(ordervec), std::end(ordervec), [&]{ return n++; }); + + //Sort by row first then column. + std::sort( std::begin(ordervec), + std::end(ordervec), + [&](int i1, int i2) { return (rows[i1] < rows[i2]) || + (rows[i1] == rows[i2] && columns[i1] < columns[i2]); } ); + + + //reorder based of index-sorting. Only swap cost no extra memory. + // @todo see if extra memory creation is fine + // https://stackoverflow.com/a/22183350 + for (size_t i = 0; i < ordervec.size(); i++) + { + //permutation swap + while (ordervec[i] != ordervec[ordervec[i]]) + { + std::swap(rows[ordervec[i]], rows[ordervec[ordervec[i]]]); + std::swap(columns[ordervec[i]], columns[ordervec[ordervec[i]]]); + std::swap(vals[ordervec[i]], vals[ordervec[ordervec[i]]]); + + //swap orderings + std::swap(ordervec[i], ordervec[ordervec[i]]); + } + + } +} + +template +inline COO_Matrix::COO_Matrix(std::vector r, std::vector c, std::vector v, Intdx m, Intdx n) +{ + this->values_ = v; + this->row_indexes_ = r; + this->column_indexes_ = c; + this->rows_size_ = m; + this->columns_size_ = n; + this->sorted_ = false; +} + +template +inline COO_Matrix::COO_Matrix(Intdx m, Intdx n) +{ + this->rows_size_ = m; + this->columns_size_ = n; + this->values_ = std::vector(); + this->row_indexes_ = std::vector(); + this->column_indexes_ = std::vector(); + this->sorted_ = false; +} + +template +inline COO_Matrix::COO_Matrix() +{ + this->rows_size_ = 0; + this->columns_size_ = 0; + this->values_ = std::vector(); + this->row_indexes_ = std::vector(); + this->column_indexes_ = std::vector(); + this->sorted_ = false; +} + +template +COO_Matrix::~COO_Matrix() +{ + +} diff --git a/SystemModel.hpp b/SystemModel.hpp index 34c7ff12..14b7e4e9 100644 --- a/SystemModel.hpp +++ b/SystemModel.hpp @@ -117,6 +117,7 @@ class SystemModel : public ModelEvaluatorImpl // Set system model tolerances rtol_ = 1e-7; atol_ = 1e-9; + this->max_steps_=2000; } /** @@ -183,6 +184,17 @@ class SystemModel : public ModelEvaluatorImpl return 0; } + /** + * @brief Assume that jacobian is not avalible + * + * @return true + * @return false + */ + bool hasJacobian() + { + return false; + } + /** * @brief Initialize buses first, then all the other components. * diff --git a/SystemSteadyStateModel.hpp b/SystemSteadyStateModel.hpp index c32dd5ec..ca1af16f 100644 --- a/SystemSteadyStateModel.hpp +++ b/SystemSteadyStateModel.hpp @@ -426,4 +426,4 @@ class SystemSteadyStateModel : public ModelEvaluatorImpl }; // class SystemSteadyStateModel -} // namespace ModelLib +} // namespace ModelLib \ No newline at end of file