diff --git a/CHANGELOG.md b/CHANGELOG.md index eb263a055..54c67fab0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,7 +20,7 @@ - Added clang formatting pre-commit. - CMake fixes. - Fixed warnings, memory leaks, and failed asserts. -- Added Tgov1 example. +- Added Tgov1 example. - Added 10 generator example. - Improved data structures. - Removed dead code. @@ -50,6 +50,7 @@ - Added phasor dynamics application to generalize examples - Added LoadZIP model component type. - Added component model developer checklist to a README file. +- Added IEEEST Stabilizer Model ## v0.1 diff --git a/GridKit/Model/PhasorDynamics/CMakeLists.txt b/GridKit/Model/PhasorDynamics/CMakeLists.txt index 025255776..b83192684 100644 --- a/GridKit/Model/PhasorDynamics/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/CMakeLists.txt @@ -37,6 +37,7 @@ add_subdirectory(Exciter) add_subdirectory(Governor) add_subdirectory(Load) add_subdirectory(LoadZIP) +add_subdirectory(Stabilizer) add_subdirectory(SynchronousMachine) add_subdirectory(SignalNode) diff --git a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp index 8b0c4281c..e74e1ac08 100644 --- a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp +++ b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp @@ -8,5 +8,6 @@ #include #include #include +#include #include #include diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp index c71adf733..507ee52b1 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp @@ -60,6 +60,7 @@ namespace GridKit OMEGA, ///< Generator speed deviation VREAL, ///< Real bus voltage VIMAG, ///< Imaginary bus voltage + VS, ///< Stabilizer output signal MAXIMUM, }; diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Data.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Data.hpp index 0d4553f82..99d0febc5 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Data.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Data.hpp @@ -40,6 +40,7 @@ namespace GridKit bus, ///< Unique ID of the terminal bus speed, ///< Unique ID of the generator speed signal efd, ///< Unique ID of the output efd signal + vs, ///< Unique ID of the stabilizer output signal (optional) }; /// Variables able to be monitored for a IEEET1 Exciter model diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp index 3284b292b..2cc266605 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp @@ -129,10 +129,12 @@ namespace GridKit wb_.resize(2); // Resize signal variable data - ws_.resize(1); - ws_indices_.resize(1); + ws_.resize(2); + ws_indices_.resize(2); ws_[0] = 0.0; ws_indices_[0] = INVALID_INDEX; + ws_[1] = 0.0; + ws_indices_[1] = INVALID_INDEX; // Default variable and residual index mapping to local index for (IdxT j = 0; j < size_; ++j) @@ -157,6 +159,7 @@ namespace GridKit int Ieeet1::verify() const { static constexpr auto OMEGA = Ieeet1ExternalVariables::OMEGA; + static constexpr auto VS = Ieeet1ExternalVariables::VS; int ret = 0; @@ -169,6 +172,15 @@ namespace GridKit } } + if (signals_.template isAttached()) + { + if (!signals_.template isLinked()) + { + Log::error() << "Ieeet1: VS signal attached with no linked source\n"; + ret += 1; + } + } + return ret; } @@ -202,24 +214,29 @@ namespace GridKit ScalarT vimag = bus_->Vi(); ScalarT Ec = std::sqrt(vreal * vreal + vimag * vimag); - // Derived from External initial values - ScalarT vr = Ke_ * efd0; + // Saturation at the initial operating point + ScalarT efd_sat = (efd0 - SA_) * Math::sigmoid(efd0 - SA_); + ScalarT ksat0 = SB_ * efd_sat * efd_sat; + ScalarT ve0 = ksat0 * efd0; + + // Derived from External initial values (includes saturation) + ScalarT vr = Ke_ * efd0 + ve0; ScalarT vfx = Kf_ / Tf_ * efd0; - ScalarT vtr = Ke_ / Ka_ * efd0; + ScalarT vtr = vr / Ka_; // Vref (setpoint = terminal + error) vref_ = Ec + vtr; // IVP for Internal Variables - y_[0] = Ec; // y0 - vts - Sensed term volt - y_[1] = vr; // y1 - vr - Voltage reg - y_[2] = efd0; // y2 - efdp - Efd pre mult - y_[3] = vfx; // y3 - vfx - Exciter feedback - y_[4] = vtr; // y4 - vtr - Term Volt Err - y_[5] = 0; // y5 - vf - Feedback volt - y_[6] = 0; // y6 - ve - Excit. Cntrl Volt - y_[7] = efd0; // y7 - efd - Efd - y_[8] = 0; // y8 - ksat - Saturation + y_[0] = Ec; // y0 - vts - Sensed term volt + y_[1] = vr; // y1 - vr - Voltage reg + y_[2] = efd0; // y2 - efdp - Efd pre mult + y_[3] = vfx; // y3 - vfx - Exciter feedback + y_[4] = vtr; // y4 - vtr - Term Volt Err + y_[5] = 0; // y5 - vf - Feedback volt + y_[6] = ve0; // y6 - ve - Excit. Cntrl Volt + y_[7] = efd0; // y7 - efd - Efd + y_[8] = ksat0; // y8 - ksat - Saturation // Steady State Conditions yp_[0] = 0.0; @@ -294,7 +311,8 @@ namespace GridKit ScalarT vfx_dot = yp[3]; // Set signal variable aliases - ScalarT omega = ws[0]; + ScalarT omega = ws[0]; + ScalarT vs_signal = ws[1]; // The 'pre-limit' derivative of Pv ScalarT func = -vr + Ka_ * vtr; @@ -308,7 +326,7 @@ namespace GridKit f[3] = -vfx_dot + vf / Tf_; // Internal Algebraic Equations - f[4] = -vts + vref_ + vUEL_ + vOEL_ + vS_ - vtr - vf; + f[4] = -vts + vref_ + vUEL_ + vOEL_ + vs_signal - vtr - vf; f[5] = -vf + (efdp * Kf_) / Tf_ - vfx; f[6] = -ve + ksat * efdp; f[7] = -efd + efdp + omega * efdp * Ispdlim_; @@ -341,6 +359,15 @@ namespace GridKit ws_indices_[0] = signals_.template readExternalVariableIndex(); } + // VS signal (stabilizer output, optional) + ws_[1] = 0.0; + ws_indices_[1] = INVALID_INDEX; + if (signals_.template isAttached()) + { + ws_[1] = signals_.template readExternalVariable(); + ws_indices_[1] = signals_.template readExternalVariableIndex(); + } + // Bus voltages wb_[0] = bus_->Vr(); wb_[1] = bus_->Vi(); diff --git a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md index e033f1e7a..bccf383c1 100644 --- a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md +++ b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md @@ -85,7 +85,7 @@ a bus and has the following fields: `class` | A string indicating the class of node. See the table below for more information `name` | Optional string containing the name of the node. This may be empty or non-unique `init` | Optional object mapping string variable names to floating point values, specifying default voltages or signal values. The available initialization variables are dependent upon the node class. Any variables missing will be given default values, which are specified beneath the table below. If this object is missing, all variables will be given default values. See the table below for more information - `v_base` | Optional floating point value giving the voltage base in volts (V). + `v_base` | Optional floating point value giving the voltage base in volts (V). `mon` | Optional field, which is an array specifying variables to monitor the value of in an output channel. Available variables include all the initialization variables, along with others as determined by the node class. See the table below for more information `freq_base` | Optional field to override the system frequency base at this bus `va_base` | Optional field to override the system power base at this bus @@ -146,7 +146,8 @@ are specified: `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` `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` | `Tr`, `Ka`, `Ta`, `Ke`, `Te`, `Kf`, `Tf`, `Vrmin`, `Vrmax`, `E1`, `E2`, `Se1`, `Se2`, `Ispdlim` | `efd`, `ksat` + `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` + `Ieeest` | the IEEEST stabilizer model | `input`, `output`, `cutout`\* | `A1`, `A2`, `A3`, `A4`, `A5`, `A6`, `T1`, `T2`, `T3`, `T4`, `T5`, `T6`, `Ks`, `Lsmin`, `Lsmax`, `Vcl`, `Vcu`, `Tdelay` | `vs` `BusFault` | simple impedance-based fault at a bus | `bus`, `status`\* | `state0`, `R`, `X` | `state`, `ir`, `ii` Ports marked with \* are optional and, if missing, will be assumed to be diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Stabilizer/CMakeLists.txt new file mode 100644 index 000000000..e0d818124 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Stabilizer/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory(IEEEST) diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/CMakeLists.txt new file mode 100644 index 000000000..6fb7da2fa --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/CMakeLists.txt @@ -0,0 +1,50 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +set(_install_headers + Ieeest.hpp + IeeestData.hpp) + +if(GRIDKIT_ENABLE_ENZYME) + gridkit_add_library(phasor_dynamics_stabilizer_ieeest + SOURCES + IeeestEnzyme.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) +else() + gridkit_add_library(phasor_dynamics_stabilizer_ieeest + SOURCES + Ieeest.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_stabilizer_ieeest_dependency_tracking + SOURCES + IeeestDependencyTracking.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_stabilizer_ieeest) +target_link_libraries(phasor_dynamics_components_dependency_tracking + INTERFACE GridKit::phasor_dynamics_stabilizer_ieeest_dependency_tracking) diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.cpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.cpp new file mode 100644 index 000000000..dad12b933 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.cpp @@ -0,0 +1,35 @@ +/** + * @file Ieeest.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Jacobian stub and template instantiations for IEEEST Stabilizer. + */ + +#include "IeeestImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Stabilizer + { + /** + * @brief Jacobian evaluation not implemented yet + * + * @tparam ScalarT - Scalar data type + * @tparam IdxT - Index data type + * @return int - error code, 0 = success + */ + template + int Ieeest::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Ieeest..." << std::endl; + Log::misc() << "Jacobian evaluation not implemented!" << std::endl; + return 0; + } + + // Available template instantiations + template class Ieeest; + template class Ieeest; + } // namespace Stabilizer + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp new file mode 100644 index 000000000..7f57322aa --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp @@ -0,0 +1,178 @@ +/** + * @file Ieeest.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Declaration of the IEEEST Power System Stabilizer. + */ + +#pragma once + +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Stabilizer + { + template + struct IeeestData; + } // namespace Stabilizer + + template + class SignalNode; + + } // namespace PhasorDynamics +} // namespace GridKit + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Stabilizer + { + /// Internal variables of a `Ieeest` + enum class IeeestInternalVariables : size_t + { + X1, ///< Notch filter state 1 + X2, ///< Notch filter state 2 + X3, ///< Notch filter state 3 + X4, ///< Notch filter state 4 + X5, ///< Lead-lag 1 state + X6, ///< Lead-lag 2 state + X7, ///< Washout state + V4, ///< Notch filter output + V5, ///< Lead-lag 1 output + V6, ///< Lead-lag 2 output + V7, ///< Unlimited stabilizer signal + VSS, ///< Limited stabilizer signal + VS, ///< Stabilizer output + MAXIMUM, + }; + + /// External variables of a `Ieeest` + enum class IeeestExternalVariables : size_t + { + U, ///< Stabilizer input signal + VCT, ///< Cutout signal + MAXIMUM, + }; + + template + class Ieeest : 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::variable_indices_; + using Component::residual_indices_; + + public: + using RealT = typename Component::RealT; + using model_data_type = IeeestData; + using signal_type = SignalNode; + using MonitorT = Model::VariableMonitor; + + Ieeest(); + Ieeest(const model_data_type& data); + ~Ieeest(); + + 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; + int evaluateJacobian() override final; + + /// Get the `ComponentSignals` from this `Ieeest` + auto getSignals() + -> ComponentSignals& + { + return signals_; + } + + const Model::VariableMonitorBase* getMonitor() const override; + + __attribute__((always_inline)) inline int evaluateInternalResidual( + ScalarT*, + ScalarT*, + [[maybe_unused]] ScalarT*, + ScalarT*, + ScalarT*); + + private: + RealT A1_{0}; + RealT A2_{0}; + RealT A3_{0}; + RealT A4_{0}; + RealT A5_{0}; + RealT A6_{0}; + RealT T1_{0}; + RealT T2_{1}; + RealT T3_{0}; + RealT T4_{1}; + RealT T5_{0}; + RealT T6_{1}; + RealT Ks_{1}; + RealT Lsmin_{-0.1}; + RealT Lsmax_{0.1}; + RealT Vcl_{0}; + RealT Vcu_{1.5}; + RealT Tdelay_{0}; + + RealT a0_{1}; + RealT a1_{0}; + RealT a2_{0}; + RealT a3_{0}; + RealT a4_{0}; + + // Precomputed masks and safe inverse coefficients for branch-free degenerate paths. + RealT use_4th_order_{0}; + RealT use_3rd_order_{0}; + RealT use_2nd_order_{0}; + RealT safe_inv_a4_{0}; + RealT safe_inv_a3_{0}; + RealT safe_inv_a2_{0}; + RealT use_T2_block_{1}; + RealT bypass_T2_block_{0}; + RealT safe_inv_T2_{1}; + RealT use_T4_block_{1}; + RealT bypass_T4_block_{0}; + RealT safe_inv_T4_{1}; + RealT use_T6_block_{1}; + RealT bypass_T6_block_{0}; + RealT safe_inv_T6_{1}; + RealT use_cutout_{1}; + RealT bypass_cutout_{0}; + + ComponentSignals signals_; + + std::unique_ptr monitor_; + + void initializeParameters(const model_data_type& data); + void initializeMonitor(); + + std::vector ws_; + std::vector ws_indices_; + }; + + } // namespace Stabilizer + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestData.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestData.hpp new file mode 100644 index 000000000..00d728a58 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestData.hpp @@ -0,0 +1,82 @@ +/** + * @file IeeestData.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Modeling data for IEEEST Stabilizer + */ + +#pragma once + +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Stabilizer + { + /** + * @brief Parameter keys for IEEEST Stabilizer model. + */ + enum class IeeestParameters + { + A1, ///< Notch filter denominator coefficient + A2, ///< Notch filter denominator coefficient + A3, ///< Notch filter denominator coefficient + A4, ///< Notch filter denominator coefficient + A5, ///< Notch filter numerator coefficient + A6, ///< Notch filter numerator coefficient + T1, ///< Lead-lag 1 numerator time constant + T2, ///< Lead-lag 1 denominator time constant + T3, ///< Lead-lag 2 numerator time constant + T4, ///< Lead-lag 2 denominator time constant + T5, ///< Washout numerator time constant + T6, ///< Washout denominator time constant + Ks, ///< Stabilizer gain + Lsmin, ///< Minimum stabilizer output limit + Lsmax, ///< Maximum stabilizer output limit + Vcl, ///< Lower input cutout threshold + Vcu, ///< Upper input cutout threshold + Tdelay, ///< Input time delay (not modeled) + }; + + /** + * @brief Port keys for IEEEST Stabilizer model. + */ + enum class IeeestPorts + { + input, ///< Unique ID of the stabilizer input signal + cutout, ///< Unique ID of the cutout signal + output, ///< Unique ID of the stabilizer output signal + }; + + /** + * @brief Monitorable variables for IEEEST Stabilizer model. + */ + enum class IeeestMonitorableVariables + { + vs, ///< Stabilizer output + }; + + /** + * @brief Contains modeling data for a IEEEST Stabilizer model. + * + * @tparam RealT Real parameter data type + * @tparam IdxT Integer parameter data type + */ + template + struct IeeestData : public ComponentData + { + IeeestData() = default; + + using Parameters = IeeestParameters; + using Ports = IeeestPorts; + using MonitorableVariables = IeeestMonitorableVariables; + }; + + } // namespace Stabilizer + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestDependencyTracking.cpp new file mode 100644 index 000000000..b61781f56 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestDependencyTracking.cpp @@ -0,0 +1,35 @@ +/** + * @file IeeestDependencyTracking.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief DependencyTracking Jacobian stub and template instantiations for IEEEST Stabilizer. + */ + +#include "IeeestImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Stabilizer + { + /** + * @brief Jacobian evaluation not implemented yet + * + * @tparam ScalarT - Scalar data type + * @tparam IdxT - Index data type + * @return int - error code, 0 = success + */ + template + int Ieeest::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Ieeest..." << std::endl; + Log::misc() << "Jacobian evaluation not implemented!" << std::endl; + return 0; + } + + // Available template instantiations + template class Ieeest; + template class Ieeest; + } // namespace Stabilizer + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp new file mode 100644 index 000000000..71c8a086e --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp @@ -0,0 +1,84 @@ +/** + * @file IeeestEnzyme.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Enzyme-based sparse Jacobian for IEEEST Stabilizer. + */ + +#include + +#include "IeeestImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Stabilizer + { + /** + * @brief Jacobian evaluation experimental + * + * @tparam ScalarT - Scalar data type + * @tparam IdxT - Index data type + * @return int - error code, 0 = success + */ + template + int Ieeest::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Ieeest..." << 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::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_); + + return 0; + } + + // Available template instantiations + template class Ieeest; + template class Ieeest; + + } // namespace Stabilizer + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp new file mode 100644 index 000000000..1aba2b73c --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp @@ -0,0 +1,358 @@ +#pragma once + +/** + * @file IeeestImpl.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Definition of the IEEEST Power System Stabilizer. + */ + +#include +#include + +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Stabilizer + { + using Log = ::GridKit::Utilities::Logger; + + template + Ieeest::Ieeest() + { + size_ = 13; + } + + template + Ieeest::Ieeest(const model_data_type& data) + : monitor_(std::make_unique(data)) + { + initializeParameters(data); + initializeMonitor(); + size_ = 13; + } + + template + Ieeest::~Ieeest() + { + } + + template + void Ieeest::initializeParameters(const model_data_type& data) + { + if (data.parameters.contains(model_data_type::Parameters::A1)) + { + A1_ = std::get(data.parameters.at(model_data_type::Parameters::A1)); + } + if (data.parameters.contains(model_data_type::Parameters::A2)) + { + A2_ = std::get(data.parameters.at(model_data_type::Parameters::A2)); + } + if (data.parameters.contains(model_data_type::Parameters::A3)) + { + A3_ = std::get(data.parameters.at(model_data_type::Parameters::A3)); + } + if (data.parameters.contains(model_data_type::Parameters::A4)) + { + A4_ = std::get(data.parameters.at(model_data_type::Parameters::A4)); + } + if (data.parameters.contains(model_data_type::Parameters::A5)) + { + A5_ = std::get(data.parameters.at(model_data_type::Parameters::A5)); + } + if (data.parameters.contains(model_data_type::Parameters::A6)) + { + A6_ = std::get(data.parameters.at(model_data_type::Parameters::A6)); + } + if (data.parameters.contains(model_data_type::Parameters::T1)) + { + T1_ = std::get(data.parameters.at(model_data_type::Parameters::T1)); + } + if (data.parameters.contains(model_data_type::Parameters::T2)) + { + T2_ = std::get(data.parameters.at(model_data_type::Parameters::T2)); + } + if (data.parameters.contains(model_data_type::Parameters::T3)) + { + T3_ = std::get(data.parameters.at(model_data_type::Parameters::T3)); + } + if (data.parameters.contains(model_data_type::Parameters::T4)) + { + T4_ = std::get(data.parameters.at(model_data_type::Parameters::T4)); + } + if (data.parameters.contains(model_data_type::Parameters::T5)) + { + T5_ = std::get(data.parameters.at(model_data_type::Parameters::T5)); + } + if (data.parameters.contains(model_data_type::Parameters::T6)) + { + T6_ = std::get(data.parameters.at(model_data_type::Parameters::T6)); + } + if (data.parameters.contains(model_data_type::Parameters::Ks)) + { + Ks_ = std::get(data.parameters.at(model_data_type::Parameters::Ks)); + } + if (data.parameters.contains(model_data_type::Parameters::Lsmin)) + { + Lsmin_ = std::get(data.parameters.at(model_data_type::Parameters::Lsmin)); + } + if (data.parameters.contains(model_data_type::Parameters::Lsmax)) + { + Lsmax_ = std::get(data.parameters.at(model_data_type::Parameters::Lsmax)); + } + if (data.parameters.contains(model_data_type::Parameters::Vcl)) + { + Vcl_ = std::get(data.parameters.at(model_data_type::Parameters::Vcl)); + } + if (data.parameters.contains(model_data_type::Parameters::Vcu)) + { + Vcu_ = std::get(data.parameters.at(model_data_type::Parameters::Vcu)); + } + if (data.parameters.contains(model_data_type::Parameters::Tdelay)) + { + Tdelay_ = std::get(data.parameters.at(model_data_type::Parameters::Tdelay)); + } + + a0_ = 1; + a1_ = A1_ + A3_; + a2_ = A2_ + A4_ + A1_ * A3_; + a3_ = A1_ * A4_ + A2_ * A3_; + a4_ = A2_ * A4_; + + // Precompute masks and safe inverse coefficients so the residual stays branch-free. + use_4th_order_ = static_cast(a4_ != 0.0); + use_3rd_order_ = static_cast(a4_ == 0.0 && a3_ != 0.0); + use_2nd_order_ = static_cast(a4_ == 0.0 && a3_ == 0.0 && a2_ != 0.0); + safe_inv_a4_ = use_4th_order_ / (a4_ + (1.0 - use_4th_order_)); + safe_inv_a3_ = use_3rd_order_ / (a3_ + (1.0 - use_3rd_order_)); + safe_inv_a2_ = use_2nd_order_ / (a2_ + (1.0 - use_2nd_order_)); + + use_T2_block_ = static_cast(T2_ != 0.0); + bypass_T2_block_ = 1.0 - use_T2_block_; + safe_inv_T2_ = use_T2_block_ / (T2_ + bypass_T2_block_); + + use_T4_block_ = static_cast(T4_ != 0.0); + bypass_T4_block_ = 1.0 - use_T4_block_; + safe_inv_T4_ = use_T4_block_ / (T4_ + bypass_T4_block_); + + use_T6_block_ = static_cast(T6_ != 0.0); + bypass_T6_block_ = 1.0 - use_T6_block_; + safe_inv_T6_ = use_T6_block_ / (T6_ + bypass_T6_block_); + + // Vcl = Vcu = 0 means "cutout disabled" (passthrough: vs = vss). + use_cutout_ = static_cast(Vcl_ != 0.0 || Vcu_ != 0.0); + bypass_cutout_ = 1.0 - use_cutout_; + } + + template + int Ieeest::setGridKitComponentID(IdxT component_id) + { + gridkit_component_id_ = component_id; + return 0; + } + + template + int Ieeest::allocate() + { + 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); + + ws_.resize(2); + ws_indices_.resize(2); + ws_[0] = 0.0; + ws_indices_[0] = INVALID_INDEX; + ws_[1] = 0.0; + ws_indices_[1] = INVALID_INDEX; + + for (IdxT j = 0; j < size_; ++j) + { + this->setVariableIndex(j, j); + this->setResidualIndex(j, j); + } + + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y_[12], &(this->getVariableIndex(12))); + } + + return 0; + } + + template + int Ieeest::verify() const + { + int ret = 0; + + if (signals_.template isAttached()) + { + if (!signals_.template isLinked()) + { + Log::error() << "Ieeest: input signal U attached with no linked source\n"; + ret += 1; + } + } + else + { + Log::error() << "Ieeest: required input signal U is not attached\n"; + ret += 1; + } + + if (signals_.template isAttached()) + { + if (!signals_.template isLinked()) + { + Log::error() << "Ieeest: cutout signal VCT attached with no linked source\n"; + ret += 1; + } + } + + if (a4_ == 0 && a3_ == 0 && a2_ == 0) + { + Log::error() << "Ieeest: a2, a3, and a4 are all zero - no valid notch filter\n"; + ret += 1; + } + + return ret; + } + + template + int Ieeest::initialize() + { + for (IdxT i = 0; i < size_; ++i) + { + y_[static_cast(i)] = 0.0; + yp_[static_cast(i)] = 0.0; + } + + return 0; + } + + template + int Ieeest::tagDifferentiable() + { + tag_[0] = true; + tag_[1] = true; + tag_[2] = true; + tag_[3] = true; + tag_[4] = true; + tag_[5] = true; + tag_[6] = true; + tag_[7] = false; + tag_[8] = false; + tag_[9] = false; + tag_[10] = false; + tag_[11] = false; + tag_[12] = false; + + return 0; + } + + template + __attribute__((always_inline)) inline int Ieeest::evaluateInternalResidual( + ScalarT* y, + ScalarT* yp, + [[maybe_unused]] ScalarT* wb, + ScalarT* ws, + ScalarT* f) + { + ScalarT x1 = y[0]; + ScalarT x2 = y[1]; + ScalarT x3 = y[2]; + ScalarT x4 = y[3]; + ScalarT x5 = y[4]; + ScalarT x6 = y[5]; + ScalarT x7 = y[6]; + ScalarT v4 = y[7]; + ScalarT v5 = y[8]; + ScalarT v6 = y[9]; + ScalarT v7 = y[10]; + ScalarT vss = y[11]; + ScalarT vs = y[12]; + + ScalarT x1_dot = yp[0]; + ScalarT x2_dot = yp[1]; + ScalarT x3_dot = yp[2]; + ScalarT x4_dot = yp[3]; + ScalarT x5_dot = yp[4]; + ScalarT x6_dot = yp[5]; + ScalarT x7_dot = yp[6]; + + ScalarT u = ws[0]; + ScalarT vct = ws[1]; + + f[0] = -x1_dot + x2; + f[1] = -x2_dot + (use_4th_order_ + use_3rd_order_) * x3 + + use_2nd_order_ * (-a0_ * x1 - a1_ * x2 + u) * safe_inv_a2_; + f[2] = -x3_dot + use_4th_order_ * x4 + + use_3rd_order_ * (-a0_ * x1 - a1_ * x2 - a2_ * x3 + u) * safe_inv_a3_; + f[3] = -x4_dot + use_4th_order_ * (-a0_ * x1 - a1_ * x2 - a2_ * x3 - a3_ * x4 + u) * safe_inv_a4_; + + f[4] = -x5_dot + use_T2_block_ * (v4 - x5) * safe_inv_T2_; + f[5] = -x6_dot + use_T4_block_ * (v5 - x6) * safe_inv_T4_; + f[6] = -x7_dot + use_T6_block_ * (v6 - x7) * safe_inv_T6_; + + f[7] = -v4 + x1 + A5_ * x2 + (use_4th_order_ + use_3rd_order_) * A6_ * x3; + f[8] = -v5 + bypass_T2_block_ * v4 + use_T2_block_ * (x5 + T1_ * (v4 - x5) * safe_inv_T2_); + f[9] = -v6 + bypass_T4_block_ * v5 + use_T4_block_ * (x6 + T3_ * (v5 - x6) * safe_inv_T4_); + f[10] = -v7 + bypass_T6_block_ * Ks_ * v6 + use_T6_block_ * Ks_ * T5_ * (v6 - x7) * safe_inv_T6_; + + f[11] = -vss + v7 * Math::indicator_zero(Lsmin_, Lsmax_, v7) + + Lsmin_ * Math::sigmoid(Lsmin_ - v7) + + Lsmax_ * Math::sigmoid(v7 - Lsmax_); + // Cutout: Vcl=Vcu=0 means disabled (passthrough). Otherwise, sigmoid window. + f[12] = -vs + bypass_cutout_ * vss + + use_cutout_ * vss * Math::sigmoid(vct - Vcl_) * Math::sigmoid(Vcu_ - vct); + + return 0; + } + + template + int Ieeest::evaluateResidual() + { + if (signals_.template isAttached()) + { + ws_[0] = signals_.template readExternalVariable(); + ws_indices_[0] = signals_.template readExternalVariableIndex(); + } + + ws_[1] = (Vcl_ + Vcu_) / TWO; + if (signals_.template isAttached()) + { + ws_[1] = signals_.template readExternalVariable(); + ws_indices_[1] = signals_.template readExternalVariableIndex(); + } + + evaluateInternalResidual(y_.data(), yp_.data(), wb_.data(), ws_.data(), f_.data()); + + return 0; + } + + template + const Model::VariableMonitorBase* Ieeest::getMonitor() const + { + return monitor_.get(); + } + + template + void Ieeest::initializeMonitor() + { + using Variable = typename model_data_type::MonitorableVariables; + monitor_->set(Variable::vs, [this] + { return y_[12]; }); + } + + } // namespace Stabilizer + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md new file mode 100644 index 000000000..7a5e7b57c --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md @@ -0,0 +1,125 @@ +# IEEEST + +The **IEEEST** model is a standard IEEE power system stabilizer used in transient stability simulations. +It consists of a 4th-order notch filter, two lead–lag blocks, a washout block, and an output limiter with input cutout logic. + +Notes: +- The **cutout logic uses** $V_{ct}$ (as labeled in the block diagram), not $u_d$. + +## Block Diagram + +
+ + + Figure 1: Stabilizer IEEEST model. Figure courtesy of [PowerWorld](https://www.powerworld.com/WebHelp/) +
+ + +## Model Parameters + +Symbol | Units | Description | Typical Value | Note +------ | ----- | ----------- | -------------- | ---- +$A_1$ | [s] | Notch filter denominator coefficient | 1.013 +$A_2$ | [s²] | Notch filter denominator coefficient | 0.013 +$A_3$ | [s] | Notch filter denominator coefficient | 0.0 +$A_4$ | [s²] | Notch filter denominator coefficient | 0.0 +$A_5$ | [s] | Notch filter numerator coefficient | 1.013 +$A_6$ | [s²] | Notch filter numerator coefficient | 0.113 +$T_1$ | [s] | Lead–lag 1 numerator time constant | 0.0 +$T_2$ | [s] | Lead–lag 1 denominator time constant | 0.02 +$T_3$ | [s] | Lead–lag 2 numerator time constant | 0.0 +$T_4$ | [s] | Lead–lag 2 denominator time constant | 0.0 +$T_5$ | [s] | Washout numerator time constant | 1.65 +$T_6$ | [s] | Washout denominator time constant | 1.65 +$K_s$ | [p.u.] | Stabilizer gain | 3.0 +$L_{s\min}$ | [p.u.] | Minimum stabilizer output limit | 0.1 +$L_{s\max}$ | [p.u.] | Maximum stabilizer output limit | -0.1 +$V_{cl}$ | [p.u.] | Lower input cutout threshold | 0.0 +$V_{cu}$ | [p.u.] | Upper input cutout threshold | 0.0 +$T_{delay}$ | [s] | Input time delay | 0.0 + +### Model Derived Parameters +```math +\begin{aligned} +a_0 &= 1 \\ +a_1 &= A_1 + A_3 \\ +a_2 &= A_2 + A_4 + A_1 A_3 \\ +a_3 &= A_1 A_4 + A_2 A_3 \\ +a_4 &= A_2 A_4 +\end{aligned} +``` + +## Model Variables + +### Internal Variables + +#### Differential + +Symbol | Units | Description | Note +------ | ----- | ----------- | ---- +$x_1$ | [-] | Notch filter state | +$x_2$ | [-] | Notch filter state | +$x_3$ | [-] | Notch filter state | +$x_4$ | [-] | Notch filter state | +$x_5$ | [-] | Lead–lag 1 state | +$x_6$ | [-] | Lead–lag 2 state | +$x_7$ | [-] | Washout state | + +#### Algebraic + +Symbol | Units | Description | Note +------ | ----- | ----------- | ---- +$v_4$ | [p.u.] | Notch filter output | +$v_5$ | [p.u.] | Lead–lag 1 output | +$v_6$ | [p.u.] | Lead–lag 2 output | +$v_7$ | [p.u.] | Unlimited stabilizer signal | +$V_{ss}$ | [p.u.] | Limited stabilizer signal | +$V_s$ | [p.u.] | Stabilizer output | + +### External Variables + +#### Differential + +None. + +#### Algebraic + +Symbol | Units | Description | Note +------ | ----- | ----------- | ---- +$u$ | [p.u.] | Stabilizer input signal | +$V_{ct}$ | [p.u.] | Cutout signal (compared to $V_{cl},V_{cu}$) | from the block diagram + +## Model Equations + +### Differential Equations +```math +\begin{aligned} +\dot{x}_1 &= x_2 \\ +\dot{x}_2 &= x_3 \\ +\dot{x}_3 &= x_4 \\ +\dot{x}_4 &= -\dfrac{a_0}{a_4}x_1 + -\dfrac{a_1}{a_4}x_2 + -\dfrac{a_2}{a_4}x_3 + -\dfrac{a_3}{a_4}x_4 + +\dfrac{1}{a_4}u\\ +\dot{x}_5 &= \dfrac{1}{T_2}(v_4 - x_5) \\ +\dot{x}_6 &= \dfrac{1}{T_4}(v_5 - x_6)\\ +\dot{x}_7 &= \dfrac{1}{T_6}(v_6 - x_7) +\end{aligned} +``` + +### Algebraic Equations +```math +\begin{aligned} +0 &= -v_4 + x_1 + A_5 x_2 + A_6 x_3 \\ +0 &= -v_5 + x_5 + \dfrac{T_1}{T_2}(v_4 - x_5) \\ +0 &= -v_6 + x_6 + \dfrac{T_3}{T_4}(v_5 - x_6) \\ +0 &= -v_7 + K_s \dfrac{T_5}{T_6}(v_6 - x_7) \\ +0 &= -V_{ss} + \min\!\big(\max(v_7, L_{s\min}), L_{s\max}\big) \\ +0 &= -V_s + +\begin{cases} +V_{ss}, & V_{cl} < V_{ct} < V_{cu} \\ +0, & \text{otherwise} +\end{cases} +\end{aligned} +``` diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/PSS1A/README.md b/GridKit/Model/PhasorDynamics/Stabilizer/PSS1A/README.md new file mode 100644 index 000000000..f380cc21c --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Stabilizer/PSS1A/README.md @@ -0,0 +1,102 @@ +# PSS1A + +> [!NOTE] +> This is not yet implemented + +> [!NOTE] +> The Parameters, variables, and equations need to be formatted and verified - this is WIP. + + + +
+ + + + Figure 1: Power system stabilizer PSS1A model. Figure courtesy of [PowerWorld](https://www.powerworld.com/WebHelp/) +
+ + + +## Model Parameters + + +- $I_{cs}$ - stabilizer input code, (2) +- $A_{1}$ - notch filter parameters, (0) +- $A_{2}$ - notch filter parameters, (0) +- $T_{1}$ - lead/lag time constant, sec (0.25) +- $T_{2}$ - lead/lag time constant, sec (0.03) +- $T_{3}$ - lead/lag time constant, sec (0.25) +- $T_{4}$ - lead/lag time constant, sec (0.03) +- $T_{5}$ - washout numerator time constant, sec (20) +- $T_{6}$ - washout denomirator time constant/transducer time constant, sec (0.02) +- $K_{S}$ - stabilizer gains, (10) +- $L_{smax}$ - maximum stabilizer output, pu (0.1) +- $L_{smin}$ - minimum stabilizer output, pu (-0.1) +- $V_{cu}$ - stabilizer input cutoff threshold, pu (0) +- $V_{cl}$ - stabilizer input cutoff threshold, pu (0) + + + + +## Model Variables + +These were the variables listed in the old documentation. + +1. rotor speed deviation (p.u.) +2. bus frequency deviation (p.u.) - default +3. generator electrical power in Gen MVA Base (p.u.) +4. generator accelerating power (p.u.) +5. bus voltage (p.u.) +6. derivative of p.u. bus voltage + +### Internal Variables + +#### Differential + +TBD + +#### Algebraic + +TBD + +### External Variables + +#### Differential + +None. + +#### Algebraic + +Symbol | Units | Description | Note +------ | ----- | ----------- | ---- +$u$ | [p.u.] | Stabilizer input signal | +$V_{ct}$ | [p.u.] | Cutout signal (compared to $V_{cl},V_{cu}$) | from the block diagram + + +### Differential Equations + +```math +\begin{aligned} +\dot{V_{1}} &= \dfrac{1}{ T_{6} }( V_{SI}-V_{1} ) \\ +\dot{x_{1}} &= -\dfrac{ V_{2} }{ T_{5} } \\ +\dfrac{d^{2}V_{3}}{dt^{2}}+\dfrac{A_{1}}{A_{2}}\dfrac{dV_{3}}{dt}&=\dfrac{1}{A_{2}}(V_{2}-1) \\ +\dfrac{dx_{2}}{dt}&=\dfrac{1}{T_{2}}(V_{3}-V_{4}) \\ +\dfrac{dx_{3}}{dt}&=\dfrac{1}{T_{4}}(V_{4}-V_{5}) + +\end{aligned} +``` + +### Algebraic Equations + +```math +\begin{aligned} +V_{2} &= x_{1} + K_{S} V_{1} \\ +V_{4}&=x_{2}+\dfrac{T_{1}}{T_{2}}V_{3} \\ +V_{5}&=x_{3}+\dfrac{T_{3}}{T_{4}}V_{4} \\ +V_{llout} &= \begin{cases} + L_{SMAX} &\text{if } V_{5}>V_{SMAX} \\ + L_{SMIN} &\text{if } V_{5} - - - - Figure 1: Power system stabilizer PSS1A model. Figure courtesy of [PowerWorld](https://www.powerworld.com/WebHelp/) - +## Types -## Nomenclature +The GridKit implemented stabilizers include: -### Inputs ($`I_{cs}`$) - -1. rotor speed deviation (p.u.) -2. bus frequency deviation (p.u.) - default -3. generator electrical power in Gen MVA Base (p.u.) -4. generator accelerating power (p.u.) -5. bus voltage (p.u.) -6. derivative of p.u. bus voltage - -### Parameters -- $`I_{cs}`$ - stabilizer input code, (2) -- $`A_{1}`$ - notch filter parameters, (0) -- $`A_{2}`$ - notch filter parameters, (0) -- $`T_{1}`$ - lead/lag time constant, sec (0.25) -- $`T_{2}`$ - lead/lag time constant, sec (0.03) -- $`T_{3}`$ - lead/lag time constant, sec (0.25) -- $`T_{4}`$ - lead/lag time constant, sec (0.03) -- $`T_{5}`$ - washout numerator time constant, sec (20) -- $`T_{6}`$ - washout denomirator time constant/transducer time constant, sec (0.02) -- $`K_{S}`$ - stabilizer gains, (10) -- $`L_{smax}`$ - maximum stabilizer output, pu (0.1) -- $`L_{smin}`$ - minimum stabilizer output, pu (-0.1) -- $`V_{cu}`$ - stabilizer input cutoff threshold, pu (0) -- $`V_{cl}`$ - stabilizer input cutoff threshold, pu (0) - - -## Equations -First block -```math -\dfrac{dV_{1}}{dt}=\dfrac{1}{T_{6}}(V_{SI}-V_{1}) -``` -Second block -```math -\dfrac{dx_{1}}{dt}=-\dfrac{V_{2}}{T_{5}} -``` -```math -V_{2}=x_{1}+K_{S}V_{1} -``` -Third block -```math -\dfrac{d^{2}V_{3}}{dt^{2}}+\dfrac{A_{1}}{A_{2}}\dfrac{dV_{3}}{dt}=\dfrac{1}{A_{2}}(V_{2}-1) -``` -Fourth block -```math -\dfrac{dx_{2}}{dt}=\dfrac{1}{T_{2}}(V_{3}-V_{4}) -``` -```math -V_{4}=x_{2}+\dfrac{T_{1}}{T_{2}}V_{3} -``` -Fifth block -```math -\dfrac{dx_{3}}{dt}=\dfrac{1}{T_{4}}(V_{4}-V_{5}) -``` -```math -V_{5}=x_{3}+\dfrac{T_{3}}{T_{4}}V_{4} -``` -```math -V_{llout} = \begin{cases} - L_{SMAX} &\text{if } V_{5}>V_{SMAX} \\ - L_{SMIN} &\text{if } V_{5}getSignals().template assignSignalNode(getSignal(efd)); } + if (excitedata.ports.contains(Ieeet1Data::Ports::vs)) + { + IdxT vs = excitedata.ports.at(Ieeet1Data::Ports::vs); + exciter->getSignals().template attachSignalNode(getSignal(vs)); + } + addComponent(exciter); } + // Add IEEEST stabilizers + for (const auto& stabdata : data.stabilizer) + { + auto* stabilizer = new Ieeest(stabdata); + + if (stabdata.ports.contains(IeeestPorts::input)) + { + IdxT input = stabdata.ports.at(IeeestPorts::input); + stabilizer->getSignals().template attachSignalNode(getSignal(input)); + } + + if (stabdata.ports.contains(IeeestPorts::cutout)) + { + IdxT cutout = stabdata.ports.at(IeeestPorts::cutout); + stabilizer->getSignals().template attachSignalNode(getSignal(cutout)); + } + + if (stabdata.ports.contains(IeeestPorts::output)) + { + IdxT output = stabdata.ports.at(IeeestPorts::output); + stabilizer->getSignals().template assignSignalNode(getSignal(output)); + } + + addComponent(stabilizer); + } + // Add faults for (const auto& faultdata : data.bus_fault) { diff --git a/GridKit/Model/PhasorDynamics/SystemModelData.hpp b/GridKit/Model/PhasorDynamics/SystemModelData.hpp index dabfab8c6..ff3b53e57 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelData.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelData.hpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -35,6 +36,7 @@ namespace GridKit using BusFaultDataT = BusFaultData; using Tgov1DataT = Governor::Tgov1Data; using Ieeet1DataT = Exciter::Ieeet1Data; + using IeeestDataT = Stabilizer::IeeestData; using GenrouDataT = GenrouData; using GenClassicalDataT = GenClassicalData; using LoadDataT = LoadData; @@ -88,6 +90,7 @@ namespace GridKit std::vector loadzip; ///< Loads within the model std::vector gov; ///< Governors within the model std::vector exciter; ///< Exciters within the model + std::vector stabilizer; ///< Stabilizers within the model std::vector signal; ///< Signal nodes /// Monitor sink specs diff --git a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp index 81556d6b1..35f516764 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp @@ -136,6 +136,12 @@ namespace GridKit raw_component.get_to(exciter); sm.exciter.push_back(exciter); } + else if (kind == "Ieeest") + { + typename SystemModelData::IeeestDataT stabilizer; + raw_component.get_to(stabilizer); + sm.stabilizer.push_back(stabilizer); + } else if (kind == "BusFault") { typename SystemModelData::BusFaultDataT bus_fault; diff --git a/docs/Figures/stabilizer_ieeest_diagram.png b/docs/Figures/stabilizer_ieeest_diagram.png new file mode 100644 index 000000000..7771d6d6c Binary files /dev/null and b/docs/Figures/stabilizer_ieeest_diagram.png differ diff --git a/tests/UnitTests/PhasorDynamics/CMakeLists.txt b/tests/UnitTests/PhasorDynamics/CMakeLists.txt index ddf828ce1..687dd82d9 100644 --- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt +++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt @@ -50,6 +50,13 @@ target_link_libraries(test_phasor_governortgov1 GridKit::definitions GridKit::phasor_dynamics_bus_dependency_tracking GridKit::testing) +add_executable(test_phasor_stabilizer_ieeest runStabilizerIeeestTests.cpp) +target_link_libraries(test_phasor_stabilizer_ieeest GridKit::definitions + GridKit::phasor_dynamics_stabilizer_ieeest + GridKit::phasor_dynamics_stabilizer_ieeest_dependency_tracking + GridKit::phasor_dynamics_signal + GridKit::testing) + add_executable(test_phasor_gen_classical runGenClassicalTests.cpp) target_link_libraries(test_phasor_gen_classical GridKit::definitions GridKit::phasor_dynamics_gen_classical @@ -77,6 +84,7 @@ 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 PhasorDynamicsGovernorTgov1Test COMMAND test_phasor_governortgov1) +add_test(NAME PhasorDynamicsStabilizerIeeestTest COMMAND test_phasor_stabilizer_ieeest) add_test(NAME PhasorDynamicsGenClassicalTest COMMAND test_phasor_gen_classical) add_test(NAME PhasorDynamicsLoadTest COMMAND test_phasor_load) add_test(NAME PhasorDynamicsLoadZIPTest COMMAND test_phasor_loadzip) @@ -90,6 +98,7 @@ install(TARGETS test_phasor_bus test_phasor_loadzip test_phasor_genrou test_phasor_governortgov1 + test_phasor_stabilizer_ieeest test_phasor_gen_classical test_phasor_system test_phasor_system_single_component diff --git a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp new file mode 100644 index 000000000..d35d774a6 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp @@ -0,0 +1,418 @@ +#pragma once + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + template + class StabilizerIeeestTests + { + public: + using RealT = typename PhasorDynamics::Component::RealT; + + StabilizerIeeestTests() = default; + ~StabilizerIeeestTests() = default; + + TestOutcome constructor() + { + TestStatus success = true; + + auto data = makeTestData(); + auto* stab = new PhasorDynamics::Stabilizer::Ieeest(data); + + success *= (stab != nullptr); + success *= (stab->getMonitor() != nullptr); + + delete stab; + + return success.report(__func__); + } + + /** + * @brief All states initialize to zero (stabilizer at rest). + * With u = 0, all residuals should be zero. + */ + TestOutcome zeroInitialResidual() + { + TestStatus success = true; + + // Create signal nodes for input (u) and output (Vs) + PhasorDynamics::SignalNode u_node; + PhasorDynamics::SignalNode vs_node; + ScalarT u_value{0.0}; + IdxT u_index = 13; // beyond internal variables + ScalarT vs_value{0.0}; + IdxT vs_index = INVALID_INDEX; + + // Link signal nodes to backing storage + u_node.set(&u_value, &u_index); + vs_node.set(&vs_value, &vs_index); + + auto data = makeTestData(); + PhasorDynamics::Stabilizer::Ieeest stab(data); + + // Wire: stabilizer reads u_node as input, writes vs_node as output + stab.getSignals().template attachSignalNode(&u_node); + stab.getSignals().template assignSignalNode(&vs_node); + + stab.allocate(); + success *= (stab.verify() == 0); + stab.initialize(); + stab.evaluateResidual(); + + auto tol = 10 * std::numeric_limits::epsilon(); + const auto& f = stab.getResidual(); + for (size_t i = 0; i < f.size(); ++i) + { + if (!isEqual(f[i], 0.0, tol)) + { + std::cout << "Non-zero residual at index " << i << ": " << f[i] << "\n"; + success = false; + } + } + + // Verify output signal is linked and reads the correct value + success *= vs_node.linked(); + success *= (vs_node.getVariableIndex() == 12); + success *= isEqual(vs_node.read(), static_cast(0.0), tol); + + return success.report(__func__); + } + + /** + * @brief Residual evaluation against hand-computed answer key. + * + * Sets specific y/yp values and verifies residuals match + * pre-computed expected values. See plan for derivation. + */ + TestOutcome residual() + { + TestStatus success = true; + + PhasorDynamics::SignalNode u_node; + PhasorDynamics::SignalNode vs_node; + ScalarT u_value{0.5}; + IdxT u_index = 13; + ScalarT vs_value{0.0}; + IdxT vs_index = INVALID_INDEX; + + u_node.set(&u_value, &u_index); + vs_node.set(&vs_value, &vs_index); + + auto data = makeTestData(); + PhasorDynamics::Stabilizer::Ieeest stab(data); + + stab.getSignals().template attachSignalNode(&u_node); + stab.getSignals().template assignSignalNode(&vs_node); + + stab.allocate(); + stab.initialize(); + setStatePoint(stab); + stab.evaluateResidual(); + + // Hand-computed answer key (see plan for full derivation) + const std::vector res_answer = { + 0.19, // f[0]: -x1_dot + x2 + 0.28, // f[1]: -x2_dot + x3 + 0.37, // f[2]: -x3_dot + x4 + 1.0975, // f[3]: -x4_dot + (-a0*x1 - a1*x2 - a2*x3 - a3*x4 + u) / a4 + 0.25, // f[4]: -x5_dot + (v4 - x5) / T2 + 0.24, // f[5]: -x6_dot + (v5 - x6) / T4 + -0.01, // f[6]: -x7_dot + (v6 - x7) / T6 + -0.42, // f[7]: -v4 + x1 + A5*x2 + A6*x3 + -0.25, // f[8]: -v5 + x5 + (T1/T2)*(v4 - x5) + -0.31, // f[9]: -v6 + x6 + (T3/T4)*(v5 - x6) + 1.15, // f[10]: -v7 + Ks*(T5/T6)*(v6 - x7) + 0.0, // f[11]: limiter (v7=0.05 within [-0.1, 0.1]) + 0.0, // f[12]: cutout (Vct=0.75 within [0.0, 1.5]) + }; + + // Looser tolerance for f[11],f[12] — smooth sigmoid approximations + const auto loose_tol = static_cast(1.0e-4); + auto& residual = stab.getResidual(); + + for (size_t i = 0; i < res_answer.size(); ++i) + { + auto test_tol = (i >= 11) ? loose_tol : static_cast(10 * std::numeric_limits::epsilon()); + if (!isEqual(residual[i], res_answer[i], test_tol)) + { + std::cout << "Incorrect result for residual " << i << ": " + << std::setprecision(15) << residual[i] + << " != " << res_answer[i] << "\n"; + success = false; + } + } + + // Verify output signal reads the stabilizer output + success *= isEqual(vs_node.read(), static_cast(0.05), loose_tol); + + return success.report(__func__); + } + +#ifdef GRIDKIT_ENABLE_ENZYME + /** + * @brief Compare DependencyTracking Jacobian against Enzyme Jacobian. + */ + TestOutcome jacobian() + { + TestStatus success = true; + + auto data = makeTestData(); + + std::vector + dependency_tracking_jacobian = DependencyTrackingJacobian(data); + + std::vector + enzyme_jacobian = EnzymeJacobian(data); + + // Compare DependencyTracking dependencies to Enzyme's + auto tol = 10 * std::numeric_limits::epsilon(); + 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( + PhasorDynamics::Stabilizer::IeeestData ieeestdata) + { + using DepVar = DependencyTracking::Variable; + + // Set up signal nodes with DependencyTracking scalar type + PhasorDynamics::SignalNode u_node; + PhasorDynamics::SignalNode vs_node; + DepVar u_value{0.5}; + IdxT u_index = 13; + DepVar vs_value{0.0}; + IdxT vs_index = INVALID_INDEX; + + u_node.set(&u_value, &u_index); + vs_node.set(&vs_value, &vs_index); + + PhasorDynamics::Stabilizer::Ieeest stab(ieeestdata); + stab.getSignals().template attachSignalNode(&u_node); + stab.getSignals().template assignSignalNode(&vs_node); + + stab.allocate(); + stab.initialize(); + + // --- d/dy: tag internal variables as independent --- + for (size_t i = 0; i < stab.size(); ++i) + { + stab.y()[i].setVariableNumber(i); + } + // Tag external signal u as an additional independent variable + u_value.setVariableNumber(stab.size()); + u_value.setValue(0.5); + + setStatePointDep(stab); + + stab.evaluateResidual(); + std::vector residual_y = stab.getResidual(); + + // --- d/dy': tag derivatives as independent --- + stab.initialize(); + for (size_t i = 0; i < stab.size(); ++i) + { + stab.yp()[i].setVariableNumber(i); + } + + u_value = 0.5; + setStatePointDep(stab); + + stab.evaluateResidual(); + std::vector residual_yp = stab.getResidual(); + + // Print dependencies for debugging + for (size_t i = 0; i < residual_y.size(); ++i) + { + std::cout << i << "th residual, y: "; + (residual_y[i]).print(std::cout); + std::cout << "\n"; + std::cout << i << "th residual, yp: "; + (residual_yp[i]).print(std::cout); + std::cout << "\n"; + } + + // Merge d/dy and d/dy' into a single dependency map + std::vector dependencies(residual_y.size()); + for (IdxT i = 0; i < residual_y.size(); ++i) + { + auto dependency_y = (residual_y[i]).getDependencies(); + auto dependency_yp = (residual_yp[i]).getDependencies(); + + for (const auto& pair_y : dependency_y) + { + auto it_yp = dependency_yp.find(pair_y.first); + if (it_yp != dependency_yp.end()) + { + dependencies[i].insert(std::make_pair(pair_y.first, pair_y.second + it_yp->second)); + } + else + { + dependencies[i].insert(std::make_pair(pair_y.first, pair_y.second)); + } + } + + // Insert yp dependencies that did not exist in the y dependencies + for (const auto& pair_yp : dependency_yp) + { + if (!dependency_y.contains(pair_yp.first)) + { + dependencies[i].insert(std::make_pair(pair_yp.first, pair_yp.second)); + } + } + } + + return dependencies; + } + + std::vector EnzymeJacobian( + PhasorDynamics::Stabilizer::IeeestData ieeestdata) + { + PhasorDynamics::SignalNode u_node; + PhasorDynamics::SignalNode vs_node; + ScalarT u_value{0.5}; + IdxT u_index = 13; + ScalarT vs_value{0.0}; + IdxT vs_index = INVALID_INDEX; + + u_node.set(&u_value, &u_index); + vs_node.set(&vs_value, &vs_index); + + PhasorDynamics::Stabilizer::Ieeest stab(ieeestdata); + stab.getSignals().template attachSignalNode(&u_node); + stab.getSignals().template assignSignalNode(&vs_node); + + stab.allocate(); + stab.initialize(); + setStatePoint(stab); + + stab.updateTime(0.0, 1.0); // alpha = 1.0 to verify d/dy' term + + stab.evaluateResidual(); + stab.evaluateJacobian(); + + auto model_jacobian = stab.getJacobian(); + model_jacobian.deduplicate(); + model_jacobian.printMatrix("IEEEST Jacobian"); + + return GridKit::Testing::MapFromCOO(model_jacobian); + } +#endif + + private: + static constexpr ScalarT tol_ = 10 * std::numeric_limits::epsilon(); + + /** + * @brief Standard IEEEST parameter set for all tests. + * Derived: a0=1, a1=0.4, a2=0.63, a3=0.1, a4=0.08 + */ + auto makeTestData() -> PhasorDynamics::Stabilizer::IeeestData + { + using Params = PhasorDynamics::Stabilizer::IeeestParameters; + + PhasorDynamics::Stabilizer::IeeestData data; + data.device_class = "stabilizer"; + data.disambiguation_string = "ieeest_test"; + data.monitored_variables.insert(PhasorDynamics::Stabilizer::IeeestMonitorableVariables::vs); + + data.parameters[Params::A1] = 0.1; + data.parameters[Params::A2] = 0.2; + data.parameters[Params::A3] = 0.3; + data.parameters[Params::A4] = 0.4; + data.parameters[Params::A5] = 0.5; + data.parameters[Params::A6] = 0.6; + data.parameters[Params::T1] = 0.5; + data.parameters[Params::T2] = 1.0; + data.parameters[Params::T3] = 0.3; + data.parameters[Params::T4] = 1.0; + data.parameters[Params::T5] = 2.0; + data.parameters[Params::T6] = 5.0; + data.parameters[Params::Ks] = 10.0; + data.parameters[Params::Lsmin] = -0.1; + data.parameters[Params::Lsmax] = 0.1; + data.parameters[Params::Vcl] = 0.0; + data.parameters[Params::Vcu] = 1.5; + data.parameters[Params::Tdelay] = 0.0; + + return data; + } + + /** + * @brief Set a non-trivial operating point for residual/Jacobian tests. + * Avoids zeros and ones to catch coefficient errors. + */ + void setStatePoint(PhasorDynamics::Stabilizer::Ieeest& stab) + { + stab.y()[0] = 0.1; // x1 + stab.y()[1] = 0.2; // x2 + stab.y()[2] = 0.3; // x3 + stab.y()[3] = 0.4; // x4 + stab.y()[4] = 0.5; // x5 + stab.y()[5] = 0.6; // x6 + stab.y()[6] = 0.7; // x7 + stab.y()[7] = 0.8; // v4 + stab.y()[8] = 0.9; // v5 + stab.y()[9] = 1.0; // v6 + stab.y()[10] = 0.05; // v7 (within limiter range) + stab.y()[11] = 0.05; // Vss + stab.y()[12] = 0.05; // Vs + + stab.yp()[0] = 0.01; // x1_dot + stab.yp()[1] = 0.02; // x2_dot + stab.yp()[2] = 0.03; // x3_dot + stab.yp()[3] = 0.04; // x4_dot + stab.yp()[4] = 0.05; // x5_dot + stab.yp()[5] = 0.06; // x6_dot + stab.yp()[6] = 0.07; // x7_dot + } + + /** + * @brief Set the same operating point for DependencyTracking variables. + * Uses setValue() to set the numeric value while preserving dependency info. + */ + void setStatePointDep(PhasorDynamics::Stabilizer::Ieeest& stab) + { + stab.y()[0].setValue(0.1); + stab.y()[1].setValue(0.2); + stab.y()[2].setValue(0.3); + stab.y()[3].setValue(0.4); + stab.y()[4].setValue(0.5); + stab.y()[5].setValue(0.6); + stab.y()[6].setValue(0.7); + stab.y()[7].setValue(0.8); + stab.y()[8].setValue(0.9); + stab.y()[9].setValue(1.0); + stab.y()[10].setValue(0.05); + stab.y()[11].setValue(0.05); + stab.y()[12].setValue(0.05); + + stab.yp()[0].setValue(0.01); + stab.yp()[1].setValue(0.02); + stab.yp()[2].setValue(0.03); + stab.yp()[3].setValue(0.04); + stab.yp()[4].setValue(0.05); + stab.yp()[5].setValue(0.06); + stab.yp()[6].setValue(0.07); + } + }; // class StabilizerIeeestTests + + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/runStabilizerIeeestTests.cpp b/tests/UnitTests/PhasorDynamics/runStabilizerIeeestTests.cpp new file mode 100644 index 000000000..c1999c24a --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/runStabilizerIeeestTests.cpp @@ -0,0 +1,17 @@ +#include "StabilizerIeeestTests.hpp" + +int main() +{ + GridKit::Testing::TestingResults result; + + GridKit::Testing::StabilizerIeeestTests test; + + result += test.constructor(); + result += test.zeroInitialResidual(); + result += test.residual(); +#ifdef GRIDKIT_ENABLE_ENZYME + result += test.jacobian(); +#endif + + return result.summary(); +}