diff --git a/CHANGELOG.md b/CHANGELOG.md index 0fd168618..9dfb4637c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -52,6 +52,7 @@ - Added component model developer checklist to a README file. - Added IEEEST Stabilizer Model - Added SEXS-PTI Exciter Model +- Added GENSAL Machine Model - Added 200 Bus Synthetic Illinois Case - Added node objects to `PowerElectronics` module & updated all examples to make use of them. diff --git a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp index 6feca8999..b44819ed5 100644 --- a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp +++ b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp @@ -11,4 +11,5 @@ #include #include #include +#include #include diff --git a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md index ec01498f2..527d9ef84 100644 --- a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md +++ b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md @@ -144,6 +144,7 @@ are specified: `Branch` | a basic algebraic pi model for a line or transformer | `bus1`, `bus2` | `R`, `X`, `G`, `B` | `ir1`, `ii1`, `im1`, `p1`, `q1`, `ir2`, `ii2`, `im2`, `p2`, `q2` `Load` | a basic static impedence load model | `bus` | `R`, `X` | `p`, `q` `Genrou` | 6th order machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Tdop`, `Tdopp`, `Tqopp`, `Tqop`, `Xd`, `Xdp`, `Xdpp`, `Xq`, `Xqp`, `Xqpp`, `Xl`, `S10`, `S12`, `mva_base` | `ir`, `ii`, `p`, `q`, `delta`, `omega`, `speed` + `Gensal` | 5th order salient-pole machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Tdop`, `Tdopp`, `Tqopp`, `Xd`, `Xdp`, `Xdpp`, `Xq`, `Xl`, `S10`, `S12`, `mva_base` | `ir`, `ii`, `p`, `q`, `delta`, `omega`, `speed` `GenClassical`| the classical machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Xdp`, `mva_base` | `ir`, `ii`, `p`, `q`, `delta`, `omega` `Tgov1 ` | the TGOV1 governor model | `pmech`, `speed` | `R`, `T1`, `T2`, `T3`, `Pvmax`, `Pvmin`, `Dt` | `none` `Ieeet1` | the IEEET1 exciter model | `bus`, `speed`, `efd`, `vs`\* | `Tr`, `Ka`, `Ta`, `Ke`, `Te`, `Kf`, `Tf`, `Vrmin`, `Vrmax`, `E1`, `E2`, `Se1`, `Se2`, `Ispdlim` | `efd`, `ksat` diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/CMakeLists.txt b/GridKit/Model/PhasorDynamics/SynchronousMachine/CMakeLists.txt index 96c3b2f94..331c656b1 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/CMakeLists.txt @@ -6,4 +6,5 @@ # ]] add_subdirectory(GENROUwS) +add_subdirectory(GENSALwS) add_subdirectory(GenClassical) diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp index e764c5142..8eb942d98 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/GenrouImpl.hpp @@ -401,7 +401,8 @@ namespace GridKit ScalarT ii = (p * vi - q * vr) / vm2; // Initial ksat guess from |V| - ScalarT ksat = SB_ * (vm - SA_) * (vm - SA_); + ScalarT vm_sat = vm - SA_; + ScalarT ksat = (vm_sat > ZERO) ? SB_ * vm_sat * vm_sat : ScalarT{ZERO}; ScalarT delta, id, iq, vd, vq; ScalarT psiqpp, psidpp, psipp; @@ -428,7 +429,7 @@ namespace GridKit id = ir * std::sin(delta) - ii * std::cos(delta); iq = ir * std::cos(delta) + ii * std::sin(delta); vd = vr * std::sin(delta) - vi * std::cos(delta) + id * Ra_ - iq * Xqpp_; - vq = vr * std::cos(delta) + vi * std::sin(delta) + id * Xqpp_ - iq * Ra_; + vq = vr * std::cos(delta) + vi * std::sin(delta) + id * Xqpp_ + iq * Ra_; psiqpp = -vd; psidpp = vq; Edp = (Xq1_ - Xqd_ * (Xqp_ - Xqpp_) * ksat) * iq / (ONE + Xqd_ * ksat); @@ -440,7 +441,8 @@ namespace GridKit ScalarT psiqpp_fl = -psiqp * Xq4_ - Edp * Xq5_; ScalarT psidpp_fl = psidp * Xd4_ + Eqp * Xd5_; psipp = std::sqrt(psiqpp_fl * psiqpp_fl + psidpp_fl * psidpp_fl); - ksat = SB_ * (psipp - SA_) * (psipp - SA_); + ScalarT psipp_sat = psipp - SA_; + ksat = (psipp_sat > ZERO) ? SB_ * psipp_sat * psipp_sat : ScalarT{ZERO}; if (iter == max_iter - 1) { @@ -460,8 +462,9 @@ namespace GridKit y_[5] = Edp; y_[6] = psiqpp = -psiqp * Xq4_ - Edp * Xq5_; y_[7] = psidpp = psidp * Xd4_ + Eqp * Xd5_; - y_[8] = psipp = std::sqrt(psiqpp * psiqpp + psidpp * psidpp); - y_[9] = ksat = SB_ * ((psipp - SA_) * (psipp - SA_)); + y_[8] = psipp = std::sqrt(psiqpp * psiqpp + psidpp * psidpp); + ScalarT psipp_sat = psipp - SA_; + y_[9] = ksat = (psipp_sat > ZERO) ? SB_ * psipp_sat * psipp_sat : ScalarT{ZERO}; y_[10] = vd = -psiqpp * (ONE + omega); y_[11] = vq = psidpp * (ONE + omega); y_[12] = (psidpp - id * Xdpp_) * iq - (psiqpp - iq * Xdpp_) * id; @@ -566,17 +569,18 @@ namespace GridKit f[5] = Edp_dot - (ONE / Tqop_) * (-Edp + Xqd_ * psiqpp * ksat + Xq1_ * (iq - Xq3_ * (Edp + iq * Xq2_ - psiqp))); /* 11 Genrou algebraic equations */ - f[6] = psiqpp - (-psiqp * Xq4_ - Edp * Xq5_); - f[7] = psidpp - (psidp * Xd4_ + Eqp * Xd5_); - f[8] = psipp - std::sqrt((psidpp * psidpp) + (psiqpp * psiqpp)); - f[9] = ksat - SB_ * ((psipp - SA_) * (psipp - SA_)); - f[10] = vd + psiqpp * (ONE + omega); - f[11] = vq - psidpp * (ONE + omega); - f[12] = telec - ((psidpp - id * Xdpp_) * iq - (psiqpp - iq * Xdpp_) * id); - f[13] = id - (ir * std::sin(delta) - ii * std::cos(delta)); - f[14] = iq - (ir * std::cos(delta) + ii * std::sin(delta)); - f[15] = ir + G_ * vr - B_ * vi - inr; - f[16] = ii + B_ * vr + G_ * vi - ini; + f[6] = psiqpp - (-psiqp * Xq4_ - Edp * Xq5_); + f[7] = psidpp - (psidp * Xd4_ + Eqp * Xd5_); + f[8] = psipp - std::sqrt((psidpp * psidpp) + (psiqpp * psiqpp)); + ScalarT psipp_sat = psipp - SA_; + f[9] = ksat - SB_ * psipp_sat * psipp_sat * Math::sigmoid(psipp_sat); + f[10] = vd + psiqpp * (ONE + omega); + f[11] = vq - psidpp * (ONE + omega); + f[12] = telec - ((psidpp - id * Xdpp_) * iq - (psiqpp - iq * Xdpp_) * id); + f[13] = id - (ir * std::sin(delta) - ii * std::cos(delta)); + f[14] = iq - (ir * std::cos(delta) + ii * std::sin(delta)); + f[15] = ir + G_ * vr - B_ * vi - inr; + f[16] = ii + B_ * vr + G_ * vi - ini; /* 2 Genrou current source definitions */ f[17] = inr - (G_ * (std::sin(delta) * vd + std::cos(delta) * vq) - B_ * (-std::cos(delta) * vd + std::sin(delta) * vq)); diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/README.md b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/README.md index b065a4edc..95bf1eff6 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/README.md +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROUwS/README.md @@ -131,7 +131,7 @@ Note that for implementation purposes, some of these equations may be simplified 0 &= -V_{d} -\psi''_{q}(1+\omega)\\ 0 &= -V_{q} +\psi''_{d}(1+\omega)\\ 0 &= -T_{elec} +(\psi''_{d} - I_dX_d'')I_q-(\psi''_{q} - I_qX_d'')I_d \\ - 0 &= -k_{sat} + S_B (\psi''-S_A)^2 \\ + 0 &= -k_{sat} + S_B(\psi''-S_A)^2\sigma(\psi''-S_A) \\ 0 &= -I_d + I_r \sin(\delta) - I_i \cos(\delta) \\ 0 &= -I_q + I_r \cos(\delta) + I_i \sin(\delta) \\ 0 &= -I_r + G (V_d \sin(\delta) + V_q \cos(\delta) - V_r) - B (V_d \cos(\delta) + V_q \sin(\delta) - V_i) \\ @@ -156,7 +156,11 @@ from the steady-state initial conditions. \psi^{''}_{d} &= V_q \\ \psi^{''}_{q} &= -V_d \\ \psi^{''} &= \sqrt{(\psi''_{d})^2+(\psi''_{q})^2} \\ - k_{sat} &= S_B(\psi^{''}-S_A)^2 \\ + k_{sat} &= + \begin{cases} + S_B(\psi^{''}-S_A)^2, & \psi^{''} > S_A\\ + 0, & \psi^{''} \le S_A + \end{cases}\\ T_{elec} &= (\psi''_{d} - I_dX_d^{''})I_q-(\psi''_{q} - I_qX_d^{''})I_d \\ P_{m} &= T_{elec} \\ \psi_d' &= \psi_d'' - (X_d'' - X_\ell)I_d \\ diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/CMakeLists.txt b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/CMakeLists.txt new file mode 100644 index 000000000..c4f14f07e --- /dev/null +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/CMakeLists.txt @@ -0,0 +1,45 @@ +set(_install_headers + Gensal.hpp + GensalData.hpp) + +if(GRIDKIT_ENABLE_ENZYME) + gridkit_add_library(phasor_dynamics_gensal + SOURCES + GensalEnzyme.cpp + HEADERS + ${_install_headers} + INCLUDE_DIRECTORIES + PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC GridKit::phasor_dynamics_core + PUBLIC GridKit::phasor_dynamics_signal + PRIVATE ClangEnzymeFlags + COMPILE_OPTIONS + PRIVATE -mllvm -enzyme-auto-sparsity=1 -fno-math-errno -fno-vectorize) +else() + gridkit_add_library(phasor_dynamics_gensal + SOURCES + Gensal.cpp + HEADERS + ${_install_headers} + INCLUDE_DIRECTORIES + PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC GridKit::phasor_dynamics_core + PUBLIC GridKit::phasor_dynamics_signal) +endif() + +gridkit_add_library(phasor_dynamics_gensal_dependency_tracking + SOURCES + GensalDependencyTracking.cpp + INCLUDE_DIRECTORIES + PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC GridKit::phasor_dynamics_core + PUBLIC GridKit::phasor_dynamics_signal_dependency_tracking) + +# Link to interface target for all components +target_link_libraries(phasor_dynamics_components + INTERFACE GridKit::phasor_dynamics_gensal) +target_link_libraries(phasor_dynamics_components_dependency_tracking + INTERFACE GridKit::phasor_dynamics_gensal_dependency_tracking) diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/Gensal.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/Gensal.cpp new file mode 100644 index 000000000..4ced73ebd --- /dev/null +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/Gensal.cpp @@ -0,0 +1,32 @@ +/** + * @file Gensal.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Definition of a GENSAL generator model. + */ + +#include "GensalImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + /** + * @brief Jacobian evaluation not implemented yet + * + * @tparam ScalarT - scalar data type + * @tparam IdxT - matrix index data type + * @return int - error code, 0 = success + */ + template + int Gensal::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Gensal..." << std::endl; + Log::misc() << "Jacobian evaluation not implemented!" << std::endl; + return 0; + } + + // Available template instantiations + template class Gensal; + template class Gensal; + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/Gensal.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/Gensal.hpp new file mode 100644 index 000000000..dfab5b0f1 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/Gensal.hpp @@ -0,0 +1,204 @@ +/** + * @file Gensal.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Declaration of a GENSAL generator model. + * + */ + +#pragma once + +#include +#include +#include +#include + +// Forward declarations. +namespace GridKit +{ + namespace PhasorDynamics + { + template + class BusBase; + + template + struct GensalData; + } // namespace PhasorDynamics +} // namespace GridKit + +namespace GridKit +{ + namespace PhasorDynamics + { + /// Internal variables of a `Gensal` + enum class GensalInternalVariables : size_t + { + DELTA, ///< rotor angle + OMEGA, ///< speed deviation + EPQ, ///< q-axis transient voltage + PSIPD, ///< d-axis transient flux + PSIPPQ, ///< q-axis subtransient flux + PSIPPD, ///< d-axis subtransient flux + KSAT, ///< saturation signal + VD, ///< d-axis terminal voltage + VQ, ///< q-axis terminal voltage + TE, ///< electrical torque + ID, ///< d-axis current + IQ, ///< q-axis current + IR, ///< network real current + II, ///< network imaginary current + INR, ///< Norton source real current + INI, ///< Norton source imaginary current + MAXIMUM, + }; + + /// External variables of a `Gensal` + enum class GensalExternalVariables : size_t + { + VR, ///< network real voltage + VI, ///< network imaginary voltage + PM, ///< mechanical power + EFD, ///< field voltage + MAXIMUM, + }; + + template + class Gensal : public Component + { + using Component::gridkit_component_id_; + using Component::alpha_; + using Component::f_; + using Component::nnz_; + using Component::size_; + using Component::tag_; + using Component::time_; + using Component::y_; + using Component::yp_; + using Component::wb_; + using Component::h_; + using Component::J_; + using Component::J_rows_buffer_; + using Component::J_cols_buffer_; + using Component::J_vals_buffer_; + using Component::mva_system_base_; + using Component::variable_indices_; + using Component::residual_indices_; + + public: + using RealT = typename Component::RealT; + using bus_type = BusBase; + using model_data_type = GensalData; + using MonitorT = Model::VariableMonitor; + + Gensal(bus_type* bus, const model_data_type& data); + ~Gensal(); + + int setGridKitComponentID(IdxT) override final; + int allocate() override final; + int verify() const override final; + int initialize() override final; + int tagDifferentiable() override final; + int evaluateResidual() override final; + + // Still to be implemented + int evaluateJacobian() override final; + + // Temporary access functions for governor + // Should be abstracted + ScalarT getSpeed(); + ScalarT getTorque(); + + /// Get the `ComponentSignals` from this `Gensal` + auto getSignals() + -> ComponentSignals& + { + return signals_; + } + + const Model::VariableMonitorBase* getMonitor() const override; + + private: + void initializeParameters(const model_data_type& data); + /// Associate variable getter functions with enum values + void initializeMonitor(); + void setDerivedParams(); + + ScalarT& Vr() + { + return bus_->Vr(); + } + + ScalarT& Vi() + { + return bus_->Vi(); + } + + ScalarT& Ir() + { + return bus_->Ir(); + } + + ScalarT& Ii() + { + return bus_->Ii(); + } + + public: + __attribute__((always_inline)) inline int evaluateInternalResidual(ScalarT*, ScalarT*, ScalarT*, ScalarT*, ScalarT*); + __attribute__((always_inline)) inline int evaluateBusResidual(ScalarT*, ScalarT*, ScalarT*, ScalarT*); + + private: + /* Identification */ + bus_type* bus_; + + /// Component signal extension + ComponentSignals signals_; + + /* Initial terminal conditions */ + RealT p0_{0.0}; + RealT q0_{0.0}; + + /* Input parameters */ + RealT H_{3.0}; + RealT D_{0.0}; + RealT Ra_{0.0}; + RealT Tdop_{7.0}; + RealT Tdopp_{0.04}; + RealT Tqopp_{0.05}; + RealT Xd_{2.1}; + RealT Xdp_{0.2}; + RealT Xdpp_{0.18}; + RealT Xq_{0.5}; + RealT Xl_{0.15}; + RealT S10_{0.0}; + RealT S12_{0.0}; + RealT mva_base_{100.0}; + + /* Derived parameters */ + RealT SA_; + RealT SB_; + RealT Xd1_; + RealT Xd2_; + RealT Xd3_; + RealT Xd4_; + RealT Xd5_; + RealT Xq2_; + RealT G_; + RealT B_; + + /* Setpoints for control variables (determined at initialization) */ + ScalarT pmech_set_{0.0}; // TODO remove default initialization and ensure this gets set + ScalarT efd_set_{0.0}; // TODO remove default initialization and ensure this gets set + + /* Local copies of signal variables */ + std::vector ws_; + std::vector ws_indices_; + + /// Variable monitor + std::unique_ptr monitor_; + }; + + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalData.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalData.hpp new file mode 100644 index 000000000..53585ad1b --- /dev/null +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalData.hpp @@ -0,0 +1,79 @@ +/** + * @file GensalData.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Modeling data for a GENSAL generator model. + * + */ +#pragma once + +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + /// Initial parameters for a Gensal generator model + enum class GensalParameters + { + p0, ///< Initial active power + q0, ///< Initial reactive power + H, ///< Rotor inertia + D, ///< Damping coefficient + Ra, ///< Winding resistance + Tdop, ///< Open circuit direct axis transient time + Tdopp, ///< Open circuit direct axis sub-transient time + Tqopp, ///< Open circuit quadrature axis sub-transient time + Xd, ///< Direct axis synchronous reactance + Xdp, ///< Direct axis transient reactance + Xdpp, ///< Direct axis sub-transient reactance + Xq, ///< Quadrature axis synchronous reactance + Xl, ///< Stator leakage reactance + S10, ///< Saturation factor at 1.0 pu flux + S12, ///< Saturation factor at 1.2 pu flux + mva_base, ///< MVA base of the gensal model + }; + + /// Ports for a Gensal generator model + enum class GensalPorts + { + bus, ///< Unique ID of the connecting bus + pmech, ///< Unique ID of the signal providing mechanical power + speed, ///< Unique ID of the signal receiving speed deviation + efd, ///< Unique ID of the signal providing exciter field voltage + }; + + /// Variables able to be monitored for a Gensal generator model + enum class GensalMonitorableVariables + { + ir, + ii, + p, + q, + delta, + omega, + speed + }; + + /** + * @brief Contains modeling data for a Gensal generator model. + * + * @tparam RealT Real parameter data type + * @tparam IdxT Integer parameter data type + * + * Integer parameters are of the same type as matrix and vector indices. + */ + template + struct GensalData : public ComponentData + { + GensalData() = default; + + using Parameters = GensalParameters; + using Ports = GensalPorts; + using MonitorableVariables = GensalMonitorableVariables; + }; + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalDependencyTracking.cpp new file mode 100644 index 000000000..3bff58954 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalDependencyTracking.cpp @@ -0,0 +1,26 @@ +#include "GensalImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + /** + * @brief Jacobian evaluation not implemented yet + * + * @tparam ScalarT - scalar data type + * @tparam IdxT - matrix index data type + * @return int - error code, 0 = success + */ + template + int Gensal::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Gensal..." << std::endl; + Log::misc() << "Jacobian evaluation not implemented!" << std::endl; + return 0; + } + + // Available template instantiations + template class Gensal; + template class Gensal; + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalEnzyme.cpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalEnzyme.cpp new file mode 100644 index 000000000..a70625e1c --- /dev/null +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalEnzyme.cpp @@ -0,0 +1,130 @@ +/** + * @file GensalEnzyme.cpp + * @author Luke Lowery (lukel@tamu.edu) + * + */ + +#include + +#include "GensalImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + /** + * @brief Jacobian evaluation experimental + * + * @tparam ScalarT - scalar data type + * @tparam IdxT - matrix index data type + * @return int - error code, 0 = success + */ + template + int Gensal::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Gensal..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + + J_.zeroMatrix(); + if (J_rows_buffer_ == nullptr) + { + // Reserve space for a dense matrix of size_*size_. + // Enzyme will compute the appropriate nnz from sparsification. + J_rows_buffer_ = new IdxT[static_cast(size_) * static_cast(size_)]; + J_cols_buffer_ = new IdxT[static_cast(size_) * static_cast(size_)]; + J_vals_buffer_ = new RealT[static_cast(size_) * static_cast(size_)]; + } + + GridKit::Enzyme::Sparse::DfDy, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, + ScalarT, + IdxT>::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + J_); + + GridKit::Enzyme::Sparse::DfDwb, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, + ScalarT, + IdxT>::eval(this, + f_.size(), + static_cast(bus_->size()), + (this->getResidualIndices()).data(), + (bus_->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + J_); + + GridKit::Enzyme::Sparse::DfDws, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal, + ScalarT, + IdxT>::eval(this, + f_.size(), + ws_.size(), + (this->getResidualIndices()).data(), + ws_indices_.data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + J_); + + GridKit::Enzyme::Sparse::DhDwb, + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, + ScalarT, + IdxT>::eval(this, + static_cast(bus_->size()), + static_cast(bus_->size()), + (bus_->getResidualIndices()).data(), + (bus_->getVariableIndices()).data(), + y_.data(), + yp_.data(), + (bus_->y()).data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + bus_->getJacobian()); + + GridKit::Enzyme::Sparse::DhDy, + GridKit::Enzyme::Sparse::MemberFunctions::BusResidual, + ScalarT, + IdxT>::eval(this, + static_cast(bus_->size()), + y_.size(), + (bus_->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + J_); + + return 0; + } + + // Available template instantiations + template class Gensal; + template class Gensal; + + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalImpl.hpp new file mode 100644 index 000000000..f0c6aaf65 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/GensalImpl.hpp @@ -0,0 +1,496 @@ +#pragma once + +#define _USE_MATH_DEFINES +#include +#include + +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + using Log = ::GridKit::Utilities::Logger; + + /** + * @brief Constructor for a GENSAL generator model with saturation + */ + template + Gensal::Gensal(bus_type* bus, const model_data_type& data) + : bus_(bus), + monitor_(std::make_unique(data)) + { + initializeParameters(data); + initializeMonitor(); + + size_ = 16; + setDerivedParams(); + } + + template + Gensal::~Gensal() + { + } + + /// Helper function to extract and assign model parameters from the model's associated + /// data structure. + template + void Gensal::initializeParameters(const model_data_type& data) + { + if (data.parameters.contains(model_data_type::Parameters::p0)) + { + p0_ = std::get(data.parameters.at(model_data_type::Parameters::p0)); + } + + if (data.parameters.contains(model_data_type::Parameters::q0)) + { + q0_ = std::get(data.parameters.at(model_data_type::Parameters::q0)); + } + + if (data.parameters.contains(model_data_type::Parameters::H)) + { + H_ = std::get(data.parameters.at(model_data_type::Parameters::H)); + } + + if (data.parameters.contains(model_data_type::Parameters::D)) + { + D_ = std::get(data.parameters.at(model_data_type::Parameters::D)); + } + + if (data.parameters.contains(model_data_type::Parameters::Ra)) + { + Ra_ = std::get(data.parameters.at(model_data_type::Parameters::Ra)); + } + + if (data.parameters.contains(model_data_type::Parameters::Tdop)) + { + Tdop_ = std::get(data.parameters.at(model_data_type::Parameters::Tdop)); + } + + if (data.parameters.contains(model_data_type::Parameters::Tdopp)) + { + Tdopp_ = std::get(data.parameters.at(model_data_type::Parameters::Tdopp)); + } + + if (data.parameters.contains(model_data_type::Parameters::Tqopp)) + { + Tqopp_ = std::get(data.parameters.at(model_data_type::Parameters::Tqopp)); + } + + if (data.parameters.contains(model_data_type::Parameters::Xd)) + { + Xd_ = std::get(data.parameters.at(model_data_type::Parameters::Xd)); + } + + if (data.parameters.contains(model_data_type::Parameters::Xdp)) + { + Xdp_ = std::get(data.parameters.at(model_data_type::Parameters::Xdp)); + } + + if (data.parameters.contains(model_data_type::Parameters::Xdpp)) + { + Xdpp_ = std::get(data.parameters.at(model_data_type::Parameters::Xdpp)); + } + + if (data.parameters.contains(model_data_type::Parameters::Xq)) + { + Xq_ = std::get(data.parameters.at(model_data_type::Parameters::Xq)); + } + + if (data.parameters.contains(model_data_type::Parameters::Xl)) + { + Xl_ = std::get(data.parameters.at(model_data_type::Parameters::Xl)); + } + + if (data.parameters.contains(model_data_type::Parameters::S10)) + { + S10_ = std::get(data.parameters.at(model_data_type::Parameters::S10)); + } + + if (data.parameters.contains(model_data_type::Parameters::S12)) + { + S12_ = std::get(data.parameters.at(model_data_type::Parameters::S12)); + } + + if (data.parameters.contains(model_data_type::Parameters::mva_base)) + { + mva_base_ = std::get(data.parameters.at(model_data_type::Parameters::mva_base)); + } + } + + template + const Model::VariableMonitorBase* Gensal::getMonitor() const + { + return monitor_.get(); + } + + template + void Gensal::initializeMonitor() + { + using Variable = typename model_data_type::MonitorableVariables; + monitor_->set(Variable::ir, [this] + { return y_[12]; }); + monitor_->set(Variable::ii, [this] + { return y_[13]; }); + monitor_->set(Variable::p, [this] + { return Vr() * Ir() + Vi() * Ii(); }); + monitor_->set(Variable::q, [this] + { return Vi() * Ir() - Vr() * Ii(); }); + monitor_->set(Variable::delta, [this] + { return y_[0]; }); + monitor_->set(Variable::omega, [this] + { return y_[1]; }); + monitor_->set(Variable::speed, [this] + { return 1.0 + y_[1]; }); + } + + /** + * @brief Set the component ID + */ + template + int Gensal::setGridKitComponentID(IdxT component_id) + { + gridkit_component_id_ = component_id; + return 0; + } + + /*! + * @brief allocate method computes sparsity pattern of the Jacobian. + */ + template + int Gensal::allocate() + { + // Resize component model data + auto size = static_cast(size_); + f_.resize(size); + y_.resize(size); + yp_.resize(size); + tag_.resize(size); + variable_indices_.resize(size); + residual_indices_.resize(size); + + // Resize bus data + wb_.resize(2); + h_.resize(2); + + // Resize signal variable data + ws_.resize(2); + ws_indices_.resize(2); + ws_indices_[0] = INVALID_INDEX; + ws_indices_[1] = INVALID_INDEX; + + // Default variable and residual index mapping to local index + for (IdxT j = 0; j < size_; ++j) + { + this->setVariableIndex(j, j); + this->setResidualIndex(j, j); + } + + // Set output signals + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set(&y_[1], &(this->getVariableIndex(1))); + } + + return 0; + } + + /** + * @brief verify method checks that attached signals are also linked + */ + template + int Gensal::verify() const + { + static constexpr auto PM = GensalExternalVariables::PM; + static constexpr auto EFD = GensalExternalVariables::EFD; + + int ret = 0; + + if (signals_.template isAttached()) + { + if (!signals_.template isLinked()) + { + Log::error() << "Gensal: pmech signal attached with no linked governor\n"; + ret += 1; + } + } + + if (signals_.template isAttached()) + { + if (!signals_.template isLinked()) + { + Log::error() << "Gensal: efd signal attached with no linked exciter\n"; + ret += 1; + } + } + + return ret; + } + + /** + * Initialization of the generator model + * + */ + template + int Gensal::initialize() + { + // Network frame terminal values + ScalarT vr = Vr(); + ScalarT vi = Vi(); + ScalarT p = static_cast(p0_) * mva_system_base_ / mva_base_; + ScalarT q = static_cast(q0_) * mva_system_base_ / mva_base_; + ScalarT vm2 = vr * vr + vi * vi; + ScalarT ir = (p * vr + q * vi) / vm2; + ScalarT ii = (p * vi - q * vr) / vm2; + + ScalarT Er = vr + Ra_ * ir - Xq_ * ii; + ScalarT Ei = vi + Ra_ * ii + Xq_ * ir; + ScalarT delta = std::atan2(Ei, Er); + ScalarT omega(0.0); + + ScalarT id = ir * std::sin(delta) - ii * std::cos(delta); + ScalarT iq = ir * std::cos(delta) + ii * std::sin(delta); + ScalarT psiqpp = -Xq2_ * iq; + ScalarT vd = -psiqpp * (ONE + omega); + ScalarT vq = vr * std::cos(delta) + vi * std::sin(delta) + id * Xdpp_ + iq * Ra_; + ScalarT psidpp = vq / (ONE + omega); + ScalarT psidp = psidpp - (Xdpp_ - Xl_) * id; + ScalarT Eqp = psidp + Xd2_ * id; + ScalarT Eqp_sat = Eqp - SA_; + ScalarT ksat = SB_ * Eqp_sat * Eqp_sat * Math::sigmoid(Eqp_sat); + ScalarT Te = (psidpp - id * Xdpp_) * iq - (psiqpp - iq * Xdpp_) * id; + + y_[0] = delta; + y_[1] = omega; + y_[2] = Eqp; + y_[3] = psidp; + y_[4] = psiqpp; + y_[5] = psidpp; + y_[6] = ksat; + y_[7] = vd; + y_[8] = vq; + y_[9] = Te; + y_[10] = id; + y_[11] = iq; + y_[12] = ir; + y_[13] = ii; + y_[14] = G_ * (vd * std::sin(delta) + vq * std::cos(delta)) + - B_ * (vd * -std::cos(delta) + vq * std::sin(delta)); + y_[15] = B_ * (vd * std::sin(delta) + vq * std::cos(delta)) + + G_ * (vd * -std::cos(delta) + vq * std::sin(delta)); + + pmech_set_ = Te; + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(Te); + } + + efd_set_ = Eqp + Xd1_ * (id + Xd3_ * (Eqp - psidp - Xd2_ * id)) + ksat; + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(efd_set_); + } + + for (IdxT i = 0; i < size_; ++i) + { + yp_[static_cast(i)] = 0.0; + } + + return 0; + } + + /** + * \brief Identify differential variables. + */ + template + int Gensal::tagDifferentiable() + { + for (IdxT i = 0; i < size_; ++i) + { + tag_[static_cast(i)] = i < 5; + } + return 0; + } + + /** + * @brief Internal residual + * + */ + template + __attribute__((always_inline)) inline int Gensal::evaluateInternalResidual( + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + ScalarT* ws, + ScalarT* f) + { + /* Read variables */ + ScalarT delta = y[0]; + ScalarT omega = y[1]; + ScalarT Eqp = y[2]; + ScalarT psidp = y[3]; + ScalarT psiqpp = y[4]; + ScalarT psidpp = y[5]; + ScalarT ksat = y[6]; + ScalarT vd = y[7]; + ScalarT vq = y[8]; + ScalarT telec = y[9]; + ScalarT id = y[10]; + ScalarT iq = y[11]; + ScalarT ir = y[12]; + ScalarT ii = y[13]; + ScalarT inr = y[14]; + ScalarT ini = y[15]; + + /* Read derivatives */ + ScalarT delta_dot = yp[0]; + ScalarT omega_dot = yp[1]; + ScalarT Eqp_dot = yp[2]; + ScalarT psidp_dot = yp[3]; + ScalarT psiqpp_dot = yp[4]; + + // Set coupling variable aliases + ScalarT vr = wb[0]; + ScalarT vi = wb[1]; + + // Set signal variable aliases + ScalarT pmech = ws[0]; + ScalarT efd = ws[1]; + + /* 5 Gensal differential equations */ + // TODO: Replace hard-coded 60 Hz with the model's resolved frequency base + // once frequency-base ownership/propagation is designed consistently. + f[0] = delta_dot - omega * (TWO * M_PI * 60.0); + f[1] = omega_dot - (ONE / (TWO * H_)) * ((pmech - D_ * omega) / (ONE + omega) - telec); + f[2] = Eqp_dot - (ONE / Tdop_) * (efd - (Eqp + Xd1_ * (id + Xd3_ * (Eqp - psidp - Xd2_ * id)) + ksat)); + f[3] = psidp_dot - (ONE / Tdopp_) * (Eqp - psidp - Xd2_ * id); + f[4] = psiqpp_dot - (ONE / Tqopp_) * (-psiqpp - Xq2_ * iq); + + /* 9 Gensal algebraic equations */ + f[5] = psidpp - (psidp * Xd4_ + Eqp * Xd5_); + ScalarT Eqp_sat = Eqp - SA_; + f[6] = ksat - SB_ * Eqp_sat * Eqp_sat * Math::sigmoid(Eqp_sat); + f[7] = vd + psiqpp * (ONE + omega); + f[8] = vq - psidpp * (ONE + omega); + f[9] = telec - ((psidpp - id * Xdpp_) * iq - (psiqpp - iq * Xdpp_) * id); + f[10] = id - (ir * std::sin(delta) - ii * std::cos(delta)); + f[11] = iq - (ir * std::cos(delta) + ii * std::sin(delta)); + f[12] = ir + G_ * vr - B_ * vi - inr; + f[13] = ii + B_ * vr + G_ * vi - ini; + + /* 2 Gensal current source definitions */ + f[14] = inr - (G_ * (std::sin(delta) * vd + std::cos(delta) * vq) - B_ * (-std::cos(delta) * vd + std::sin(delta) * vq)); + f[15] = ini - (B_ * (std::sin(delta) * vd + std::cos(delta) * vq) + G_ * (-std::cos(delta) * vd + std::sin(delta) * vq)); + + return 0; + } + + /** + * @brief Bus residual + * + */ + template + __attribute__((always_inline)) inline int Gensal::evaluateBusResidual( + ScalarT* y, + [[maybe_unused]] ScalarT* yp, + ScalarT* wb, + ScalarT* h) + { + ScalarT inr = y[14]; + ScalarT ini = y[15]; + ScalarT vr = wb[0]; + ScalarT vi = wb[1]; + + h[0] = (inr - vr * G_ + vi * B_) * mva_base_ / mva_system_base_; + h[1] = (ini - vr * B_ - vi * G_) * mva_base_ / mva_system_base_; + + return 0; + } + + /** + * \brief Residual evaluation and contribution to the connected bus + * + */ + template + int Gensal::evaluateResidual() + { + // Mechanical Power + ws_[0] = pmech_set_; + if (signals_.template isAttached()) + { + ws_[0] = signals_.template readExternalVariable(); + ws_indices_[0] = signals_.template readExternalVariableIndex(); + } + + // Exciter Efield + ws_[1] = efd_set_; + if (signals_.template isAttached()) + { + ws_[1] = signals_.template readExternalVariable(); + ws_indices_[1] = signals_.template readExternalVariableIndex(); + } + + // Bus voltages + wb_[0] = Vr(); + wb_[1] = Vi(); + + // Residual evaluation + evaluateInternalResidual(y_.data(), yp_.data(), wb_.data(), ws_.data(), f_.data()); + evaluateBusResidual(y_.data(), yp_.data(), wb_.data(), h_.data()); + + // Gensal contribution to bus algebraic equations + Ir() += h_[0]; + Ii() += h_[1]; + + return 0; + } + + /** + * @brief Access generator relative speed + * + * @tparam ScalarT - scalar data type + * @tparam IdxT - matrix index data type + * @return int - error code, 0 = success + */ + template + ScalarT Gensal::getSpeed() + { + return y_[1]; + } + + template + ScalarT Gensal::getTorque() + { + return y_[9]; + } + + template + void Gensal::setDerivedParams() + { + SA_ = 0; + SB_ = 0; + if (S12_ != 0) + { + RealT s112 = std::sqrt(S10_ / S12_); + + SA_ = (1.2 * s112 + ONE) / (s112 + ONE); + SB_ = (1.2 * s112 - ONE) / (s112 - ONE); + if (SB_ < SA_) + SA_ = SB_; + SB_ = S12_ / ((SA_ - 1.2) * (SA_ - 1.2)); + } + Xd1_ = Xd_ - Xdp_; + Xd2_ = Xdp_ - Xl_; + Xd3_ = (Xdp_ - Xdpp_) / (Xd2_ * Xd2_); + Xd4_ = (Xdp_ - Xdpp_) / Xd2_; + Xd5_ = (Xdpp_ - Xl_) / Xd2_; + Xq2_ = Xq_ - Xdpp_; + G_ = Ra_ / (Ra_ * Ra_ + Xdpp_ * Xdpp_); + B_ = -Xdpp_ / (Ra_ * Ra_ + Xdpp_ * Xdpp_); + } + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/README.md b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/README.md index 92b6d5a24..02aec5818 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/README.md +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSALwS/README.md @@ -1,55 +1,54 @@ -# GENSAL +# GENSAL -This synchronous machine model is 5th order and is specifically designed for salient pole machines. It is a standard model used in phasor-domain industry stability studies. -See the [General Synchronous Machine Model](../README.md) for general synchronous machine information. +This synchronous machine model is 5th order and is specifically designed for +salient-pole machines. It is a standard model used in phasor-domain industry +stability studies. +See the [General Synchronous Machine Model](../README.md) for general +synchronous machine information. Notes: -- $X''_{q}=X''_{d}$ (no subtransient saliency) -- $X''_{d}$ does not saturate -- Only d-axis affected by saturation -- $X_{q}=X'_{q}$ +- $X_q''=X_d''$ (no subtransient saliency) +- $X_q=X_q'$ - $T'_{q0}$ is neglected +- Only d-axis affected by saturation ## Block Diagram
- - - Figure 2: GENSAL. Figure courtesy of + + Figure 2: GENSAL. Figure courtesy of [PowerWorld](https://www.powerworld.com/WebHelp/)
## Model Parameters -Symbol | Units | Description | Typical Value | Note -------------|---------|---------------------------------|---------------| ------ -$\omega_0$ | [rad/s] | synchronous frequency | $2\pi60$ -$H$ | [s] | rotor inertia | 3 -$D$ | [p.u.] | damping coefficient | 0 -$R_a$ | [p.u.] | winding resistance | 0 -$X_{\ell}$ | [p.u.] | Stator leakage reactance | 0.15 | -$X_{d}$ | [p.u.] | Direct axis synchronous reactance | 2.1 | -$X'_{d}$ | [p.u.] | Direct axis transient reactance | 0.2 | -$X''_{d}$ | [p.u.] | Direct axis sub-transient reactance | 0.18 | -$X_{q}$ | [p.u.] | Quadrature axis synchronous reactance | 0.5 | -$T'_{d0}$ | [s] | Open circuit direct axis transient time const. | 7 | -$T''_{d0}$ | [s] | Open circuit direct axis sub-transient time const. | 0.04 | -$T''_{q0}$ | [s] | Open circuit quadrature axis sub-transient time const. | 0.05 | -$S_{10}$ | [p.u.] | Saturation factor at 1.0 pu flux | 0 | -$S_{12}$ | [p.u.] | Saturation factor at 1.2 pu flux | 0 | +Symbol | Units | Description | Typical Value | Note +------------|---------|--------------------------------------------|---------------| ------ +$\omega_0$ | [rad/s] | synchronous frequency | $2\pi \cdot 60$ +$H$ | [s] | rotor inertia | 3 +$D$ | [p.u.] | damping coefficient | 0 +$R_a$ | [p.u.] | winding resistance | 0 +$X_{\ell}$ | [p.u.] | Stator leakage reactance | 0.15 | +$X_d$ | [p.u.] | Direct axis synchronous reactance | 2.1 | +$X'_d$ | [p.u.] | Direct axis transient reactance | 0.2 | +$X''_d$ | [p.u.] | Direct axis sub-transient reactance | 0.18 | +$X_q$ | [p.u.] | Quadrature axis synchronous reactance | 0.5 | +$T'_{d0}$ | [s] | Open circuit direct axis transient time const. | 7 | +$T''_{d0}$ | [s] | Open circuit direct axis sub-transient time const. | 0.04 | +$T''_{q0}$ | [s] | Open circuit quadrature axis sub-transient time const. | 0.05 | +$S_{10}$ | [p.u.] | Saturation factor at 1.0 pu flux | 0 | +$S_{12}$ | [p.u.] | Saturation factor at 1.2 pu flux | 0 | ### Model Derived Parameters ``` math \begin{aligned} - G &=\dfrac{R_a}{R_a^2+(X''_q)^2}& - B &= -\dfrac{X''_q}{R_a^2+(X''_q)^2}\\ - S_A &= \dfrac{1.2\sqrt{S_{10}/S_{12}} +1}{\sqrt{S_{10}/S_{12}} +1} & - S_B &= \dfrac{1.2\sqrt{S_{10}/S_{12}} -1}{\sqrt{S_{10}/S_{12}} -1} \\ - X_{d1} &= X_d-X_d' \\ - X_{d2} &= X_d'-X_\ell & X_{q2} = X_q-X''_q \\ - X_{d3} &= (X_d'-X_d'')/X_{d2}^2 \\ - X_{d5} &= (X_d''-X_\ell)/X_{d2} \\ - X_{qd} &= (X_q-X_\ell)/(X_d-X_\ell) + G &= \dfrac{R_a}{R_a^2+(X_d'')^2} & + B &= -\dfrac{X_d''}{R_a^2+(X_d'')^2}\\ + S_A &= \dfrac{1.2\sqrt{S_{10}/S_{12}} +1}{\sqrt{S_{10}/S_{12}} +1} & + S_B &= \dfrac{1.2\sqrt{S_{10}/S_{12}} -1}{\sqrt{S_{10}/S_{12}} -1} \\ + X_{d1} &= X_d-X_d' & X_{q2} &= X_q-X_d'' \\ + X_{d2} &= X_d'-X_\ell & X_{d3} &= (X_d'-X_d'')/X_{d2}^2 \\ + X_{d4} &= (X_d'-X_d'')/X_{d2} & X_{d5} &= (X_d''-X_\ell)/X_{d2} \end{aligned} ``` @@ -59,25 +58,26 @@ $S_{12}$ | [p.u.] | Saturation factor at 1.2 pu flux | 0 | #### Differential -Symbol | Units | Description | Note -------------|---------|---------------------------------| ------ -$\delta$ | [rad] | Machine internal rotor angle | -$\omega$ | [p.u.] | Machine speed | Optionally read by governor or stabilizer component -$\psi'_d$ | [p.u.] | Direct axis subtransient flux | -$\psi''_q$ | [p.u.] | Quadrature axis subtransient flux | -$E'_q$ | [p.u.] | Quadrature axis subtransient flux | +Symbol | Units | Description | Note +-------------|--------|-----------------------------------|------- +$\delta$ | [rad] | Machine internal rotor angle | +$\omega$ | [p.u.] | Machine speed deviation | Optionally read by governor or stabilizer component +$E'_q$ | [p.u.] | Quadrature axis transient flux | +$\psi'_d$ | [p.u.] | Direct axis transient flux | +$\psi''_q$ | [p.u.] | Total q-axis subtransient flux | #### Algebraic -Symbol | Units | Description | Note -------------|---------|---------------------------------| ------ -$V_d$ | [p.u.] | Machine internal voltage, d-axis | -$V_q$ | [p.u.] | Machine internal voltage, q-axis | -$I_d$ | [p.u.] | Terminal current, d-axis | -$I_q$ | [p.u.] | Terminal current, q-axis | -$I_r$ | [p.u.] | Terminal current, real component on network reference frame | Read by bus and optionally by controllers -$I_i$ | [p.u.] | Terminal current, imaginary component on network reference frame | Read by bus and optionally by controllers -$\psi''_d$ | [p.u.] | Total d-axis subtransient flux -$T_{e}$ | [p.u.] | Electrical torque +Symbol | Units | Description | Note +------------|--------|-----------------------------------| ------ +$\psi''_d$ | [p.u.] | Total d-axis subtransient flux | +$k_{sat}$ | [p.u.] | Additive saturation signal | +$V_d$ | [p.u.] | Machine internal voltage, d-axis | +$V_q$ | [p.u.] | Machine internal voltage, q-axis | +$T_e$ | [p.u.] | Electrical torque | +$I_d$ | [p.u.] | Terminal current, d-axis | +$I_q$ | [p.u.] | Terminal current, q-axis | +$I_r$ | [p.u.] | Terminal current, real component on network reference frame | Read by bus and optionally by controllers +$I_i$ | [p.u.] | Terminal current, imaginary component on network reference frame | Read by bus and optionally by controllers ### External Variables @@ -85,43 +85,80 @@ $T_{e}$ | [p.u.] | Electrical torque None. #### Algebraic -Symbol | Units | Description | Note -------------|---------|---------------------------------| ------ -$V_r$ | [p.u.] | Terminal voltage, real component on network reference frame | owned by bus object -$V_i$ | [p.u.] | Terminal voltage, imaginary component on network reference frame | owned by bus object -$P_{m}$ | [p.u.] | Mechanical power from the prime mover | Owned by governor, constant if no governor is connected to the machine -$E_{fd}$ | [p.u.] | Field winding voltage from the excitation system | Owned by exciter, constant if no exciter is connected to the machine +Symbol | Units | Description | Note +---------|--------|---------------------------------------------------------| ------ +$V_r$ | [p.u.] | Terminal voltage, real component on network reference frame | owned by bus object +$V_i$ | [p.u.] | Terminal voltage, imaginary component on network reference frame | owned by bus object +$P_m$ | [p.u.] | Mechanical power from the prime mover | Owned by governor, constant if no governor is connected to the machine +$E_{fd}$ | [p.u.] | Field winding voltage from the excitation system | Owned by exciter, constant if no exciter is connected to the machine ## Model Equations ### Differential Equations ``` math \begin{aligned} - \dot\delta &= (\omega-1)\cdot\omega_0 \\ - \dot\omega &= \dfrac{1}{2H}\left(\dfrac{P_{mech}-D (\omega-1)}{\omega} - - T_{elec}\right)\\ - \dot{\psi}'_{d} &= \dfrac{1}{T''_{d0}}(E'_{q}-\psi'_{d}-X_{d2}I_{d})\\ - \dot{\psi}''_{q} &= \dfrac{1}{T''_{q0}}(-\psi''_{q}-X_{q2}I_{q})\\ - \dot{E}'_{q} &= \dfrac{1}{T'_{d0}} + \dot\delta &= \omega\cdot\omega_0 \\ + \dot\omega &= \dfrac{1}{2H}\left(\dfrac{P_m-D\omega}{1+\omega} + - T_e\right)\\ + \dot{E}'_q &= \dfrac{1}{T'_{d0}} \left( - E_{fd}-E'_{q}-X_{d1} - (I_{d}+X_{d3}(E'_{q}-\psi'_{d}-X_{d2}I_{d})) - - S_B (E'_q-S_A)^2 + E_{fd}-E'_q-X_{d1} + (I_d+X_{d3}(E'_q-\psi'_d-X_{d2}I_d)) + -k_{sat} \right)\\ + \dot{\psi}'_d &= \dfrac{1}{T''_{d0}}(E'_q-\psi'_d-X_{d2}I_d)\\ + \dot{\psi}''_q &= \dfrac{1}{T''_{q0}}(-\psi''_q-X_{q2}I_q) \end{aligned} ``` ### Algebraic Equations -Note that for implementation purposes, some of these equations may be simplified into functions and the internal variables eliminated. Nevertheless, for modeling clarity and conformance to typical practice, the full equations are given here. ``` math \begin{aligned} - 0 &= -V_{d} -\psi''_{q}\omega\\ - 0 &= -V_{q} +\psi''_{d}\omega\\ - 0 &= -I_d + I_r \sin(\delta) - I_i \cos(\delta) \\ - 0 &= -I_q + I_r \cos(\delta) + I_i \sin(\delta) \\ - 0 &= -I_r + G (V_d \sin(\delta) + V_q \cos(\delta) - V_r) - B (V_d \cos(\delta) + V_q \sin(\delta) - V_i) \\ - 0 &= -I_i + B (V_d \sin(\delta) + V_q \cos(\delta) - V_r) + G (V_d \cos(\delta) + V_q \sin(\delta) - V_i) \\ - 0 &= -\psi''_{d} +E'_{q}X_{d5} + \psi'_{d}X_{d4}\\ - 0 &= -T_{e} +(\psi''_{d} - I_dX_d'')I_q-(\psi''_{q} - I_qX_d'')I_d + 0 &= -\psi''_d + E'_qX_{d5}+\psi'_dX_{d4}\\ + 0 &= -k_{sat} + S_B(E'_q-S_A)^2\sigma(E'_q-S_A)\\ + 0 &= -V_d -\psi''_q(1+\omega)\\ + 0 &= -V_q +\psi''_d(1+\omega)\\ + 0 &= -T_e +(\psi''_d-I_dX_d'')I_q-(\psi''_q-I_qX_d'')I_d\\ + 0 &= -I_d + I_r \sin(\delta) - I_i \cos(\delta) \\ + 0 &= -I_q + I_r \cos(\delta) + I_i \sin(\delta) \\ + 0 &= -I_r + G (V_d \sin(\delta) + V_q \cos(\delta) - V_r) - B (-V_d \cos(\delta) + V_q \sin(\delta) - V_i) \\ + 0 &= -I_i + B (V_d \sin(\delta) + V_q \cos(\delta) - V_r) + G (-V_d \cos(\delta) + V_q \sin(\delta) - V_i) \end{aligned} ``` + +## Initialization + +Using the power-flow solution, initial currents are calculated from active and +reactive power injection. The remaining variables are initialized from the +steady-state GENSAL equations. + +``` math +\begin{aligned} + \omega &= 0 \\ + \delta &= \text{arg}\left[V_r+jV_i+(R_a+jX_q)(I_r+jI_i)\right]\\ + I_d &= I_r\sin(\delta)-I_i\cos(\delta)\\ + I_q &= I_r\cos(\delta)+I_i\sin(\delta)\\ + \psi''_q &= -X_{q2}I_q\\ + V_d &= -\psi''_q\\ + V_q &= V_r\cos(\delta)+V_i\sin(\delta)+X_d''I_d+R_aI_q\\ + \psi''_d &= V_q\\ + \psi'_d &= \psi''_d-(X_d''-X_\ell)I_d\\ + E'_q &= \psi'_d+X_{d2}I_d\\ + k_{sat} &= S_B(E'_q-S_A)^2\sigma(E'_q-S_A)\\ + T_e &= (\psi''_d-I_dX_d'')I_q-(\psi''_q-I_qX_d'')I_d\\ + P_m &= T_e\\ + E_{fd} &= E'_q+X_{d1}(I_d+X_{d3}(E'_q-\psi'_d-X_{d2}I_d))+k_{sat} +\end{aligned} +``` + +## Model Outputs + +Symbol | Units | Description | Note +-----------|--------|-----------------------------------|------ +$I_r$ | [p.u.] | Terminal current, real component on network reference frame | Oriented leaving the machine +$I_i$ | [p.u.] | Terminal current, imaginary component on network reference frame | Oriented leaving the machine +$P$ | [p.u.] | Active power, $V_rI_r+V_iI_i$ | Oriented leaving the machine +$Q$ | [p.u.] | Reactive power, $V_iI_r-V_rI_i$ | Oriented leaving the machine +$\delta$ | [rad] | Machine internal rotor angle | +$\omega$ | [p.u.] | Machine speed deviation | $\omega=0$ at synchronous speed +$speed$ | [p.u.] | Per-unit machine speed | $1+\omega$ diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/README.md b/GridKit/Model/PhasorDynamics/SynchronousMachine/README.md index 2017877fa..e7fd5fbc9 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/README.md +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/README.md @@ -1,13 +1,10 @@ # General Synchronous Machine Model -> [!NOTE] -> Only the GENROU and classical generator models have been implemented. - ## Convention
- + Figure 1: Synchronous Machine. Figure courtesy of [PowerWorld](https://www.powerworld.com/files/Synchronous-Machines.pdf/)
diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index 455bb54ad..eeb0bf438 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -185,6 +185,38 @@ namespace GridKit addComponent(gen); } + // Add GENSAL generators + for (const auto& gendata : data.gensal) + { + IdxT bus_index = 0; + if (gendata.ports.contains(GensalData::Ports::bus)) + { + bus_index = gendata.ports.at(GensalData::Ports::bus); + } + + auto* gen = new Gensal(getBus(bus_index), gendata); + + if (gendata.ports.contains(GensalData::Ports::speed)) + { + IdxT speed = gendata.ports.at(GensalData::Ports::speed); + gen->getSignals().template assignSignalNode(getSignal(speed)); + } + + if (gendata.ports.contains(GensalData::Ports::pmech)) + { + IdxT pmech = gendata.ports.at(GensalData::Ports::pmech); + gen->getSignals().template attachSignalNode(getSignal(pmech)); + } + + if (gendata.ports.contains(GensalData::Ports::efd)) + { + IdxT efd = gendata.ports.at(GensalData::Ports::efd); + gen->getSignals().template attachSignalNode(getSignal(efd)); + } + + addComponent(gen); + } + // Add classical generators for (const auto& gendata : data.genclassical) { diff --git a/GridKit/Model/PhasorDynamics/SystemModelData.hpp b/GridKit/Model/PhasorDynamics/SystemModelData.hpp index f1f181c19..855347e78 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelData.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelData.hpp @@ -17,6 +17,7 @@ #include #include #include +#include #include #include @@ -40,6 +41,7 @@ namespace GridKit using SexsPtiDataT = Exciter::SexsPtiData; using IeeestDataT = Stabilizer::IeeestData; using GenrouDataT = GenrouData; + using GensalDataT = GensalData; using GenClassicalDataT = GenClassicalData; using LoadDataT = LoadData; using LoadZIPDataT = LoadZIPData; @@ -87,6 +89,7 @@ namespace GridKit std::vector branch; ///< Branches within the model std::vector bus_fault; ///< Bus faults within the model std::vector genrou; ///< GENROU instances within the model + std::vector gensal; ///< GENSAL instances within the model std::vector genclassical; ///< Classical generator instances within the model std::vector load; ///< Loads within the model std::vector loadzip; ///< Loads within the model diff --git a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp index e536d1bdd..0a65042d5 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp @@ -106,6 +106,12 @@ namespace GridKit raw_component.get_to(genrou); sm.genrou.push_back(genrou); } + else if (kind == "Gensal") + { + typename SystemModelData::GensalDataT gensal; + raw_component.get_to(gensal); + sm.gensal.push_back(gensal); + } else if (kind == "GenClassical") { typename SystemModelData::GenClassicalDataT gen_classical; diff --git a/examples/PhasorDynamics/Large/Texas/README.md b/examples/PhasorDynamics/Large/Texas/README.md index c97efbc53..56ef96476 100644 --- a/examples/PhasorDynamics/Large/Texas/README.md +++ b/examples/PhasorDynamics/Large/Texas/README.md @@ -43,7 +43,7 @@ Have not been implemented in GridKit: - GGOV1 - HYGOV - IEEEG1 -- GENSAL +- GENSAL (In GridKit, not added to this case yet) - IEEEST (In GridKit, not added to this case yet) The following examples needs to be constructed with this case. diff --git a/tests/UnitTests/PhasorDynamics/CMakeLists.txt b/tests/UnitTests/PhasorDynamics/CMakeLists.txt index 89c5b499b..986c3a79a 100644 --- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt +++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt @@ -39,6 +39,14 @@ target_link_libraries(test_phasor_genrou GridKit::definitions GridKit::phasor_dynamics_bus_dependency_tracking GridKit::testing) +add_executable(test_phasor_gensal runGensalTests.cpp) +target_link_libraries(test_phasor_gensal GridKit::definitions + GridKit::phasor_dynamics_gensal + GridKit::phasor_dynamics_gensal_dependency_tracking + GridKit::phasor_dynamics_bus + GridKit::phasor_dynamics_bus_dependency_tracking + GridKit::testing) + add_executable(test_phasor_governortgov1 runGovernorTgov1Tests.cpp) target_link_libraries(test_phasor_governortgov1 GridKit::definitions GridKit::phasor_dynamics_genrou @@ -89,6 +97,7 @@ target_link_libraries(test_phasor_system_single_component GridKit::phasor_dynami add_test(NAME PhasorDynamicsBusTest COMMAND test_phasor_bus) add_test(NAME PhasorDynamicsBranchTest COMMAND test_phasor_branch) add_test(NAME PhasorDynamicsGenrouTest COMMAND test_phasor_genrou) +add_test(NAME PhasorDynamicsGensalTest COMMAND test_phasor_gensal) add_test(NAME PhasorDynamicsGovernorTgov1Test COMMAND test_phasor_governortgov1) add_test(NAME PhasorDynamicsExciterSexsPtiTest COMMAND test_phasor_exciter_sexspti) add_test(NAME PhasorDynamicsStabilizerIeeestTest COMMAND test_phasor_stabilizer_ieeest) @@ -104,6 +113,7 @@ install(TARGETS test_phasor_bus test_phasor_load test_phasor_loadzip test_phasor_genrou + test_phasor_gensal test_phasor_governortgov1 test_phasor_exciter_sexspti test_phasor_stabilizer_ieeest diff --git a/tests/UnitTests/PhasorDynamics/GenrouTests.hpp b/tests/UnitTests/PhasorDynamics/GenrouTests.hpp index 5dec7e78c..fc2476fff 100644 --- a/tests/UnitTests/PhasorDynamics/GenrouTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GenrouTests.hpp @@ -146,7 +146,7 @@ namespace GridKit 0.21, -0.07, -0.19223748416156686, - 1.8896749891587163, + 2.0, 1.4, 0.31, 2.211, diff --git a/tests/UnitTests/PhasorDynamics/GensalTests.hpp b/tests/UnitTests/PhasorDynamics/GensalTests.hpp new file mode 100644 index 000000000..8c6f16586 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/GensalTests.hpp @@ -0,0 +1,381 @@ +#define _USE_MATH_DEFINES +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + + template + class GensalTests + { + private: + using RealT = typename PhasorDynamics::Component::RealT; + using GensalDataT = PhasorDynamics::GensalData; + static constexpr ScalarT tol_ = 1000 * std::numeric_limits::epsilon(); + + static GensalDataT makeGensalData() + { + using Parameter = typename GensalDataT::Parameters; + using Port = typename GensalDataT::Ports; + + GensalDataT data; + data.device_class = "Gensal"; + data.disambiguation_string = "1"; + data.ports[Port::bus] = 1; + data.parameters[Parameter::p0] = RealT{1.0}; + data.parameters[Parameter::q0] = RealT{0.05013}; + data.parameters[Parameter::H] = RealT{3.0}; + data.parameters[Parameter::D] = RealT{0.0}; + data.parameters[Parameter::Ra] = RealT{0.0}; + data.parameters[Parameter::Tdop] = RealT{7.0}; + data.parameters[Parameter::Tdopp] = RealT{0.04}; + data.parameters[Parameter::Tqopp] = RealT{0.05}; + data.parameters[Parameter::Xd] = RealT{2.1}; + data.parameters[Parameter::Xdp] = RealT{0.2}; + data.parameters[Parameter::Xdpp] = RealT{0.18}; + data.parameters[Parameter::Xq] = RealT{0.5}; + data.parameters[Parameter::Xl] = RealT{0.15}; + data.parameters[Parameter::S10] = RealT{0.0}; + data.parameters[Parameter::S12] = RealT{0.0}; + + return data; + } + + public: + GensalTests() = default; + ~GensalTests() = default; + + TestOutcome constructor() + { + TestStatus success = true; + + auto* bus = new PhasorDynamics::Bus(1.0, 0.0); + auto data = makeGensalData(); + + PhasorDynamics::Component* machine = + new PhasorDynamics::Gensal(bus, data); + + success *= (machine != nullptr); + + if (machine) + { + delete machine; + } + delete bus; + + return success.report(__func__); + } + + /** + * @brief Checks residual evaluation at initialized steady state. + */ + TestOutcome residual() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + auto data = makeGensalData(); + PhasorDynamics::Gensal gen(&bus, data); + + bus.allocate(); + bus.initialize(); + bus.evaluateResidual(); + + gen.allocate(); + gen.initialize(); + gen.evaluateResidual(); + + const std::vector& f = gen.getResidual(); + for (const auto& f_val : f) + { + if (!isEqual(f_val, 0.0, tol_)) + { + success = false; + break; + } + } + + return success.report(__func__); + } + + /** + * @brief Checks initialized steady state with nonzero armature resistance. + */ + TestOutcome residual_nonzero_ra() + { + TestStatus success = true; + + using Parameter = typename GensalDataT::Parameters; + + auto data = makeGensalData(); + data.parameters[Parameter::p0] = RealT{0.8}; + data.parameters[Parameter::q0] = RealT{0.2}; + data.parameters[Parameter::Ra] = RealT{0.05}; + data.parameters[Parameter::S10] = RealT{0.0}; + data.parameters[Parameter::S12] = RealT{0.0}; + + PhasorDynamics::Bus bus(1.0, 0.1); + PhasorDynamics::Gensal gen(&bus, data); + + bus.allocate(); + bus.initialize(); + bus.evaluateResidual(); + + gen.allocate(); + gen.initialize(); + gen.evaluateResidual(); + + const std::vector& f = gen.getResidual(); + for (const auto& f_val : f) + { + if (!isEqual(f_val, 0.0, tol_)) + { + success = false; + break; + } + } + + return success.report(__func__); + } + + /** + * @brief Verifies residual equations against hard-coded values. + */ + TestOutcome hard_coded_residual() + { + TestStatus success = true; + + using Parameter = typename GensalDataT::Parameters; + + auto data = makeGensalData(); + data.parameters[Parameter::p0] = RealT{1.0}; + data.parameters[Parameter::q0] = RealT{1.0}; + data.parameters[Parameter::H] = RealT{0.5}; + data.parameters[Parameter::D] = RealT{-1.0}; + data.parameters[Parameter::Ra] = RealT{0.5}; + data.parameters[Parameter::Tdop] = RealT{2.0}; + data.parameters[Parameter::Tdopp] = RealT{4.0}; + data.parameters[Parameter::Tqopp] = RealT{5.0}; + data.parameters[Parameter::Xd] = RealT{2.0}; + data.parameters[Parameter::Xdp] = RealT{1.0}; + data.parameters[Parameter::Xdpp] = RealT{0.5}; + data.parameters[Parameter::Xq] = RealT{1.5}; + data.parameters[Parameter::Xl] = RealT{0.25}; + data.parameters[Parameter::S10] = RealT{0.0}; + data.parameters[Parameter::S12] = RealT{0.0}; + + ScalarT Vr1{1.0}; + ScalarT Vi1{1.0}; + + PhasorDynamics::Bus bus(Vr1, Vi1); + PhasorDynamics::Gensal gen(&bus, data); + + const std::vector res_answer = { + 0.0, + 0.0, + 2.1083333333333334, + -1.028125, + 0.65, + 0.0, + 0.2, + -1.1, + -1.4, + 1.8125, + 0.5, + 0.25, + 2.65, + -0.05, + 0.3, + -1.2}; + + bus.allocate(); + bus.initialize(); + + gen.allocate(); + + gen.y()[0] = M_PI; // delta + gen.y()[1] = 1.0; // omega + gen.y()[2] = 2.0; // Eqp + gen.y()[3] = 0.5; // psidp + gen.y()[4] = -0.75; // psiqpp + gen.y()[5] = 1.0; // psidpp + gen.y()[6] = 0.2; // ksat + gen.y()[7] = 0.4; // vd + gen.y()[8] = 0.6; // vq + gen.y()[9] = 1.5; // telec + gen.y()[10] = 0.25; // id + gen.y()[11] = -0.5; // iq + gen.y()[12] = 0.75; // ir + gen.y()[13] = -0.25; // ii + gen.y()[14] = 0.1; // inr + gen.y()[15] = -0.2; // ini + + gen.yp()[0] = 2 * M_PI * 60.0; // delta_dot + gen.yp()[1] = -1.0; // omega_dot + gen.yp()[2] = 0.3; // Eqp_dot + gen.yp()[3] = -0.7; // psidp_dot + gen.yp()[4] = 0.9; // psiqpp_dot + + gen.evaluateResidual(); + std::vector& residual = gen.getResidual(); + + for (size_t i = 0; i < res_answer.size(); ++i) + { + if (!isEqual(residual[i], res_answer[i], tol_)) + { + std::cout << "Incorrect result for residual " << i << ": " + << residual[i] << " != " << res_answer[i] << "\n"; + success = false; + break; + } + } + + return success.report(__func__); + } + +#ifdef GRIDKIT_ENABLE_ENZYME + /** + * @brief Checks Jacobian evaluation. + */ + TestOutcome jacobian() + { + TestStatus success = true; + + auto tol = 10 * std::numeric_limits::epsilon(); + + std::vector dependency_tracking_jacobian = DependencyTrackingJacobian(); + std::vector enzyme_jacobian = EnzymeJacobian(); + + for (size_t i = 0; i < dependency_tracking_jacobian.size(); ++i) + { + success *= (GridKit::Testing::isEqual(dependency_tracking_jacobian[i], enzyme_jacobian[i], tol)); + } + + return success.report(__func__); + } + + private: + std::vector DependencyTrackingJacobian() + { + DependencyTracking::Variable Vr1{1.0}; + DependencyTracking::Variable Vi1{0.0}; + PhasorDynamics::Bus bus(Vr1, Vi1); + auto data = makeGensalData(); + PhasorDynamics::Gensal gen(&bus, data); + + bus.allocate(); + gen.allocate(); + + bus.initialize(); + gen.initialize(); + + for (size_t i = 0; i < gen.size(); ++i) + { + gen.y()[i].setVariableNumber(i); + } + for (size_t i = 0; i < bus.size(); ++i) + { + bus.y()[i].setVariableNumber(i + gen.size()); + } + + bus.evaluateResidual(); + gen.evaluateResidual(); + std::vector residual_y = gen.getResidual(); + + bus.initialize(); + gen.initialize(); + + for (size_t i = 0; i < gen.size(); ++i) + { + gen.yp()[i].setVariableNumber(i); + } + + bus.evaluateResidual(); + gen.evaluateResidual(); + std::vector residual_yp = gen.getResidual(); + + std::vector dependencies(residual_y.size()); + for (IdxT i = 0; i < residual_y.size(); ++i) + { + DependencyTracking::Variable::DependencyMap dependency_y = (residual_y[i]).getDependencies(); + DependencyTracking::Variable::DependencyMap dependency_yp = (residual_yp[i]).getDependencies(); + + for (const auto& pair_y : dependency_y) + { + auto index_y = pair_y.first; + auto value_y = pair_y.second; + auto it_yp = dependency_yp.find(index_y); + if (it_yp != dependency_yp.end()) + { + auto value_yp = it_yp->second; + dependencies[i].insert(std::make_pair(index_y, value_y + value_yp)); + } + else + { + dependencies[i].insert(std::make_pair(index_y, value_y)); + } + } + + for (const auto& pair_yp : dependency_yp) + { + auto index_yp = pair_yp.first; + auto value_yp = pair_yp.second; + auto it_y = dependency_y.find(index_yp); + if (it_y == dependency_y.end()) + { + dependencies[i].insert(std::make_pair(index_yp, value_yp)); + } + } + } + + return dependencies; + } + + std::vector EnzymeJacobian() + { + ScalarT Vr1{1.0}; + ScalarT Vi1{0.0}; + PhasorDynamics::Bus bus(Vr1, Vi1); + auto data = makeGensalData(); + PhasorDynamics::Gensal gen(&bus, data); + + bus.allocate(); + gen.allocate(); + + bus.initialize(); + gen.initialize(); + + gen.updateTime(0.0, 1.0); + + for (size_t i = 0; i < bus.size(); ++i) + { + bus.setVariableIndex(i, i + gen.size()); + bus.setResidualIndex(i, i + gen.size()); + } + + bus.evaluateResidual(); + gen.evaluateResidual(); + + bus.evaluateJacobian(); + gen.evaluateJacobian(); + GridKit::LinearAlgebra::COO_Matrix& model_jacobian = gen.getJacobian(); + model_jacobian.deduplicate(); + + return GridKit::Testing::MapFromCOO(model_jacobian); + } +#endif + }; // class GensalTests + + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/runGensalTests.cpp b/tests/UnitTests/PhasorDynamics/runGensalTests.cpp new file mode 100644 index 000000000..c712c3e80 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/runGensalTests.cpp @@ -0,0 +1,18 @@ +#include "GensalTests.hpp" + +int main() +{ + GridKit::Testing::TestingResults result; + + GridKit::Testing::GensalTests test; + + result += test.constructor(); + result += test.hard_coded_residual(); + result += test.residual(); + result += test.residual_nonzero_ra(); + +#ifdef GRIDKIT_ENABLE_ENZYME + result += test.jacobian(); +#endif + return result.summary(); +} diff --git a/tests/UnitTests/Utilities/CaseFormatTests.hpp b/tests/UnitTests/Utilities/CaseFormatTests.hpp index d4da03e2d..b00ab2cf6 100644 --- a/tests/UnitTests/Utilities/CaseFormatTests.hpp +++ b/tests/UnitTests/Utilities/CaseFormatTests.hpp @@ -9,6 +9,7 @@ #include #include #include +#include #include #include #include @@ -63,6 +64,8 @@ namespace GridKit { "class": "Branch", "ports": {"bus1":1, "bus2":2}, "id": "1", "params": {"R":0.0, "X":0.1, "G":0.0, "B":0.0} }, { "class": "Genrou", "ports": {"bus":1}, "id": "1", "params": {"p0":1.0, "q0":0.05013, "H":3.0, "D":0.0, "Ra":0.0, "Tdop":7.0, "Tdopp":0.04, "Tqopp":0.05, "Tqop":0.75, "Xd":2.1, "Xdp":0.2, "Xdpp":0.18, "Xq":0.5, "Xqp": 0.0, "Xqpp":0.18, "Xl":0.15, "S10":0.0, "S12":0.0}, "mon": ["delta", "omega"] }, + { "class": "Gensal", "ports": {"bus":1}, "id": "2", "params": {"p0":1.0, "q0":0.05013, "H":3.0, "D":0.0, "Ra":0.0, "Tdop":7.0, "Tdopp":0.04, "Tqopp":0.05, + "Xd":2.1, "Xdp":0.2, "Xdpp":0.18, "Xq":0.5, "Xl":0.15, "S10":0.0, "S12":0.0}, "mon": ["delta", "omega"] }, { "class": "BusFault", "ports": {"bus":1}, "id": "1", "params": {"state0": false, "R":0.0, "X":1e-3} } ] })"; @@ -88,6 +91,7 @@ namespace GridKit success *= result.branch.size() == 1; success *= result.bus_fault.size() == 1; success *= result.genrou.size() == 1; + success *= result.gensal.size() == 1; success *= result.load.size() == 0; success *= result.bus[0].bus_id == 1; @@ -138,6 +142,15 @@ namespace GridKit success *= result.genrou[0].monitored_variables.contains(GenrouMonitorableVariables::delta); success *= result.genrou[0].monitored_variables.contains(GenrouMonitorableVariables::omega); + success *= std::get(result.gensal[0].parameters[GensalParameters::p0]) == 1.0; + success *= std::get(result.gensal[0].parameters[GensalParameters::q0]) == 0.05013; + success *= std::get(result.gensal[0].parameters[GensalParameters::Xd]) == 2.1; + success *= std::get(result.gensal[0].parameters[GensalParameters::Xq]) == 0.5; + success *= result.gensal[0].ports[GensalPorts::bus] == 1; + success *= result.gensal[0].disambiguation_string == "2"; + success *= result.gensal[0].monitored_variables.contains(GensalMonitorableVariables::delta); + success *= result.gensal[0].monitored_variables.contains(GensalMonitorableVariables::omega); + success *= std::get(result.bus_fault[0].parameters[BusFaultParameters::R]) == 0.0; success *= std::get(result.bus_fault[0].parameters[BusFaultParameters::X]) == 1e-3; success *= !std::get(result.bus_fault[0].parameters[BusFaultParameters::state0]);