diff --git a/CMake/FindSuiteSparse.cmake b/CMake/FindSuiteSparse.cmake index 2d73feb76..5a1026935 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 19048d375..7a6e651e0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -81,11 +81,14 @@ 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) +#Provide string to source +add_compile_definitions("SOURCE_ROOT=${CMAKE_SOURCE_DIR}") + list(APPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/CMake) # https://gitlab.kitware.com/cmake/community/-/wikis/doc/cmake/RPATH-handling#always-full-rpath @@ -110,7 +113,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 +146,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 000000000..e54f5a150 --- /dev/null +++ b/CircuitGraph.hpp @@ -0,0 +1,119 @@ + + +#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 should replace N and E with Node and Component classes respectively + * + * @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/BusFactory.hpp b/ComponentLib/Bus/BusFactory.hpp index 68fc5578e..232171d49 100644 --- a/ComponentLib/Bus/BusFactory.hpp +++ b/ComponentLib/Bus/BusFactory.hpp @@ -97,4 +97,4 @@ namespace ModelLib { } }; -} // namespace ModelLib \ No newline at end of file +} // namespace ModelLib diff --git a/ComponentLib/Bus/BusPV.cpp b/ComponentLib/Bus/BusPV.cpp index d5a05dc98..13ff37512 100644 --- a/ComponentLib/Bus/BusPV.cpp +++ b/ComponentLib/Bus/BusPV.cpp @@ -95,7 +95,7 @@ template 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; diff --git a/ComponentLib/CMakeLists.txt b/ComponentLib/CMakeLists.txt index 0c422d5b9..29a5a7988 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 c34df84ad..5e3be2041 100644 --- a/ComponentLib/Generator/GeneratorPV.hpp +++ b/ComponentLib/Generator/GeneratorPV.hpp @@ -124,18 +124,14 @@ 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 000000000..3fba4be7f --- /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(DiscreteGenerator) +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 000000000..62efc1795 --- /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 000000000..b6287b359 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/Capacitor/Capacitor.cpp @@ -0,0 +1,137 @@ + + + +#include +#include +#include +#include "Capacitor.hpp" + +namespace ModelLib { + +/*! + * @brief Constructor for a constant load model + * + * Calls default ModelEvaluatorImpl constructor. + */ + +template +Capacitor::Capacitor(IdxT id, ScalarT C) + : C_(C) +{ + this->size_ = 3; + this->n_intern = 1; + this->n_extern = 2; + this->extern_indices = {0,1}; + this->idc_ = id; +} + +template +Capacitor::~Capacitor() +{ +} + +/*! + * @brief allocate method computes sparsity pattern of the Jacobian. + */ +template +int Capacitor::allocate() +{ + this->y_.resize(this->size_); + this->yp_.resize(this->size_); + this->f_.resize(this->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 Contributes to the bus residual. + * + * Must be connected to a PQ bus. + */ +template +int Capacitor::evaluateResidual() +{ + //input + this->f_[0] = C_ * this->yp_[2]; + //output + this->f_[1] = -C_ * this->yp_[2]; + + //internal + this->f_[2] = -C_ * this->yp_[2] + this->y_[0] - this->y_[1] - this->y_[2]; + return 0; +} + +template +int Capacitor::evaluateJacobian() +{ + this->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}; + this->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_, -this->C_}; + COO_Matrix Jacder = COO_Matrix(rcordder, ccordder, valsder,3,3); + + //Perform dF/dy + \alpha dF/dy' + this->J_.AXPY(this->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 000000000..e1a9974c3 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/Capacitor/Capacitor.hpp @@ -0,0 +1,65 @@ + + +#ifndef _CAP_HPP_ +#define _CAP_HPP_ + +#include +#include +#include + + +namespace ModelLib +{ + template class BaseBus; +} + + +namespace ModelLib +{ + /*! + * @brief Declaration of a passive 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_; + + + 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 000000000..fe2efdf96 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/CircuitComponent.hpp @@ -0,0 +1,139 @@ + + +#ifndef _CIRCCOMP_HPP_ +#define _CIRCCOMP_HPP_ + +#include +#include +#include +#include + + +namespace ModelLib +{ + /*! + * @brief Declaration of a passive 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 this->n_extern; + } + + size_t getInternalSize() + { + return this->n_intern; + } + + std::set getExternIndices() + { + return this->extern_indices; + } + + bool setExternalConnectionNodes(size_t index, IdxT id) + { + this->connection_nodes[index] = id; + return true; + } + + IdxT getNodeConnection(size_t index) + { + return this->connection_nodes.at(index); + } + + inline std::vector parkTransformMatrix(ScalarT angle) + { + std::vector result(9); + ScalarT anpim = angle - (2.0/3.0)*M_PI; + ScalarT anpip = angle + (2.0/3.0)*M_PI; + result[0] = (2.0/3.0)*cos(angle); + result[1] = (2.0/3.0)*cos(anpim); + result[2] = (2.0/3.0)*cos(anpip); + result[3] = (2.0/3.0)*sin(angle); + result[4] = (2.0/3.0)*sin(anpim); + result[5] = (2.0/3.0)*sin(anpip); + result[6] = 1.0/3.0; + result[7] = 1.0/3.0; + result[8] = 1.0/3.0; + return result; + } + + inline std::vector parkTransformMatrixDerivative(ScalarT angle) + { + std::vector result(9); + ScalarT anpim = angle - (2.0/3.0)*M_PI; + ScalarT anpip = angle + (2.0/3.0)*M_PI; + result[0] = (2.0/3.0)*sin(angle); + result[1] = (2.0/3.0)*sin(anpim); + result[2] = (2.0/3.0)*sin(anpip); + result[3] = (2.0/3.0)*-cos(angle); + result[4] = (2.0/3.0)*-cos(anpim); + result[5] = (2.0/3.0)*-cos(anpip); + result[6] = 0; + result[7] = 0; + result[8] = 0; + return result; + } + + inline std::vector inverseParkTransformMatrix(ScalarT angle) + { + std::vector result(9); + ScalarT anpim = angle - (2.0/3.0)*M_PI; + ScalarT anpip = angle + (2.0/3.0)*M_PI; + result[0] = cos(angle); + result[1] = sin(angle); + result[2] = 1.0; + result[3] = cos(anpim); + result[4] = sin(anpim); + result[5] = 1.0; + result[6] = cos(anpip); + result[7] = sin(anpip); + result[8] = 1.0; + return result; + } + + inline std::vector inverseParkTransformMatrixDerivative(ScalarT angle) + { + std::vector result(9); + ScalarT anpim = angle - (2.0/3.0)*M_PI; + ScalarT anpip = angle + (2.0/3.0)*M_PI; + result[0] = sin(angle); + result[1] = -cos(angle); + result[2] = 0.0; + result[3] = sin(anpim); + result[4] = -cos(anpim); + result[5] = 0.0; + result[6] = sin(anpip); + result[7] = -cos(anpip); + result[8] = 0.0; + return result; + } + + 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/DiscreteGenerator/CMakeLists.txt b/ComponentLib/PowerElectronicsComponents/DiscreteGenerator/CMakeLists.txt new file mode 100644 index 000000000..c0ed2c7a9 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/DiscreteGenerator/CMakeLists.txt @@ -0,0 +1,8 @@ + + + +gridkit_add_library(powerelec_disgen + SOURCES + DiscreteGenerator.cpp + OUTPUT_NAME + gridkit_powerelec_disgen) \ No newline at end of file diff --git a/ComponentLib/PowerElectronicsComponents/DiscreteGenerator/DiscreteGenerator.cpp b/ComponentLib/PowerElectronicsComponents/DiscreteGenerator/DiscreteGenerator.cpp new file mode 100644 index 000000000..038a24ea2 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/DiscreteGenerator/DiscreteGenerator.cpp @@ -0,0 +1,348 @@ + + + +#include +#include +#include +#include "DiscreteGenerator.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 +DiscreteGenerator::DiscreteGenerator(IdxT id, DiscreteGeneratorParameters 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] + this->size_ = 16; + this->n_intern = 13; + this->n_extern = 3; + this->extern_indices = {0,1,2}; + this->idc_ = id; +} + +template +DiscreteGenerator::~DiscreteGenerator() +{ +} + +/*! + * @brief allocate method computes sparsity pattern of the Jacobian. + */ +template +int DiscreteGenerator::allocate() +{ + this->y_.resize(this->size_); + this->yp_.resize(this->size_); + this->f_.resize(this->size_); + + return 0; +} + +/** + * Initialization of the grid model + */ +template +int DiscreteGenerator::initialize() +{ + return 0; +} + +/* + * \brief Identify differential variables + */ +template +int DiscreteGenerator::tagDifferentiable() +{ + return 0; +} + +/** + * @brief Contributes to the bus residual. + * + * Must be connected to a PQ bus. + */ +template +int DiscreteGenerator::evaluateResidual() +{ + // ### Externals Componenets ### + + ScalarT omega = wb_ - mp_ * y_[4]; + //ref common ref motor angle + 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 ## + // Rotator 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 DiscreteGenerator 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 DiscreteGenerator::evaluateJacobian() +{ + this->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])}; + this->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])}; + this->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_}; + this->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_}; + this->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]}; + this->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]}; + this->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}; + this->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}; + this->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_}; + this->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_}; + this->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_}; + this->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_}; + this->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_}; + this->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_}; + this->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]}; + this->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_}; + this->J_.setValues(rtemp, ctemp, valtemp); + + + //Perform dF/dy + \alpha dF/dy' + + this->J_.AXPY(this->alpha_, Jacder); + + return 0; +} + +template +int DiscreteGenerator::evaluateIntegrand() +{ + return 0; +} + +template +int DiscreteGenerator::initializeAdjoint() +{ + return 0; +} + +template +int DiscreteGenerator::evaluateAdjointResidual() +{ + return 0; +} + +template +int DiscreteGenerator::evaluateAdjointIntegrand() +{ + return 0; +} + + + + +// Available template instantiations +template class DiscreteGenerator; +template class DiscreteGenerator; + + +} //namespace ModelLib + diff --git a/ComponentLib/PowerElectronicsComponents/DiscreteGenerator/DiscreteGenerator.hpp b/ComponentLib/PowerElectronicsComponents/DiscreteGenerator/DiscreteGenerator.hpp new file mode 100644 index 000000000..bca984460 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/DiscreteGenerator/DiscreteGenerator.hpp @@ -0,0 +1,99 @@ + + +#ifndef _CAP_HPP_ +#define _CAP_HPP_ + +#include +#include +#include + + +namespace ModelLib +{ + template class BaseBus; + + template + struct DiscreteGeneratorParameters + { + 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 passive DiscreteGenerator class. + * + */ + template + class DiscreteGenerator : 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_; + + + public: + DiscreteGenerator(IdxT id, DiscreteGeneratorParameters parm, bool reference_frame); + virtual ~DiscreteGenerator(); + + 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 000000000..95ef42436 --- /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 000000000..75782b435 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/InductionMotor/InductionMotor.cpp @@ -0,0 +1,135 @@ + + + +#include +#include +#include +#include "InductionMotor.hpp" + + +namespace ModelLib { + + + +/*! + * @brief Constructor for a constant InductionMotor model + * + * Calls default ModelEvaluatorImpl constructor. + */ + +template +InductionMotor::InductionMotor(IdxT id, ScalarT Lls, ScalarT Rs, ScalarT Llr, ScalarT Rr, ScalarT Lms, ScalarT J, ScalarT P) + : Lls_(Lls), + Rs_(Rs), + Llr_(Llr), + Rr_(Rr), + Lms_(Lms), + J_(J), + P_(P) +{ + this->size_ = 10; + this->n_intern = 5; + this->n_extern = 5; + this->extern_indices = {0,1,2,3,4}; + this->idc_ = id; +} + +template +InductionMotor::~InductionMotor() +{ +} + +/*! + * @brief allocate method computes sparsity pattern of the Jacobian. + */ +template +int InductionMotor::allocate() +{ + this->y_.resize(this->size_); + this->yp_.resize(this->size_); + this->f_.resize(this->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 bus residual. + * + * Must be connected to a PQ bus. + */ +template +int InductionMotor::evaluateResidual() +{ + //Note this leaves induction lumped into y. Perhaps would be better to seperate volatge and induction into seperate vectors + // for easier development + this->f_[0] = y_[5] + y_[7]; + this->f_[1] = (-1.0/2.0) * y_[5] - (sqrt(3.0)/2.0)*y_[6] + y_[7]; + this->f_[2] = (-1.0/2.0) * y_[5] + (sqrt(3.0)/2.0)*y_[6] + y_[7]; + this->f_[3] = J_ * yp_[3] - (3.0/4.0)*P_*Lms_ * (y_[5]*y_[9] - y_[6]*y_[8]); + this->f_[4] = yp_[4] - y_[3]; + this->f_[5] = (1.0/3.0)*(2.0* y_[0] - y_[1] - y_[2]) - Rs_*y_[5] - (Lls_ + Lms_) * yp_[5] - Lms_ * yp_[6]; + this->f_[6] = (1.0/sqrt(3.0))*(-y_[1] + y_[2]) - Rs_*y_[6] - (Lls_ + Lms_) * yp_[6] - Lms_ * yp_[5]; + this->f_[7] = (y_[0] + y_[1] + y_[2])/3.0 - Rs_*y_[7] - Lls_ * yp_[7]; + this->f_[8] = Rr_*y_[8] + (Llr_ + Lms_)*yp_[8] + Lms_ * yp_[5] - (P_/2)*y_[3]*((Llr_+Lms_)*y_[9] + Lms_*y_[6]); + this->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; +} + +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 000000000..2547754e2 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/InductionMotor/InductionMotor.hpp @@ -0,0 +1,69 @@ + + +#ifndef _IMOTOR_HPP_ +#define _IMOTOR_HPP_ + +#include +#include +#include + + +namespace ModelLib +{ + template class BaseBus; +} + + +namespace ModelLib +{ + /*! + * @brief Declaration of a passive 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::param_; + using CircuitComponent::idc_; + + public: + InductionMotor(IdxT id, ScalarT Lls, ScalarT Rs, ScalarT Llr, ScalarT Rr, ScalarT Lms, ScalarT J, 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 J_; + ScalarT P_; + }; +} + +#endif diff --git a/ComponentLib/PowerElectronicsComponents/Inductor/CMakeLists.txt b/ComponentLib/PowerElectronicsComponents/Inductor/CMakeLists.txt new file mode 100644 index 000000000..4c620030f --- /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 000000000..885153d41 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/Inductor/Inductor.cpp @@ -0,0 +1,138 @@ + + + +#include +#include +#include +#include "Inductor.hpp" + +namespace ModelLib { + +/*! + * @brief Constructor for a constant load model + * + * Calls default ModelEvaluatorImpl constructor. + */ + +template +Inductor::Inductor(IdxT id, ScalarT L) + : L_(L) +{ + this->size_ = 3; + this->n_intern = 1; + this->n_extern = 2; + this->extern_indices = {0,1}; + this->idc_ = id; +} + +template +Inductor::~Inductor() +{ +} + +/*! + * @brief allocate method computes sparsity pattern of the Jacobian. + */ +template +int Inductor::allocate() +{ + + this->y_.resize(this->size_); + this->yp_.resize(this->size_); + this->f_.resize(this->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 Contributes to the bus residual. + * + * Must be connected to a PQ bus. + */ +template +int Inductor::evaluateResidual() +{ + //input + this->f_[0] = -this->y_[2]; + //output + this->f_[1] = this->y_[2]; + //internal + this->f_[2] = -this->L_ * this->yp_[2] + this->y_[1] - this->y_[0] ; + return 0; +} + +template +int Inductor::evaluateJacobian() +{ + this->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}; + this->J_.setValues(rcord, ccord, vals); + + //Create dF/dy' + std::vector rcordder{2}; + std::vector ccordder{2}; + std::vector valsder{-this->L_}; + COO_Matrix Jacder = COO_Matrix(rcordder, ccordder, valsder,3,3); + + //Perform dF/dy + \alpha dF/dy' + this->J_.AXPY(this->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 000000000..2f2c1085b --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/Inductor/Inductor.hpp @@ -0,0 +1,67 @@ + + +#ifndef _IND_HPP_ +#define _IND_HPP_ + +#include +#include +#include + + +namespace ModelLib +{ + template class BaseBus; +} + + +namespace ModelLib +{ + /*! + * @brief Declaration of a passive 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_; + + + + 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 000000000..eaf95c043 --- /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 000000000..4caaeae66 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/LinearTransformer/LinearTransformer.cpp @@ -0,0 +1,124 @@ + + + +#include +#include +#include +#include "LinearTransformer.hpp" + +namespace ModelLib { + +/*! + * @brief Constructor for a constant LinearTransformer model + * + * Calls default ModelEvaluatorImpl constructor. + */ + +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) +{ + this->size_ = 4; + this->n_intern = 2; + this->n_extern = 2; + this->extern_indices = {0,1}; + this->idc_ = id; +} + +template +LinearTransformer::~LinearTransformer() +{ +} + +/*! + * @brief allocate method computes sparsity pattern of the Jacobian. + */ +template +int LinearTransformer::allocate() +{ + this->y_.resize(this->size_); + this->yp_.resize(this->size_); + this->f_.resize(this->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 Contributes to the bus residual. + * + * Must be connected to a PQ bus. + */ +template +int LinearTransformer::evaluateResidual() +{ + //Note this leaves induction lumped into y. Perhaps would be better to seperate volatge and induction into seperate vectors + // for easier development + this->f_[0] = this->y_[2]; + this->f_[1] = this->y_[3]; + this->f_[2] = this->y_[0] - this->R1_ * this->y_[2] - this->L1_ * this->yp_[2] - this->M_ * this->yp_[3]; + this->f_[2] = this->y_[1] - this->R2_ * this->y_[3] - this->M_ * this->yp_[2] - this->L2_ * this->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 000000000..cf56adfa2 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/LinearTransformer/LinearTransformer.hpp @@ -0,0 +1,67 @@ + + +#ifndef _LTRANS_HPP_ +#define _LTRANS_HPP_ + +#include +#include +#include + + +namespace ModelLib +{ + template class BaseBus; +} + + +namespace ModelLib +{ + /*! + * @brief Declaration of a passive 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::param_; + using CircuitComponent::idc_; + + 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 000000000..1bfb8dc3d --- /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 000000000..206108221 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/MicrogridBusDQ/MicrogridBusDQ.cpp @@ -0,0 +1,132 @@ + +#include +#include +#include +#include "MicrogridBusDQ.hpp" + +namespace ModelLib { + +/*! + * @brief Constructor for a constant MicrogridBusDQ model + * + * Calls default ModelEvaluatorImpl constructor. + * + * Each microgrid line has a virtual resistance RN + */ + +template +MicrogridBusDQ::MicrogridBusDQ(IdxT id, ScalarT RN) + : RN_(RN) +{ + // externals [vbus_d, vbus_q] + this->size_ = 2; + this->n_intern = 0; + this->n_extern = 2; + this->extern_indices = {0,1}; + this->idc_ = id; +} + +template +MicrogridBusDQ::~MicrogridBusDQ() +{ +} + +/*! + * @brief allocate method computes sparsity pattern of the Jacobian. + */ +template +int MicrogridBusDQ::allocate() +{ + this->y_.resize(this->size_); + this->yp_.resize(this->size_); + this->f_.resize(this->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" on both ends with parameters RN1 and RN2 + * + */ +template +int MicrogridBusDQ::evaluateResidual() +{ + //bus voltage + f_[0] = -y_[0] / RN_; + f_[1] = -y_[1] / RN_; + + return 0; +} + +/** + * @brief Generate Jacobian for Transmission Line + * + * @tparam ScalarT + * @tparam IdxT + * @return int + */ +template +int MicrogridBusDQ::evaluateJacobian() +{ + this->J_.zeroMatrix(); + + //Create dF/dy + std::vector rtemp{0,1}; + std::vector ctemp{0,1}; + std::vector vals{-1.0 / RN_,-1.0 / RN_}; + this->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 000000000..44154da38 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/MicrogridBusDQ/MicrogridBusDQ.hpp @@ -0,0 +1,66 @@ + + +#ifndef _VIRBUSDQ_HPP_ +#define _VIRBUSDQ_HPP_ + +#include +#include +#include + + +namespace ModelLib +{ + template class BaseBus; +} + + +namespace ModelLib +{ + /*! + * @brief Declaration of a passive 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_; + + + 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 000000000..d59cebb9c --- /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 000000000..e040ceb70 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/MicrogridLine/MicrogridLine.cpp @@ -0,0 +1,172 @@ + +#include +#include +#include +#include "MicrogridLine.hpp" + +namespace ModelLib { + +/*! + * @brief Constructor for a constant MicrogridLine model + * + * Calls default ModelEvaluatorImpl constructor. + * + * Each microgrid line has a virtual resistance RN + */ + +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] + this->size_ = 7; + this->n_intern = 2; + this->n_extern = 5; + this->extern_indices = {0,1,2,3,4}; + this->idc_ = id; +} + +template +MicrogridLine::~MicrogridLine() +{ +} + +/*! + * @brief allocate method computes sparsity pattern of the Jacobian. + */ +template +int MicrogridLine::allocate() +{ + this->y_.resize(this->size_); + this->yp_.resize(this->size_); + this->f_.resize(this->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 + * This model has "Virtual resistors" on both ends with parameters RN1 and RN2 + * + */ +template +int MicrogridLine::evaluateResidual() +{ + //ref motor + this->f_[0] = 0.0; + + //input + this->f_[1] = -y_[5] ; + this->f_[2] = -y_[6] ; + + //output + this->f_[3] = y_[5] ; + this->f_[4] = y_[6] ; + + //Internal variables + this->f_[5] = -yp_[5] - (R_ / L_) * y_[5] + y_[0]*y_[6] + (y_[1] - y_[3])/L_; + this->f_[6] = -yp_[6] - (R_ / L_) * y_[6] - y_[0]*y_[5] + (y_[2] - y_[4])/L_; + + + return 0; +} + +/** + * @brief Generate Jacobian for Transmission Line + * + * @tparam ScalarT + * @tparam IdxT + * @return int + */ +template +int MicrogridLine::evaluateJacobian() +{ + this->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}; + this->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]}; + this->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_)}; + this->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' + this->J_.AXPY(this->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 000000000..11ef78617 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/MicrogridLine/MicrogridLine.hpp @@ -0,0 +1,67 @@ + + +#ifndef _TRANLINE_HPP_ +#define _TRANLINE_HPP_ + +#include +#include +#include + + +namespace ModelLib +{ + template class BaseBus; +} + + +namespace ModelLib +{ + /*! + * @brief Declaration of a passive 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_; + + + 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 000000000..313d47688 --- /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 000000000..c39f53388 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/MicrogridLoad/MicrogridLoad.cpp @@ -0,0 +1,169 @@ + +#include +#include +#include +#include "MicrogridLoad.hpp" + +namespace ModelLib { + +/*! + * @brief Constructor for a constant MicrogridLoad 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 + */ + +template +MicrogridLoad::MicrogridLoad(IdxT id, ScalarT R,ScalarT L) + : R_(R), L_(L) +{ + // internals [id, iq] + // externals [\omegaref, vbd_out, vbq_out] + this->size_ = 5; + this->n_intern = 2; + this->n_extern = 3; + this->extern_indices = {0,1,2}; + this->idc_ = id; + +} + +template +MicrogridLoad::~MicrogridLoad() +{ +} + +/*! + * @brief allocate method computes sparsity pattern of the Jacobian. + */ +template +int MicrogridLoad::allocate() +{ + this->y_.resize(this->size_); + this->yp_.resize(this->size_); + this->f_.resize(this->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 + this->f_[0] = 0.0; + + //only input for loads + + //input + this->f_[1] = -y_[3] ; + this->f_[2] = -y_[4] ; + + //Internal variables + this->f_[3] = -yp_[3] - (R_ / L_) * y_[3] + y_[0]*y_[4] + y_[1] / L_; + this->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() +{ + this->J_.zeroMatrix(); + + //Create dF/dy + std::vector rtemp{1,2}; + std::vector ctemp{3,4}; + std::vector vals{-1.0,-1.0}; + this->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]}; + this->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_)}; + this->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' + this->J_.AXPY(this->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 000000000..22694f614 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/MicrogridLoad/MicrogridLoad.hpp @@ -0,0 +1,67 @@ + + +#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_; + + + 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 000000000..9386bda83 --- /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 000000000..45d9a34f3 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/Resistor/Resistor.cpp @@ -0,0 +1,127 @@ + + + +#include +#include +#include +#include "Resistor.hpp" + +namespace ModelLib { + +/*! + * @brief Constructor for a constant resistor model + * + * Calls default ModelEvaluatorImpl constructor. + */ + +template +Resistor::Resistor(IdxT id, ScalarT R) + : R_(R) +{ + this->size_ = 2; + this->n_intern = 0; + this->n_extern = 2; + this->extern_indices = {0,1}; + this->idc_ = id; +} + +template +Resistor::~Resistor() +{ +} + +/*! + * @brief allocate method computes sparsity pattern of the Jacobian. + */ +template +int Resistor::allocate() +{ + this->y_.resize(this->size_); + this->yp_.resize(this->size_); + this->f_.resize(this->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 Contributes to the bus residual. + * + * Must be connected to a PQ bus. + */ +template +int Resistor::evaluateResidual() +{ + //input + this->f_[0] = (this->y_[0] - this->y_[1])/this->R_ ; + //ouput + this->f_[1] = (this->y_[1] - this->y_[0])/this->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 / this->R_, -1.0 / this->R_, -1.0 / this->R_, 1.0 / this->R_}; + this->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 000000000..984d8304b --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/Resistor/Resistor.hpp @@ -0,0 +1,66 @@ + + +#ifndef _RES_HPP_ +#define _RES_HPP_ + +#include +#include +#include + + +namespace ModelLib +{ + template class BaseBus; +} + + +namespace ModelLib +{ + /*! + * @brief Declaration of a passive 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_; + + + 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 000000000..bfaeace4b --- /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 000000000..05b1e59ac --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/SynchronousMachine/SynchronousMachine.cpp @@ -0,0 +1,151 @@ + + + +#include +#include +#include +#include "SynchronousMachine.hpp" + + +namespace ModelLib { + + + +/*! + * @brief Constructor for a constant SynchronousMachine model + * + * Calls default ModelEvaluatorImpl constructor. + */ + +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 J, ScalarT P, ScalarT mub) + : Lls_(Lls), + Llkq_(Llkq), + Llfd_(Llfd), + Llkd_(Llkd), + Lmq_(Lmq), + Lmd_(Lmd), + Rs_(Rs), + Rkq_(Rkq), + Rfd_(Rfd), + Rkd_(Rkd), + J_(J), + P_(P), + mub_(mub) +{ + this->size_ = 13; + this->n_intern = 6; + this->n_extern = 7; + this->extern_indices = {0,1,2,3,4}; + this->idc_ = id; +} + +template +SynchronousMachine::~SynchronousMachine() +{ +} + +/*! + * @brief allocate method computes sparsity pattern of the Jacobian. + */ +template +int SynchronousMachine::allocate() +{ + this->y_.resize(this->size_); + this->yp_.resize(this->size_); + this->f_.resize(this->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 Contributes to the bus residual. + * + * Must be connected to a PQ bus. + */ +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); + + this->f_[0] = y_[6]*cos1 + y_[7]*sin1 + y_[8]; + this->f_[1] = y_[6]*cos23m + y_[7]*sin23m + y_[8]; + this->f_[2] = y_[6]*cos23p + y_[7]*sin23p + y_[8]; + this->f_[3] = J_ * 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])); + this->f_[4] = yp_[5] - y_[4]; + this->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]); + this->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]); + this->f_[7] = (-1.0/3.0)*(y_[0] + y_[1] + y_[2]) + Rs_*y_[8] + Lls_*yp_[8]; + this->f_[8] = rkq1*y_[9] + (llkq1 + Lmq_)*yp_[9] + Lmq_*yp_[6] + Lmq_*yp_[10]; + this->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 000000000..a18eca9c9 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/SynchronousMachine/SynchronousMachine.hpp @@ -0,0 +1,76 @@ + + +#ifndef _SYNMACH_HPP_ +#define _SYNMACH_HPP_ + +#include +#include +#include +#include + + +namespace ModelLib +{ + template class BaseBus; +} + + +namespace ModelLib +{ + /*! + * @brief Declaration of a passive 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::param_; + using CircuitComponent::idc_; + + 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 J, 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 J_; + ScalarT P_; + ScalarT mub_; + }; +} + +#endif diff --git a/ComponentLib/PowerElectronicsComponents/TransmissionLine/CMakeLists.txt b/ComponentLib/PowerElectronicsComponents/TransmissionLine/CMakeLists.txt new file mode 100644 index 000000000..6e2cea4f4 --- /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 000000000..62965c68e --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/TransmissionLine/TransmissionLine.cpp @@ -0,0 +1,201 @@ + +#include +#include +#include +#include "TransmissionLine.hpp" + +namespace ModelLib { + +/*! + * @brief Constructor for a constant 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 + */ + +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] + this->size_ = 12; + this->n_intern = 4; + this->n_extern = 8; + this->extern_indices = {0,1,2,3,4,5,6,7}; + this->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() +{ + this->y_.resize(this->size_); + this->yp_.resize(this->size_); + this->f_.resize(this->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 + this->f_[0] = y_[8] ; + this->f_[1] = y_[9] ; + + this->f_[2] = y_[10] ; + this->f_[3] = y_[11] ; + //ouput + this->f_[4] = -y_[8] ; + this->f_[5] = -y_[9] ; + + this->f_[6] = -y_[10] ; + this->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 + this->f_[8] = YReMat_ * (V1re - V2re) - (YImMatDi_ * V1im + YImMatOff_ * V2im) - y_[8] ; + this->f_[9] = YReMat_ * (V1im - V2im) + (YImMatDi_ * V1re + YImMatOff_ * V2re) - y_[9] ; + + //row2 + this->f_[10] = YReMat_ * (V2re - V1re) - (YImMatOff_ * V1im + YImMatDi_ * V2im) - y_[10]; + this->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}; + this->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_}; + this->J_.setValues(rtemp, ctemp, vals); + + + std::fill(rcord.begin(), rcord.end(), 9); + vals = {YImMatDi_ ,YReMat_, YImMatOff_, -YReMat_,-YImMatDi_ ,-YReMat_, -YImMatOff_, YReMat_}; + this->J_.setValues(rtemp, ctemp, vals); + + + std::fill(rcord.begin(), rcord.end(), 10); + vals = {-YReMat_, -YImMatDi_ ,YReMat_, -YImMatOff_,YReMat_, YImMatDi_ ,-YReMat_, YImMatOff_}; + this->J_.setValues(rtemp, ctemp, vals); + + + std::fill(rcord.begin(), rcord.end(), 11); + vals = {YImMatDi_ ,-YReMat_, YImMatOff_, YReMat_,-YImMatDi_ ,YReMat_, -YImMatOff_, -YReMat_}; + this->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 000000000..0f729d449 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/TransmissionLine/TransmissionLine.hpp @@ -0,0 +1,71 @@ + + +#ifndef _TRANLINE_HPP_ +#define _TRANLINE_HPP_ + +#include +#include +#include + + +namespace ModelLib +{ + template class BaseBus; +} + + +namespace ModelLib +{ + /*! + * @brief Declaration of a passive TransmissionLine class. + * + */ + 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_; + + + 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 000000000..7196f4d43 --- /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 000000000..e58e50556 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/VoltageSource/VoltageSource.cpp @@ -0,0 +1,129 @@ + + + +#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) +{ + this->size_ = 3; + this->n_intern = 1; + this->n_extern = 2; + this->extern_indices = {0,1}; + this->idc_ = id; +} + +template +VoltageSource::~VoltageSource() +{ +} + +/*! + * @brief allocate method computes sparsity pattern of the Jacobian. + */ +template +int VoltageSource::allocate() +{ + this->y_.resize(this->size_); + this->yp_.resize(this->size_); + this->f_.resize(this->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 Contributes to the bus residual. + * + * Must be connected to a PQ bus. + */ +template +int VoltageSource::evaluateResidual() +{ + //Note this leaves induction lumped into y. Perhaps would be better to seperate volatge and induction into seperate vectors + // for easier development + //input + this->f_[0] = -this->y_[2]; + //ouput + this->f_[1] = this->y_[2]; + //internal + this->f_[2] = this->y_[1] - this->y_[0] - this->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}; + this->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 000000000..3ac1c8699 --- /dev/null +++ b/ComponentLib/PowerElectronicsComponents/VoltageSource/VoltageSource.hpp @@ -0,0 +1,64 @@ + + +#ifndef _VOSO_HPP_ +#define _VOSO_HPP_ + +#include +#include +#include + + +namespace ModelLib +{ + template class BaseBus; +} + + +namespace ModelLib +{ + /*! + * @brief Declaration of a passive 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_; + + 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 feec6623a..3df418ba2 100644 --- a/Examples/CMakeLists.txt +++ b/Examples/CMakeLists.txt @@ -56,6 +56,11 @@ # add_subdirectory(MatPowerTesting) +add_subdirectory(RLCircuit) +add_subdirectory(Microgrid) +add_subdirectory(Scale_Microgrid) +add_subdirectory(SparseTest) +add_subdirectory(DiscreteGeneratorTest) if(TARGET SUNDIALS::kinsol) add_subdirectory(Grid3Bus) diff --git a/Examples/DiscreteGeneratorTest/CMakeLists.txt b/Examples/DiscreteGeneratorTest/CMakeLists.txt new file mode 100644 index 000000000..9238fb346 --- /dev/null +++ b/Examples/DiscreteGeneratorTest/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 DiscreteGeneratorTest COMMAND $) +install(TARGETS dgtest RUNTIME DESTINATION bin) diff --git a/Examples/DiscreteGeneratorTest/DGTest.cpp b/Examples/DiscreteGeneratorTest/DGTest.cpp new file mode 100644 index 000000000..d3dd2c68b --- /dev/null +++ b/Examples/DiscreteGeneratorTest/DGTest.cpp @@ -0,0 +1,78 @@ + + +#include +#include +#include +#include +#include +#include + +#include + + +int main(int argc, char const *argv[]) +{ + + ModelLib::DiscreteGeneratorParameters 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::DiscreteGenerator *dg = new ModelLib::DiscreteGenerator(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/Grid3Bus/3bus.mat b/Examples/Grid3Bus/3bus.mat new file mode 100644 index 000000000..4663eb33c --- /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 02985ff1d..654b7eb7e 100644 --- a/Examples/Grid3Bus/Grid3BusSys.cpp +++ b/Examples/Grid3Bus/Grid3BusSys.cpp @@ -89,6 +89,11 @@ static const std::string BUS3_DATA_STRING = R"( function mpc = case5 % Created by Reid Gomillion +//Note: This was traced from the subsequent calls +static const std::string BUS3_DATA_STRING = R"( +function mpc = case5 +% Created by Reid Gomillion + % MATPOWER %% MATPOWER Case Format : Version 2 @@ -144,6 +149,7 @@ 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 @@ -181,12 +187,12 @@ int monolithic_case() kinsol->printFinalStats(); int retval1 = 0; - retval1 += !isEqual(th2, theta2_ref, 1e-4); - retval1 += !isEqual(V2, V2_ref , 1e-4); - retval1 += !isEqual(th3, theta3_ref, 1e-4); + retval1 += !isEqual(th2, -4.87979, 1e-4); + retval1 += !isEqual(V2, 1.08281, 1e-4); + retval1 += !isEqual(th3, 1.46241, 1e-4); if(retval1 == 0) - std::cout << "\nSuccess!\n\n\n"; + std::cout << "\nSucess!\n\n\n"; else std::cout << "\nFailed!\n\n\n"; @@ -235,20 +241,20 @@ int parser_case() std::cout << "Solution:\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"; + 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"; // Print solver performance statistics kinsol->printFinalStats(); int retval2 = 0; - retval2 += !isEqual(th2, theta2_ref, 1e-4); - retval2 += !isEqual(V2, V2_ref , 1e-4); - retval2 += !isEqual(th3, theta3_ref, 1e-4); + retval2 += !isEqual(th2, -4.87979, 1e-4); + retval2 += !isEqual(V2, 1.08281, 1e-4); + retval2 += !isEqual(th3, 1.46241, 1e-4); if(retval2 == 0) - std::cout << "\nSuccess!\n\n\n"; + std::cout << "\nSucess!\n\n\n"; else std::cout << "\nFailed!\n\n\n"; @@ -273,19 +279,19 @@ int hardwired_case() SystemSteadyStateModel* sysmodel = new SystemSteadyStateModel(); // Next create and add buses ... - // Create a slack bus, fix V=1, theta=0, bus ID = 1 + // Create a slack bus, fix V=1, theta=0, bus ID = 1" << std::endl; BusData bd1; bd1.bus_i = 1; bd1.type = 3; bd1.Vm = 1.0; bd1.Va = 0.0; auto* bus1 = BusFactory::create(bd1); sysmodel->addBus(bus1); - //Create a PQ bus, initialize V=1, theta=0, bus ID = 2 + //Create a PQ bus, initialize V=1, theta=0, bus ID = 2" << std::endl; BusData bd2; bd2.bus_i = 2; bd2.type = 1; bd2.Vm = 1.0; bd2.Va = 0.0; auto* bus2 = BusFactory::create(bd2); sysmodel->addBus(bus2); - // Create a PV bus, fix V=1.1, initialize theta=0, and set power injection Pg=2 + // Create a PV bus, fix V=1.1, initialize theta=0, and set power injection Pg=2" << std::endl; BusData bd3; bd3.bus_i = 3; bd3.type = 2; bd3.Vm = 1.1; bd3.Va = 0.0; auto* bus3 = BusFactory::create(bd3); @@ -366,12 +372,12 @@ int hardwired_case() kinsol->printFinalStats(); int retval2 = 0; - retval2 += !isEqual(th2, theta2_ref, 1e-4); - retval2 += !isEqual(V2, V2_ref , 1e-4); - retval2 += !isEqual(th3, theta3_ref, 1e-4); + retval2 += !isEqual(th2, -4.87979, 1e-4); + retval2 += !isEqual(V2, 1.08281, 1e-4); + retval2 += !isEqual(th3, 1.46241, 1e-4); if(retval2 == 0) - std::cout << "\nSuccess!\n\n\n"; + std::cout << "\nSucess!\n\n\n"; else std::cout << "\nFailed!\n\n\n"; @@ -386,6 +392,7 @@ int hardwired_case() int main() { //return the results of each case + //swapping orders of test causes memory error, specifically hardware <-> parser int resolve = 0; std::cout << std::string(32,'-') << std::endl; resolve += monolithic_case(); diff --git a/Examples/Microgrid/CMakeLists.txt b/Examples/Microgrid/CMakeLists.txt new file mode 100644 index 000000000..0958144c4 --- /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 000000000..275974b19 --- /dev/null +++ b/Examples/Microgrid/Microgrid.cpp @@ -0,0 +1,444 @@ + + +#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; + size_t max_step_amount = 3000; + bool usejac = true; + + //TODO:setup as named parameters + //Create circuit 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::DiscreteGeneratorParameters 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::DiscreteGeneratorParameters 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::DiscreteGenerator *dg1 = new ModelLib::DiscreteGenerator(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::DiscreteGenerator *dg2 = new ModelLib::DiscreteGenerator(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::DiscreteGenerator *dg3 = new ModelLib::DiscreteGenerator(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::DiscreteGenerator *dg4 = new ModelLib::DiscreteGenerator(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("%u : %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-14 + 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("%u : %e ,\n", i, abs(true_vec[i] - yfinial[i]) / abs(true_vec[i])); + } + + return 0; +} diff --git a/Examples/RLCircuit/CMakeLists.txt b/Examples/RLCircuit/CMakeLists.txt new file mode 100644 index 000000000..476eb5e95 --- /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 000000000..55fe476b9 --- /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/Scale_Microgrid/CMakeLists.txt b/Examples/Scale_Microgrid/CMakeLists.txt new file mode 100644 index 000000000..8ffc91d81 --- /dev/null +++ b/Examples/Scale_Microgrid/CMakeLists.txt @@ -0,0 +1,13 @@ + + + + +add_executable(scalemicrogrid Scale_Microgrid.cpp) +target_link_libraries(scalemicrogrid GRIDKIT::powerelec_disgen + GRIDKIT::powerelec_mircoline + GRIDKIT::powerelec_microload + GRIDKIT::solvers_dyn + GRIDKIT::powerelec_mircobusdq) + +add_test(NAME ScaleMicrogrid COMMAND $) +install(TARGETS scalemicrogrid RUNTIME DESTINATION bin) diff --git a/Examples/Scale_Microgrid/Scale_Microgrid.cpp b/Examples/Scale_Microgrid/Scale_Microgrid.cpp new file mode 100644 index 000000000..9f1baef64 --- /dev/null +++ b/Examples/Scale_Microgrid/Scale_Microgrid.cpp @@ -0,0 +1,317 @@ + + +#include +#include +#include +#include +#include + + +#include +#include +#include +#include + +#include +#include +#include + +//macro to get macro +#define STRING(x) #x +#define XSTRING(x) STRING(x) + +int main(int argc, char const *argv[]) +{ + double abstol = 1.0e-8; + double reltol = 1.0e-8; + size_t max_step_amount = 3000; + bool usejac = true; + + //TODO:setup as named parameters + //Create circuit model + ModelLib::PowerElectronicsModel* sysmodel = new ModelLib::PowerElectronicsModel(reltol, abstol, usejac, max_step_amount); + + + //The amount of DG line load cobinations to generate for scale + size_t Nsize = 8; + + //Set to false if the parameters choosen are not used to generate data files + bool orgparams = true; + + //Modeled after the problem in the paper + // Every Bus has the same virtual resistance. This is due to the numerical stability as mentioned in the paper. + double RN = 1.0e4; + + //DG Params Vector + //All DGs have the same set of parameters except for the first two. + ModelLib::DiscreteGeneratorParameters DG_parms1; + DG_parms1.wb = 2.0*M_PI*50.0; + DG_parms1.wc = 31.41; + DG_parms1.mp = 9.4e-5; + DG_parms1.Vn = 380.0; + DG_parms1.nq = 1.3e-3; + DG_parms1.F = 0.75; + DG_parms1.Kiv = 420.0; + DG_parms1.Kpv = 0.1; + DG_parms1.Kic = 2.0e4; + DG_parms1.Kpc = 15.0; + DG_parms1.Cf = 5.0e-5; + DG_parms1.rLf = 0.1; + DG_parms1.Lf = 1.35e-3; + DG_parms1.rLc = 0.03; + DG_parms1.Lc = 0.35e-3; + + ModelLib::DiscreteGeneratorParameters DG_parms2; + DG_parms2.wb = 2.0*M_PI*50.0; + DG_parms2.wc = 31.41; + DG_parms2.mp = 12.5e-5; + DG_parms2.Vn = 380.0; + DG_parms2.nq = 1.5e-3; + DG_parms2.F = 0.75; + DG_parms2.Kiv = 390.0; + DG_parms2.Kpv = 0.05; + DG_parms2.Kic = 16.0e3; + DG_parms2.Kpc = 10.5; + DG_parms2.Cf = 50.0e-6; + DG_parms2.rLf = 0.1; + DG_parms2.Lf = 1.35e-3; + DG_parms2.rLc = 0.03; + DG_parms2.Lc = 0.35e-3; + + std::vector> DGParams_list(2*Nsize, DG_parms2); + + DGParams_list[0] = DG_parms1; + DGParams_list[1] = DG_parms1; + + //line vector params + //Every odd line has the same parameters and every even line has the same parameters + double rline1 = 0.23; + double Lline1 = 0.1 / (2.0 * M_PI * 50.0); + double rline2 = 0.35; + double Lline2 = 0.58 / (2.0 * M_PI * 50.0); + std::vector rline_list(2*Nsize-1, 0.0); + std::vector Lline_list(2*Nsize-1, 0.0); + for (size_t i = 0; i < rline_list.size(); i++) + { + rline_list[i] = (i % 2) ? rline2 : rline1; + Lline_list[i] = (i % 2) ? Lline2 : Lline1; + } + + + //load parms + //Only the first load has the same paramaters. + double rload1 = 3.0; + double Lload1 = 2.0 / (2.0 * M_PI * 50.0); + double rload2 = 2.0; + double Lload2 = 1.0 / (2.0 * M_PI * 50.0); + + std::vector rload_list(Nsize, rload2); + std::vector Lload_list(Nsize, Lload2); + rload_list[0] = rload1; + Lload_list[0] = Lload1; + + // 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); + + std::vector vdqbus_index(2*Nsize,0); + vdqbus_index[0] = vec_size_internals + 1; + for (size_t i = 1; i < vdqbus_index.size(); i++) + { + vdqbus_index[i] = vdqbus_index[i-1] + 2; + } + + //Total size of the vector setup + size_t vec_size_total = vec_size_internals + vec_size_externals; + + + //Create the reference DG + ModelLib::DiscreteGenerator *dg_ref = new ModelLib::DiscreteGenerator(0, DGParams_list[0], true); + //ref motor + dg_ref->setExternalConnectionNodes(0,vec_size_internals); + //outputs + dg_ref->setExternalConnectionNodes(1,vdqbus_index[0]); + dg_ref->setExternalConnectionNodes(2,vdqbus_index[0] + 1); + //"grounding" of the difference + dg_ref->setExternalConnectionNodes(3,-1); + //internal connections + for (size_t i = 0; i < 12; i++) + { + + dg_ref->setExternalConnectionNodes(4 + i,i); + } + sysmodel->addComponent(dg_ref); + + //Keep track of models and index location + size_t indexv = 12; + size_t model_id = 1; + //Add all other DGs + for (size_t i = 1; i < 2*Nsize; i++) + { + + //current DG to add + ModelLib::DiscreteGenerator *dg = new ModelLib::DiscreteGenerator(model_id++, DGParams_list[i], false); + //ref motor + dg->setExternalConnectionNodes(0,vec_size_internals); + //outputs + dg->setExternalConnectionNodes(1,vdqbus_index[i]); + dg->setExternalConnectionNodes(2,vdqbus_index[i] + 1); + //internal connections + for (size_t j = 0; j < 13; j++) + { + + dg->setExternalConnectionNodes(3 + j,indexv + j); + } + indexv += 13; + sysmodel->addComponent(dg); + } + + // Load all the Line compoenents + for (size_t i = 0; i < 2*Nsize - 1; i++) + { + //line + ModelLib::MicrogridLine *line_model = new ModelLib::MicrogridLine(model_id++, rline_list[i], Lline_list[i]); + //ref motor + line_model->setExternalConnectionNodes(0,vec_size_internals); + //input connections + line_model->setExternalConnectionNodes(1,vdqbus_index[i]); + line_model->setExternalConnectionNodes(2,vdqbus_index[i] + 1); + //output connections + line_model->setExternalConnectionNodes(3,vdqbus_index[i+1]); + line_model->setExternalConnectionNodes(4,vdqbus_index[i+1] + 1); + //internal connections + for (size_t j = 0; j < 2; j++) + { + line_model->setExternalConnectionNodes(5 + j,indexv + j); + } + indexv += 2; + sysmodel->addComponent(line_model); + } + + // Load all the Load components + for (size_t i = 0; i < Nsize; i++) + { + ModelLib::MicrogridLoad *load_model = new ModelLib::MicrogridLoad(model_id++, rload_list[i], Lload_list[i]); + //ref motor + load_model->setExternalConnectionNodes(0,vec_size_internals); + //input connections + load_model->setExternalConnectionNodes(1,vdqbus_index[2*i]); + load_model->setExternalConnectionNodes(2,vdqbus_index[2*i] + 1); + //internal connections + for (size_t j = 0; j < 2; j++) + { + + load_model->setExternalConnectionNodes(3 + j,indexv + j); + } + indexv += 2; + sysmodel->addComponent(load_model); + } + + //Add all the microgrid Virtual DQ Buses + for (size_t i = 0; i < 2*Nsize; i++) + { + ModelLib::MicrogridBusDQ *virDQbus_model = new ModelLib::MicrogridBusDQ( + model_id++, RN); + + virDQbus_model->setExternalConnectionNodes(0, vdqbus_index[i]); + virDQbus_model->setExternalConnectionNodes(1, vdqbus_index[i] + 1); + sysmodel->addComponent(virDQbus_model); + } + + //allocate all the intial conditions + 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. Every state is to specified to the zero intially + 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 + for (size_t i = 0; i < 2*Nsize; i++) + { + sysmodel->yp()[13*i - 1 + 3] = DGParams_list[i].Vn; + sysmodel->yp()[13*i - 1 + 5] = DGParams_list[i].Kpv * DGParams_list[i].Vn; + sysmodel->yp()[13*i - 1 + 7] = (DGParams_list[i].Kpc * DGParams_list[i].Kpv * DGParams_list[i].Vn) / DGParams_list[i].Lf; + } + + //since the intial P_com = 0, the set the intial vector to the reference frame + sysmodel->y()[vec_size_internals] = DG_parms1.wb; + + sysmodel->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("%u : %.16e \n", i, fres[i]); + } + std::cout << "}\n"; + + sysmodel->updateTime(0.0, 1.0e-8); + sysmodel->evaluateJacobian(); + std::cout << "Intial Jacobian with alpha:\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++) + { + printf("%u: % 2.16e\n", i, yfinial[i]); + } + + //if proper MATLAB results are avalible and model is using origional parameters + std::vector compDataSizes = {2,4,8}; + //All data generated with abstol=1e-12 and reltol=1e-12 with ode23tb + if (orgparams && compDataSizes.end() != std::find(compDataSizes.begin(), compDataSizes.end(), Nsize)) + { + //get the cmake source dir from the cmake configuration + std::string base_path = XSTRING(SOURCE_ROOT); + + base_path += "/Examples/Scale_Microgrid/microgrid_N"; + base_path += std::to_string(Nsize); + base_path += ".txt"; + std::cout << base_path << std::endl; + + //load vector from data file + std::vector true_vec(vec_size_total, 0.0); + double val_cur = 0.0; + std::ifstream data_file(base_path); + size_t i = 0; + while (data_file >> val_cur && i < true_vec.size()) + { + true_vec[i++] = val_cur; + } + data_file.close(); + + //check relative error + std::cout << "Test the Relative Error\n"; + for (size_t i = 0; i < true_vec.size(); i++) + { + printf("%u: % 2.16e\n", i, abs(true_vec[i] - yfinial[i]) / abs(true_vec[i])); + } + } + + + return 0; +} diff --git a/Examples/Scale_Microgrid/microgrid_N2.txt b/Examples/Scale_Microgrid/microgrid_N2.txt new file mode 100644 index 000000000..d0edf2a51 --- /dev/null +++ b/Examples/Scale_Microgrid/microgrid_N2.txt @@ -0,0 +1,70 @@ +22975.4182636905 +12753.1173017451 +0.0376305671131306 +-0.0209815421114537 +0.0184828563235406 +-0.00015632917065729 +63.2193617412383 +-29.4226556659582 +363.420924977355 +-3.01001804456668e-06 +63.2193618667302 +-35.0920143829141 +-0.00755583999797958 +22975.8457840701 +8742.0172166684 +0.0371009878165673 +-0.0142112776774797 +0.0187407972932474 +-9.8913852531819e-05 +62.3296548917654 +-17.9650261562766 +368.635407345931 +3.79041261513032e-05 +62.3296560841534 +-23.7157298601546 +-0.0827401584095097 +17277.712521392 +16493.7578328327 +0.0311645357391787 +-0.0298601153771274 +0.0225001003469904 +-0.000264388302121386 +48.6166482764357 +-40.8862612232969 +355.25957014696 +-0.000167008045114112 +48.6166481095298 +-46.4283022775928 +-0.0844566209033113 +17277.2493364959 +9182.29479881977 +0.0302503981389185 +-0.0161722538457184 +0.023180586071601 +-0.000129590157483103 +47.1906392248137 +-19.3574255862892 +366.226354188735 +0.000121054462396047 +47.19063792117 +-25.0705499054599 +-18.8125403122072 +21.1471334522508 +43.2997340692497 +-3.03798323340816 +-44.8715356440453 +28.9585224749495 +81.99613295014 +-56.2385627562355 +132.749774614687 +-82.280664729262 +311.999576042192 +357.692287974629 +-5.85078929333349 +364.119335219897 +-8.84631033126304 +347.214535767435 +-32.7241067379802 +360.411028950125 +-34.9283280833745 diff --git a/Examples/Scale_Microgrid/microgrid_N4.txt b/Examples/Scale_Microgrid/microgrid_N4.txt new file mode 100644 index 000000000..29a0b71b8 --- /dev/null +++ b/Examples/Scale_Microgrid/microgrid_N4.txt @@ -0,0 +1,142 @@ +27828.3291148094 +10602.3611530864 +0.0452452394028709 +-0.0173410584955486 +0.018686818367312 +-0.000129595027000909 +76.0120809944985 +-23.2367034227098 +366.216582730654 +-0.000167782960864797 +76.0120710728117 +-28.9413214479885 +-0.0154989336357971 +27832.0772154743 +7838.80352597054 +0.0448129806183038 +-0.0127296381144373 +0.0188642014642206 +-9.04502393904097e-05 +75.2859236720964 +-15.4316980534341 +369.809174828127 +-5.17440472732025e-05 +75.285919083362 +-21.1922719962368 +-0.142596248886844 +20950.1519748442 +16334.3924580166 +0.0377787488551238 +-0.0295717781600636 +0.0225780601074885 +-0.000265567743682681 +58.9348270364006 +-40.4081583028089 +355.498460539903 +-8.66167274613643e-05 +58.9348242715342 +-45.9457664944992 +-0.164628885385433 +20956.2846481359 +12345.7855779785 +0.0371673110634281 +-0.0220171188824754 +0.0229486760225139 +-0.000191111142238558 +57.9810673844105 +-28.5265491599081 +361.481100853144 +8.24481547667206e-05 +57.9810719040216 +-34.1573355283559 +-0.245126576498858 +20978.9206932923 +17875.6825958711 +0.0380785081131898 +-0.0325621468970297 +0.0224354615928625 +-0.000295085318936756 +59.4025243779293 +-45.1101946009202 +353.186117339098 +5.52721744209165e-05 +59.4025238741703 +-50.6117204294148 +-0.253672269778405 +20990.6142554844 +11779.5685373245 +0.0371243965775962 +-0.0209716906203726 +0.0230016494389548 +-0.000180846578022512 +57.9140614794815 +-26.8817831662193 +362.331254349606 +0.000209247259473183 +57.9140782926947 +-32.5257339504705 +-0.288747285297145 +21004.5418902462 +16686.981801999 +0.0379131316203435 +-0.0302602142538338 +0.0225460339603269 +-0.000272422675391038 +59.1444740212228 +-41.4902940591195 +354.969993387113 +8.14114753778542e-05 +59.1444846638236 +-47.0195494892644 +-0.285876486857154 +21005.895228757 +8274.13913126834 +0.0366267978764156 +-0.0145584209340773 +0.0233276948336748 +-0.000117627382397687 +57.1378568500881 +-16.7923279871216 +367.58883605301 +0.000263351556352378 +57.1378679159653 +-22.518140603782 +-6.57009481597077 +28.1053117219546 +68.341829672125 +5.75010361666728 +-7.27398138223966 +42.3413009908014 +44.2901915856089 +-0.850144510998897 +-26.9214526328474 +37.7939474726438 +20.9421230166698 +-8.2155394685309 +-48.4341275019927 +37.706248875936 +82.5460878639016 +-57.0458911301263 +127.388509162578 +-90.4406093871607 +116.522583626035 +-102.149734009655 +112.650118587972 +-107.825883544914 +311.543402422187 +360.780248808368 +-7.42039816792328 +365.078323407317 +-13.2331839185884 +344.46378568747 +-54.5546252592288 +350.335287885457 +-63.572021962256 +334.342198258035 +-88.749585714325 +344.281369936645 +-94.7726681838823 +332.2248206198 +-103.942854393111 +347.103721199815 +-107.812479611979 diff --git a/Examples/Scale_Microgrid/microgrid_N8.txt b/Examples/Scale_Microgrid/microgrid_N8.txt new file mode 100644 index 000000000..3e1ba779b --- /dev/null +++ b/Examples/Scale_Microgrid/microgrid_N8.txt @@ -0,0 +1,286 @@ +29266.6517718661 +9842.93617289812 +0.0475360674334622 +-0.0160323945316821 +0.0187557247105805 +-0.000119821248231538 +79.8612117213119 +-21.0144010498971 +367.201497494684 +-0.00135963018489767 +79.8611351716573 +-26.7318789917162 +-0.018074876389402 +29296.9645459408 +7415.51099954903 +0.0471789367393934 +-0.0120180219897291 +0.0189117902285582 +-8.57914733703415e-05 +79.2614266664798 +-14.2194852461262 +370.357116139458 +-0.00065261697325357 +79.261388452934 +-19.9860446882469 +-0.161178046356866 +22178.7337755581 +15750.48869542 +0.0399575496599213 +-0.028431003000158 +0.0226540351215628 +-0.000255811047442045 +62.3338080750298 +-38.6061338491152 +356.372404460464 +-0.000888421043595241 +62.3337589017493 +-44.1546192236143 +-0.191206982234175 +22245.3587939821 +12476.0312983393 +0.0395296431152022 +-0.0222832972915751 +0.0229591349358969 +-0.000195357283522886 +61.6667253463805 +-28.9363062720176 +361.284131747149 +0.000606091049765824 +61.6667589580497 +-34.5611135316627 +-0.300475535038271 +22471.354929649 +17410.9400696732 +0.0407475331761201 +-0.0316496339533813 +0.0225047347006168 +-0.000287941144857336 +63.5662158978541 +-43.6654709736205 +353.882194971394 +-0.00050508442766251 +63.5661867600733 +-49.1745383652835 +-0.321406735230742 +22562.1063953332 +12532.040739496 +0.0400730650460856 +-0.0224103188019356 +0.0229590093469276 +-0.000197040279593967 +62.5145400797163 +-29.1337602252211 +361.200953962968 +0.00123813270414561 +62.5146115608051 +-34.7565572443792 +-0.398765881347133 +22822.9625175864 +17526.9123376124 +0.0413769691613727 +-0.0318846988884131 +0.0224996948233429 +-0.000290724959665147 +64.5480523645119 +-44.0326735095454 +353.709220392926 +-0.000255137918253834 +64.548036917491 +-49.5382706457328 +-0.409294621792759 +22917.4640781459 +11700.1980576741 +0.0405370737861658 +-0.0208717910251266 +0.0230419643488667 +-0.000182321342412399 +63.2383663569188 +-26.7111380754575 +362.449323338104 +0.00143605100750772 +63.2384503029344 +-32.3525635247183 +-0.460038982115987 +23175.0870418607 +17429.0372806116 +0.0419717095772374 +-0.0317010913432995 +0.0225145433123769 +-0.000289381352876061 +65.4757968347883 +-43.7415548084125 +353.856643943042 +-0.000116590766441278 +65.4757901489451 +-49.2486679790579 +-0.461848591689221 +23264.6732405303 +10919.0360806903 +0.0409947737541995 +-0.0194331296392293 +0.0231200668519233 +-0.000168575863901554 +63.9522672718067 +-24.4455560147541 +363.621937642842 +0.00154556746783549 +63.9523565934322 +-30.1044436917781 +-0.491529580804623 +23488.6137650364 +17313.8654584652 +0.0424952351974662 +-0.0314883547842804 +0.0225303364455447 +-0.000287714847080755 +66.2924442955506 +-43.4048501801219 +354.030311475628 +0.000169934086992562 +66.2924545883603 +-48.9139721676174 +-0.486548093474613 +23557.6136117665 +10241.0790360359 +0.0413774884672506 +-0.0181860120867723 +0.0231877453638413 +-0.000156639599075328 +64.5492156984859 +-22.4816084883194 +364.639179755733 +0.00142420279185394 +64.5492984953198 +-28.1556597465189 +-0.499726325044563 +23722.4168131373 +16958.7765974266 +0.0428487548646152 +-0.0308090857078185 +0.0225672437721518 +-0.000281355782039133 +66.8440182379463 +-42.3350991467661 +354.562942832682 +0.000573219240561573 +66.844044569212 +-47.8519918014414 +-0.490219448053285 +23769.8747006891 +9267.5206785042 +0.0415639300026818 +-0.0164109903099032 +0.0232814268213158 +-0.000139385907988285 +64.8399200066663 +-19.6874999818575 +366.100410278813 +0.00128332197988872 +64.8400012051758 +-25.3838027935624 +-0.491586583798841 +23854.0240861864 +15160.7060140903 +0.0427385845553072 +-0.0273659210853948 +0.0227360281885303 +-0.00024756822257102 +66.6722121992594 +-36.9193366680176 +357.260483774461 +0.000734666517031364 +66.6722602990682 +-42.4779079553274 +-0.482907936930927 +23870.5548740189 +6166.60742794455 +0.0412253604361869 +-0.0108406830677542 +0.0235710268609798 +-8.45378382385454e-05 +64.3117067426907 +-10.9219886629635 +370.751596332923 +0.00130490570318174 +64.3117765573399 +-16.6904284244068 +-2.9267653635869 +30.5878043390112 +75.9238910016987 +9.17392754248112 +4.3701560041079 +48.7765765968327 +58.3101300775815 +3.13323809951331 +-6.37573932293514 +46.1701408549949 +41.9244199648313 +-6.5432868117422 +-17.3691293511209 +41.8611610718482 +27.7383173553277 +-12.9709065074016 +-27.5322418670623 +38.9032728624808 +16.2733736528389 +-16.5280348065045 +-36.4800564569504 +37.0627808819117 +7.38198504755967 +-17.9899033675727 +-44.0731708800056 +36.9539474625741 +1.14781892918688 +-15.9507601441003 +-49.1751072411967 +44.6278268007059 +82.7517114852138 +-57.3188928805757 +125.959171148664 +-93.1822145893453 +110.816672146529 +-108.811406393296 +99.5107131550674 +-119.104827911048 +92.0427992280762 +-125.056967118294 +88.0811872722293 +-127.985464658329 +87.1649116220972 +-128.957730829954 +89.0189233593932 +-129.480383815729 +311.408200092424 +361.890500304156 +-7.90450151715305 +365.594288034004 +-14.6504114199963 +344.278948204657 +-61.517809887448 +348.104769580715 +-73.1715704285651 +329.472137794467 +-107.797643650827 +335.509251816504 +-117.787067889008 +317.044841017058 +-139.602945170765 +325.182412285621 +-147.51178633583 +307.989533967543 +-158.91998673915 +318.171184823086 +-165.141414351062 +302.951769665185 +-168.712626964975 +315.009688292198 +-173.624287308414 +302.070219389208 +-171.572667424371 +315.86470341008 +-175.706947618597 +306.288884765777 +-170.785573848562 +322.018589824395 +-176.179397251772 diff --git a/Examples/SparseTest/CMakeLists.txt b/Examples/SparseTest/CMakeLists.txt new file mode 100644 index 000000000..d0b6f5e05 --- /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 000000000..81e5ac7dd --- /dev/null +++ b/Examples/SparseTest/SparseTest.cpp @@ -0,0 +1,86 @@ + + +#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--; + } + + std::cout << failval << std::endl; + + return failval; +} diff --git a/ModelEvaluator.hpp b/ModelEvaluator.hpp index 5eda481de..ad56fbd48 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 a4ef43a00..0f8969f0b 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 000000000..04508607f --- /dev/null +++ b/PowerElectronicsModel.hpp @@ -0,0 +1,331 @@ + + +#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 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 not 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 + this->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; jsize(); ++j) + { + if(component->getNodeConnection(j) != -1) + { + 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; jsize(); ++j) + { + //@todo should do a different grounding check + if (component->getNodeConnection(j) != -1) + { + 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() + { + this->J_.zeroMatrix(); + this->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]) != -1 && component->getNodeConnection(c[i]) != -1) + { + rgr.push_back(component->getNodeConnection(r[i])); + cgr.push_back(component->getNodeConnection(c[i])); + vgr.push_back(v[i]); + } + } + + //AXPY to Global Jacobian + this->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); + } + this->time_ = t; + this->alpha_ = a; + } + + void addComponent(component_type* component) + { + this->components_.push_back(component); + } + +private: + std::vector components_; + bool usejac_; + +}; // class PowerElectronicsModel + +} // namespace ModelLib diff --git a/Solver/Dynamic/CMakeLists.txt b/Solver/Dynamic/CMakeLists.txt index 759829967..85c7ea920 100644 --- a/Solver/Dynamic/CMakeLists.txt +++ b/Solver/Dynamic/CMakeLists.txt @@ -66,6 +66,11 @@ gridkit_add_library(solvers_dyn LINK_LIBRARIES PUBLIC SUNDIALS::nvecserial PUBLIC SUNDIALS::idas + PUBLIC SUNDIALS::sunlinsolklu + PUBLIC suitesparseconfig + PUBLIC klu + PUBLIC amd + PUBLIC colamd OUTPUT_NAME gridkit_solvers_dyn) diff --git a/Solver/Dynamic/Ida.cpp b/Solver/Dynamic/Ida.cpp index c71a8fd07..2be2913a9 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" @@ -103,6 +106,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 +133,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 +157,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 +166,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"); + + 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"); - 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; } @@ -400,7 +422,7 @@ namespace Sundials checkOutput(retval, "IDASetUserDataB"); /// \todo Need to set max number of steps based on user input! - retval = IDASetMaxNumStepsB(solver_, backwardID_, 2000); + retval = IDASetMaxNumStepsB(solver_, backwardID_, 200); checkOutput(retval, "IDASetMaxNumSteps"); // Set up linear solver @@ -525,8 +547,49 @@ 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) + { - template + 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) { ModelLib::ModelEvaluator* model = static_cast*>(user_data); diff --git a/Solver/Dynamic/Ida.hpp b/Solver/Dynamic/Ida.hpp index 3fa3368c9..a4cea9561 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 000000000..b51d669ec --- /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 000000000..716748fa8 --- /dev/null +++ b/SparseMatrix/COO_Matrix.hpp @@ -0,0 +1,741 @@ + + +#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 &values); + +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 + * + * @todo fails, fix + * + * @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]); + aiter++; + + this->checkIncreaseSize(r[aiter], c[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; +} + +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"; + } +} + +/** + * @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 &values) +{ + + //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(values[ordervec[i]], values[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 34c7ff120..77bbdd20b 100644 --- a/SystemModel.hpp +++ b/SystemModel.hpp @@ -183,6 +183,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. *