From 73cd093005091b41ef1bd89c18e1e166bababae8 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Sat, 8 Aug 2026 18:01:09 -0500 Subject: [PATCH 1/6] IEEEST Stabilizer Improvements/Corrections --- .../Stabilizer/IEEEST/CMakeLists.txt | 5 +- .../Stabilizer/IEEEST/Ieeest.hpp | 129 +- .../Stabilizer/IEEEST/IeeestData.hpp | 38 +- .../Stabilizer/IEEEST/IeeestEnzyme.cpp | 100 +- .../Stabilizer/IEEEST/IeeestImpl.hpp | 540 +++++--- .../Stabilizer/IEEEST/README.md | 267 +++- .../Model/PhasorDynamics/SystemModelImpl.hpp | 45 +- .../IEEEST/diagram.png} | Bin .../ComponentConnectionTests.hpp | 65 + .../PhasorDynamics/StabilizerIeeestTests.hpp | 1115 ++++++++++++----- .../runComponentConnectionTests.cpp | 1 + .../runStabilizerIeeestTests.cpp | 8 +- 12 files changed, 1575 insertions(+), 738 deletions(-) rename docs/Figures/{stabilizer_ieeest_diagram.png => PhasorDynamics/IEEEST/diagram.png} (100%) diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/CMakeLists.txt index 936a29ea1..0608d1772 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/CMakeLists.txt @@ -20,8 +20,9 @@ if(GRIDKIT_ENABLE_ENZYME) ClangEnzymeFlags COMPILE_OPTIONS PRIVATE - -mllvm - -enzyme-auto-sparsity=1 + "SHELL:-mllvm -enzyme-auto-sparsity=1" + # Preserve runtime-order stores for Enzyme sparsification. + "SHELL:-mllvm -simplifycfg-sink-common=false" -fno-math-errno) else() gridkit_add_library( diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp index f77935130..14e4494da 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp @@ -6,92 +6,85 @@ #pragma once +#include +#include +#include + +#include #include #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` + /// 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 (model output) + X1, ///< \f$x_1\f$ Notch-filter signal state + X2, ///< \f$x_2\f$ First derivative of the filtered signal + X3, ///< \f$x_3\f$ Second derivative of the filtered signal + X4, ///< \f$x_4\f$ Third derivative of the filtered signal + X5, ///< \f$x_5\f$ Lead-lag 1 state + X6, ///< \f$x_6\f$ Lead-lag 2 state + X7, ///< \f$x_7\f$ Washout state + V4, ///< \f$v_4\f$ Notch-filter output + V5, ///< \f$v_5\f$ Lead-lag 1 output + V6, ///< \f$v_6\f$ Lead-lag 2 output + V7, ///< \f$v_7\f$ Unlimited stabilizer signal + VSS, ///< \f$V_{\mathrm{ss}}\f$ Limited stabilizer signal and model output MAXIMUM, }; - /// External variables of a `Ieeest` + /// External variables of a `Ieeest`. enum class IeeestExternalVariables : size_t { - U, ///< Stabilizer input signal + U, ///< \f$u\f$ Stabilizer input signal MAXIMUM, }; template class Ieeest : public Component { - using Component::gridkit_component_id_; + using Component::abs_tol_; + using Component::allocated_; using Component::alpha_; using Component::f_; + using Component::gridkit_component_id_; + using Component::J_cols_buffer_; + using Component::J_rows_buffer_; + using Component::J_vals_buffer_; using Component::nnz_; + using Component::residual_indices_; using Component::size_; using Component::tag_; - using Component::abs_tol_; - using Component::time_; + using Component::variable_indices_; + using Component::wb_; using Component::y_; using Component::yp_; - using Component::wb_; - using Component::h_; - using Component::J_rows_buffer_; - using Component::J_cols_buffer_; - using Component::J_vals_buffer_; - using Component::variable_indices_; - using Component::residual_indices_; - using Component::allocated_; public: - using ScalarT = scalar_type; - using IdxT = index_type; - using RealT = typename Component::RealT; - using ModelDataT = IeeestData; - using SignalT = SignalNode; - using MonitorT = Model::VariableMonitor; + using ScalarT = scalar_type; + using IdxT = index_type; + using RealT = typename Component::RealT; + using SignalT = SignalNode; + using ModelDataT = IeeestData; + using MonitorT = Model::VariableMonitor; + using InternalVariablesT = IeeestInternalVariables; + using ExternalVariablesT = IeeestExternalVariables; Ieeest(); - Ieeest(const ModelDataT& data); + explicit Ieeest(const ModelDataT& data); ~Ieeest(); - int setGridKitComponentID(IdxT) override final; + int setGridKitComponentID(IdxT component_id) override final; int allocate() override final; int verify() const override final; int initialize() override final; @@ -100,7 +93,6 @@ namespace GridKit int evaluateResidual() override final; int evaluateJacobian() override final; - /// Get the `ComponentSignals` from this `Ieeest` auto getSignals() -> ComponentSignals(1.0e-3); + RealT A1_{0}; RealT A2_{0}; RealT A3_{0}; @@ -139,34 +137,17 @@ namespace GridKit RealT Vcu_{0}; RealT Tdelay_{0}; - RealT a0_{1}; + IdxT order_{0}; + 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_notch_{0}; - RealT bypass_notch_{1}; - 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 use_T4_block_{1}; - RealT bypass_T4_block_{0}; - RealT use_T6_block_{1}; - RealT bypass_T6_block_{0}; + IdxT parameter_error_count_{0}; ComponentSignals signals_; - - std::unique_ptr monitor_; - - void initializeParameters(const ModelDataT& data); - void initializeMonitor(); + std::unique_ptr monitor_; std::vector ws_; std::vector ws_indices_; diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestData.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestData.hpp index 139e928af..ba8a1466b 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestData.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestData.hpp @@ -19,24 +19,24 @@ namespace GridKit */ 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 (not modeled) - Vcu, ///< Upper input cutout threshold (not modeled) - Tdelay, ///< Input time delay (not modeled) + A1, ///< \f$A_1\f$ Notch-filter denominator coefficient + A2, ///< \f$A_2\f$ Notch-filter denominator coefficient + A3, ///< \f$A_3\f$ Notch-filter denominator coefficient + A4, ///< \f$A_4\f$ Notch-filter denominator coefficient + A5, ///< \f$A_5\f$ Notch-filter numerator coefficient + A6, ///< \f$A_6\f$ Notch-filter numerator coefficient + T1, ///< \f$T_1\f$ Lead-lag 1 numerator time constant + T2, ///< \f$T_2\f$ Lead-lag 1 denominator time constant + T3, ///< \f$T_3\f$ Lead-lag 2 numerator time constant + T4, ///< \f$T_4\f$ Lead-lag 2 denominator time constant + T5, ///< \f$T_5\f$ Washout numerator time constant + T6, ///< \f$T_6\f$ Washout denominator time constant + Ks, ///< \f$K_s\f$ Stabilizer gain + Lsmin, ///< \f$L_s^{\min}\f$ Minimum stabilizer output limit + Lsmax, ///< \f$L_s^{\max}\f$ Maximum stabilizer output limit + Vcl, ///< \f$V_{\mathrm{cl}}\f$ Unsupported lower input-cutout threshold + Vcu, ///< \f$V_{\mathrm{cu}}\f$ Unsupported upper input-cutout threshold + Tdelay, ///< \f$T_{\mathrm{delay}}\f$ Unsupported input delay }; /** @@ -70,7 +70,7 @@ namespace GridKit */ enum class IeeestMonitorableVariables { - vss, ///< Stabilizer output (limited signal) + vss, ///< \f$V_{\mathrm{ss}}\f$ Limited stabilizer signal and model output }; /** diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp index 4c1af73d9..673dfa1e5 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp @@ -1,7 +1,7 @@ /** * @file IeeestEnzyme.cpp * @author Luke Lowery (lukel@tamu.edu) - * @brief Enzyme-based sparse Jacobian for IEEEST Stabilizer. + * @brief Enzyme sparse Jacobian for the IEEEST stabilizer model. */ #include @@ -14,13 +14,6 @@ namespace GridKit { 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() { @@ -29,9 +22,6 @@ namespace GridKit if (J_rows_buffer_ == nullptr) { - // Reserve space for the dense blocks. - // The size of the buffer is the sum of maximum capacities of the blocks. - // Enyme will compute the appropriate nnz from sparsification. auto size = static_cast(size_); auto signal_size = static_cast(ws_.size()); auto buffer_size = 2 * size * size + size * signal_size; @@ -40,63 +30,61 @@ namespace GridKit J_vals_buffer_ = new RealT[buffer_size]; } + using ModelT = GridKit::PhasorDynamics::Stabilizer::Ieeest; + using Fn = GridKit::Enzyme::Sparse::MemberFunctions; + nnz_ = 0; - GridKit::Enzyme::Sparse::DfDy, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(y_.getSize()), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.data(), - ws_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDy::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); - GridKit::Enzyme::Sparse::DfDyp, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - static_cast(y_.getSize()), - (this->getResidualIndices()).data(), - (this->getVariableIndices()).data(), - y_.getData(), - yp_.getData(), - wb_.data(), - ws_.data(), - alpha_, - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDyp::eval(this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.getData(), + yp_.getData(), + wb_.data(), + ws_.data(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); - GridKit::Enzyme::Sparse::DfDws, - GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this, - static_cast(f_.getSize()), - ws_.size(), - (this->getResidualIndices()).data(), - ws_indices_.data(), - y_.getData(), - yp_.getData(), - wb_.data(), - ws_.data(), - J_rows_buffer_, - J_cols_buffer_, - J_vals_buffer_, - nnz_); + GridKit::Enzyme::Sparse::DfDws::eval(this, + static_cast(f_.getSize()), + ws_.size(), + (this->getResidualIndices()).data(), + ws_indices_.data(), + y_.getData(), + yp_.getData(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); this->constructCoo(); 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 index 970907ab0..1c2b3fc6f 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp @@ -6,7 +6,9 @@ * @brief Definition of the IEEEST Power System Stabilizer. */ -#include +#include +#include +#include #include #include @@ -25,7 +27,8 @@ namespace GridKit template Ieeest::Ieeest() { - size_ = 12; + size_ = static_cast(IeeestInternalVariables::MAXIMUM); + setDerivedParameters(); } template @@ -34,7 +37,7 @@ namespace GridKit { initializeParameters(data); initializeMonitor(); - size_ = 12; + size_ = static_cast(IeeestInternalVariables::MAXIMUM); } template @@ -42,110 +45,6 @@ namespace GridKit { } - template - void Ieeest::initializeParameters(const ModelDataT& data) - { - using Parameter = typename ModelDataT::Parameters; - if (data.parameters.contains(Parameter::A1)) - { - A1_ = std::get(data.parameters.at(Parameter::A1)); - } - if (data.parameters.contains(Parameter::A2)) - { - A2_ = std::get(data.parameters.at(Parameter::A2)); - } - if (data.parameters.contains(Parameter::A3)) - { - A3_ = std::get(data.parameters.at(Parameter::A3)); - } - if (data.parameters.contains(Parameter::A4)) - { - A4_ = std::get(data.parameters.at(Parameter::A4)); - } - if (data.parameters.contains(Parameter::A5)) - { - A5_ = std::get(data.parameters.at(Parameter::A5)); - } - if (data.parameters.contains(Parameter::A6)) - { - A6_ = std::get(data.parameters.at(Parameter::A6)); - } - if (data.parameters.contains(Parameter::T1)) - { - T1_ = std::get(data.parameters.at(Parameter::T1)); - } - if (data.parameters.contains(Parameter::T2)) - { - T2_ = std::get(data.parameters.at(Parameter::T2)); - } - if (data.parameters.contains(Parameter::T3)) - { - T3_ = std::get(data.parameters.at(Parameter::T3)); - } - if (data.parameters.contains(Parameter::T4)) - { - T4_ = std::get(data.parameters.at(Parameter::T4)); - } - if (data.parameters.contains(Parameter::T5)) - { - T5_ = std::get(data.parameters.at(Parameter::T5)); - } - if (data.parameters.contains(Parameter::T6)) - { - T6_ = std::get(data.parameters.at(Parameter::T6)); - } - if (data.parameters.contains(Parameter::Ks)) - { - Ks_ = std::get(data.parameters.at(Parameter::Ks)); - } - if (data.parameters.contains(Parameter::Lsmin)) - { - Lsmin_ = std::get(data.parameters.at(Parameter::Lsmin)); - } - if (data.parameters.contains(Parameter::Lsmax)) - { - Lsmax_ = std::get(data.parameters.at(Parameter::Lsmax)); - } - if (data.parameters.contains(Parameter::Vcl)) - { - Vcl_ = std::get(data.parameters.at(Parameter::Vcl)); - } - if (data.parameters.contains(Parameter::Vcu)) - { - Vcu_ = std::get(data.parameters.at(Parameter::Vcu)); - } - if (data.parameters.contains(Parameter::Tdelay)) - { - Tdelay_ = std::get(data.parameters.at(Parameter::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_notch_ = static_cast(a2_ != 0.0 || a3_ != 0.0 || a4_ != 0.0); - bypass_notch_ = 1.0 - use_notch_; - - 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_; - - use_T4_block_ = static_cast(T4_ != 0.0); - bypass_T4_block_ = 1.0 - use_T4_block_; - - use_T6_block_ = static_cast(T6_ != 0.0); - bypass_T6_block_ = 1.0 - use_T6_block_; - } - template int Ieeest::setGridKitComponentID(IdxT component_id) { @@ -160,28 +59,30 @@ namespace GridKit { this->allocateVectors(size_); } - auto size = static_cast(size_); - - tag_.resize(size); + const auto size = static_cast(size_); + tag_.assign(size, false); variable_indices_.resize(size); residual_indices_.resize(size); + + wb_.clear(); + + const auto signal_size = static_cast(IeeestExternalVariables::MAXIMUM); + ws_.assign(signal_size, ScalarT{0}); + ws_indices_.assign(signal_size, INVALID_INDEX); + for (IdxT j = 0; j < size_; ++j) { this->setVariableIndex(j, j); this->setResidualIndex(j, j); } - ws_.resize(1); - ws_indices_.resize(1); - ws_[0] = 0.0; - ws_indices_[0] = INVALID_INDEX; - if (signals_.template isAssigned()) { auto* y = y_.getData(); signals_.template getSignalNode()->set( - &y[11], &(this->getVariableIndex(11))); + &y[static_cast(IeeestInternalVariables::VSS)], + &(this->getVariableIndex(static_cast(IeeestInternalVariables::VSS)))); } allocated_ = true; @@ -191,25 +92,69 @@ namespace GridKit template int Ieeest::verify() const { - int ret = 0; + int ret = static_cast(parameter_error_count_); - if (signals_.template isAttached()) + auto check = [&](bool condition, const char* message) { - if (!signals_.template isLinked()) + if (!condition) { - Log::error() << "Ieeest: input signal U attached with no linked source\n"; + Log::error() << "Ieeest: " << message << '\n'; ret += 1; } - } - else + }; + + if (!signals_.template isAttached()) { Log::error() << "Ieeest: required input signal U is not attached\n"; ret += 1; } + else if (!signals_.template isLinked()) + { + Log::error() << "Ieeest: input signal U attached with no linked source\n"; + ret += 1; + } + + check(std::isfinite(a1_) && std::isfinite(a2_) + && std::isfinite(a3_) && std::isfinite(a4_), + "expanded denominator coefficients must be finite"); + + const IdxT numerator_order = A6_ != ZERO + ? static_cast(2) + : A5_ != ZERO + ? static_cast(1) + : static_cast(0); + check(numerator_order <= order_, + "numerator order must not exceed denominator order"); + check(Lsmin_ < Lsmax_, "Lsmin must be less than Lsmax"); + + RealT leading_coefficient = ONE; + const char* leading_name = nullptr; + switch (order_) + { + case 0: + break; + case 1: + leading_coefficient = a1_; + leading_name = "a1"; + break; + case 2: + leading_coefficient = a2_; + leading_name = "a2"; + break; + case 3: + leading_coefficient = a3_; + leading_name = "a3"; + break; + case 4: + leading_coefficient = a4_; + leading_name = "a4"; + break; + } - if (a4_ == 0 && a3_ == 0 && a2_ == 0 && a1_ != 0) + if (leading_name != nullptr && leading_coefficient == ZERO) { - Log::error() << "Ieeest: a2, a3, and a4 are all zero - no valid notch filter\n"; + Log::error() << "Ieeest: active leading denominator coefficient '" + << leading_name << "' must be nonzero\n"; ret += 1; } @@ -219,15 +164,59 @@ namespace GridKit template int Ieeest::initialize() { + const auto X1 = static_cast(IeeestInternalVariables::X1); + const auto X5 = static_cast(IeeestInternalVariables::X5); + const auto X6 = static_cast(IeeestInternalVariables::X6); + const auto X7 = static_cast(IeeestInternalVariables::X7); + const auto V4 = static_cast(IeeestInternalVariables::V4); + const auto V5 = static_cast(IeeestInternalVariables::V5); + const auto V6 = static_cast(IeeestInternalVariables::V6); + const auto V7 = static_cast(IeeestInternalVariables::V7); + const auto VSS = static_cast(IeeestInternalVariables::VSS); + const auto U = static_cast(IeeestExternalVariables::U); + + if (verify() != 0) + { + Log::error() << "Ieeest: cannot initialize with invalid configuration\n"; + return 1; + } + + const ScalarT u = signals_.template readExternalVariable(); + if (!std::isfinite(static_cast(u))) + { + Log::error() << "Ieeest: initial input signal U must be finite\n"; + return 1; + } + auto* y = y_.getData(); auto* yp = yp_.getData(); + std::fill_n(y, static_cast(size_), ScalarT{ZERO}); + std::fill_n(yp, static_cast(size_), ScalarT{ZERO}); - for (IdxT i = 0; i < size_; ++i) + ws_[U] = u; + ws_indices_[U] = signals_.template readExternalVariableIndex(); + + switch (order_) { - y[static_cast(i)] = 0.0; - yp[static_cast(i)] = 0.0; + case 0: + break; + case 1: + case 2: + case 3: + case 4: + y[X1] = u; + break; } + y[X5] = u; + y[X6] = u; + y[X7] = u; + y[V4] = u; + y[V5] = u; + y[V6] = u; + y[V7] = ZERO; + y[VSS] = Math::clamp(y[V7], Lsmin_, Lsmax_); + y_.setDataUpdated(); yp_.setDataUpdated(); @@ -237,18 +226,9 @@ namespace GridKit template int Ieeest::tagDifferentiable() { - tag_[0] = true; - tag_[1] = true; - tag_[2] = true; - tag_[3] = true; - tag_[4] = (T2_ != 0.0); - tag_[5] = (T4_ != 0.0); - tag_[6] = (T6_ != 0.0); - tag_[7] = false; - tag_[8] = false; - tag_[9] = false; - tag_[10] = false; - tag_[11] = false; + std::fill(tag_.begin(), tag_.end(), false); + constexpr auto X7 = static_cast(IeeestInternalVariables::X7); + std::fill_n(tag_.begin(), X7 + 1, true); return 0; } @@ -280,43 +260,106 @@ namespace GridKit const 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 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]; - - f[0] = -x1_dot + use_notch_ * 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] = -T2_ * x5_dot - x5 + v4; - f[5] = -T4_ * x6_dot - x6 + v5; - f[6] = -T6_ * x7_dot - x7 + v6; - f[7] = -v4 + bypass_notch_ * u + use_notch_ * (x1 + A5_ * x2 + (use_4th_order_ + use_3rd_order_) * A6_ * x3); - f[8] = use_T2_block_ * (-T2_ * (v5 - x5) + T1_ * (v4 - x5)) + bypass_T2_block_ * (v4 - v5); - f[9] = use_T4_block_ * (-T4_ * (v6 - x6) + T3_ * (v5 - x6)) + bypass_T4_block_ * (v5 - v6); - f[10] = use_T6_block_ * (-T6_ * v7 + Ks_ * T5_ * (v6 - x7)) + bypass_T6_block_ * (Ks_ * v6 - v7); - f[11] = -vss + Math::clamp(v7, Lsmin_, Lsmax_); + const auto X1 = static_cast(IeeestInternalVariables::X1); + const auto X2 = static_cast(IeeestInternalVariables::X2); + const auto X3 = static_cast(IeeestInternalVariables::X3); + const auto X4 = static_cast(IeeestInternalVariables::X4); + const auto X5 = static_cast(IeeestInternalVariables::X5); + const auto X6 = static_cast(IeeestInternalVariables::X6); + const auto X7 = static_cast(IeeestInternalVariables::X7); + const auto V4 = static_cast(IeeestInternalVariables::V4); + const auto V5 = static_cast(IeeestInternalVariables::V5); + const auto V6 = static_cast(IeeestInternalVariables::V6); + const auto V7 = static_cast(IeeestInternalVariables::V7); + const auto VSS = static_cast(IeeestInternalVariables::VSS); + const auto U = static_cast(IeeestExternalVariables::U); + + const ScalarT x1 = y[X1]; + const ScalarT x2 = y[X2]; + const ScalarT x3 = y[X3]; + const ScalarT x4 = y[X4]; + const ScalarT x5 = y[X5]; + const ScalarT x6 = y[X6]; + const ScalarT x7 = y[X7]; + const ScalarT v4 = y[V4]; + const ScalarT v5 = y[V5]; + const ScalarT v6 = y[V6]; + const ScalarT v7 = y[V7]; + const ScalarT vss = y[VSS]; + + const ScalarT x1_dot = yp[X1]; + const ScalarT x2_dot = yp[X2]; + const ScalarT x3_dot = yp[X3]; + const ScalarT x4_dot = yp[X4]; + const ScalarT x5_dot = yp[X5]; + const ScalarT x6_dot = yp[X6]; + const ScalarT x7_dot = yp[X7]; + + const ScalarT u = ws[U]; + + switch (order_) + { + case 0: + f[X1] = -x1_dot; + f[X2] = -x2_dot; + f[X3] = -x3_dot; + f[X4] = -x4_dot; + f[V4] = -v4 + u; + break; + case 1: + { + const ScalarT x1_rhs = (u - x1) / a1_; + f[X1] = -x1_dot + x1_rhs; + f[X2] = -x2_dot; + f[X3] = -x3_dot; + f[X4] = -x4_dot; + f[V4] = -v4 + x1 + A5_ * x1_rhs; + break; + } + case 2: + { + const ScalarT x2_rhs = (u - x1 - a1_ * x2) / a2_; + f[X1] = -x1_dot + x2; + f[X2] = -x2_dot + x2_rhs; + f[X3] = -x3_dot; + f[X4] = -x4_dot; + f[V4] = -v4 + x1 + A5_ * x2 + A6_ * x2_rhs; + break; + } + case 3: + { + const ScalarT x3_rhs = (u - x1 - a1_ * x2 - a2_ * x3) / a3_; + f[X1] = -x1_dot + x2; + f[X2] = -x2_dot + x3; + f[X3] = -x3_dot + x3_rhs; + f[X4] = -x4_dot; + f[V4] = -v4 + x1 + A5_ * x2 + A6_ * x3; + break; + } + case 4: + { + const ScalarT x4_rhs = + (u - x1 - a1_ * x2 - a2_ * x3 - a3_ * x4) / a4_; + f[X1] = -x1_dot + x2; + f[X2] = -x2_dot + x3; + f[X3] = -x3_dot + x4; + f[X4] = -x4_dot + x4_rhs; + f[V4] = -v4 + x1 + A5_ * x2 + A6_ * x3; + break; + } + } + + const ScalarT x5_rhs = (v4 - x5) / T2_; + const ScalarT x6_rhs = (v5 - x6) / T4_; + const ScalarT x7_rhs = (v6 - x7) / T6_; + + f[X5] = -x5_dot + x5_rhs; + f[X6] = -x6_dot + x6_rhs; + f[X7] = -x7_dot + x7_rhs; + f[V5] = -v5 + x5 + T1_ * x5_rhs; + f[V6] = -v6 + x6 + T3_ * x6_rhs; + f[V7] = -v7 + Ks_ * T5_ * x7_rhs; + f[VSS] = -vss + Math::clamp(v7, Lsmin_, Lsmax_); return 0; } @@ -324,16 +367,16 @@ namespace GridKit template int Ieeest::evaluateResidual() { - if (signals_.template isAttached()) - { - ws_[0] = signals_.template readExternalVariable(); - ws_indices_[0] = signals_.template readExternalVariableIndex(); - } + const auto U = static_cast(IeeestExternalVariables::U); + + ws_[U] = signals_.template readExternalVariable(); + ws_indices_[U] = signals_.template readExternalVariableIndex(); const auto* y = y_.getData(); const auto* yp = yp_.getData(); auto* f = f_.getData(); evaluateInternalResidual(y, yp, wb_.data(), ws_.data(), f); + f_.setDataUpdated(); return 0; @@ -345,12 +388,145 @@ namespace GridKit return monitor_.get(); } + // + // Private methods + // + + template + void Ieeest::initializeParameters(const ModelDataT& data) + { + using Params = typename ModelDataT::Parameters; + + parameter_error_count_ = 0; + + auto load_real = [&](auto key, RealT& target, const char* name) -> bool + { + if (!data.parameters.contains(key)) + { + return false; + } + + const auto& value = data.parameters.at(key); + RealT parsed_value{}; + if (const auto* real_value = std::get_if(&value)) + { + parsed_value = *real_value; + } + else if (const auto* index_value = std::get_if(&value)) + { + parsed_value = static_cast(*index_value); + } + else + { + Log::error() << "Ieeest: parameter '" << name << "' must be numeric\n"; + ++parameter_error_count_; + return false; + } + + if (!std::isfinite(parsed_value)) + { + Log::error() << "Ieeest: parameter '" << name << "' must be finite\n"; + ++parameter_error_count_; + return false; + } + + target = parsed_value; + return true; + }; + + load_real(Params::A1, A1_, "A1"); + load_real(Params::A2, A2_, "A2"); + load_real(Params::A3, A3_, "A3"); + load_real(Params::A4, A4_, "A4"); + load_real(Params::A5, A5_, "A5"); + load_real(Params::A6, A6_, "A6"); + load_real(Params::T1, T1_, "T1"); + load_real(Params::T2, T2_, "T2"); + load_real(Params::T3, T3_, "T3"); + load_real(Params::T4, T4_, "T4"); + load_real(Params::T5, T5_, "T5"); + load_real(Params::T6, T6_, "T6"); + load_real(Params::Ks, Ks_, "Ks"); + load_real(Params::Lsmin, Lsmin_, "Lsmin"); + load_real(Params::Lsmax, Lsmax_, "Lsmax"); + if (load_real(Params::Vcl, Vcl_, "Vcl") && Vcl_ != ZERO) + { + Log::warning() << "Ieeest: nonzero Vcl requests a lower input cutout, " + "but input cutout is not implemented and Vcl is ignored\n"; + } + if (load_real(Params::Vcu, Vcu_, "Vcu") && Vcu_ != ZERO) + { + Log::warning() << "Ieeest: nonzero Vcu requests an upper input cutout, " + "but input cutout is not implemented and Vcu is ignored\n"; + } + if (load_real(Params::Tdelay, Tdelay_, "Tdelay") && Tdelay_ != ZERO) + { + Log::warning() << "Ieeest: nonzero Tdelay requests an input delay, " + "but input delay is not implemented and Tdelay is ignored\n"; + } + + setDerivedParameters(); + } + template void Ieeest::initializeMonitor() { using Variable = typename ModelDataT::MonitorableVariables; + + constexpr auto VSS = static_cast(IeeestInternalVariables::VSS); + monitor_->set(Variable::vss, [this] - { return y_.getData()[11]; }); + { return y_.getData()[VSS]; }); + } + + template + void Ieeest::setDerivedParameters() + { + auto check_non_negative = [&](RealT value, const char* name) + { + if (value < ZERO) + { + Log::error() << "Ieeest: " << name << " must be non-negative\n"; + ++parameter_error_count_; + } + }; + + check_non_negative(T2_, "T2"); + check_non_negative(T4_, "T4"); + check_non_negative(T6_, "T6"); + + if (T2_ < TIME_CONSTANT_MINIMUM + || T4_ < TIME_CONSTANT_MINIMUM + || T6_ < TIME_CONSTANT_MINIMUM) + { + Log::warning() << "Ieeest: T2, T4, and T6 below " + << TIME_CONSTANT_MINIMUM + << " s are raised to preserve Hessenberg form\n"; + } + + T2_ = std::max(T2_, TIME_CONSTANT_MINIMUM); + T4_ = std::max(T4_, TIME_CONSTANT_MINIMUM); + T6_ = std::max(T6_, TIME_CONSTANT_MINIMUM); + + a1_ = A1_ + A3_; + a2_ = A2_ + A4_ + A1_ * A3_; + a3_ = A1_ * A4_ + A2_ * A3_; + a4_ = A2_ * A4_; + + auto factor_order = [](RealT linear, RealT quadratic) -> IdxT + { + if (quadratic != ZERO) + { + return static_cast(2); + } + if (linear != ZERO) + { + return static_cast(1); + } + return static_cast(0); + }; + + order_ = factor_order(A1_, A2_) + factor_order(A3_, A4_); } } // namespace Stabilizer diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md index e5765ba1a..d5280ead6 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md @@ -1,94 +1,192 @@ # **IEEE Stabilizer Model (IEEEST)** -Standard IEEE power system stabilizer: 4th-order notch filter, two lead–lag -blocks, washout, and output limiter. +IEEEST is a standard IEEE power system stabilizer with a derived-order notch +filter, two lead-lag blocks, washout, and an output limiter. + +## Notes + +> [!NOTE] +> $V_{\mathrm{cl}}$, $V_{\mathrm{cu}}$, and $T_{\mathrm{delay}}$ are accepted for +> input-format compatibility but are not modeled. A finite nonzero value logs a +> warning and is ignored. + +- The notch-filter order is derived from the two configured denominator + factors. It is not an input parameter. +- The numerator must be proper: its order cannot exceed the derived denominator + order. ## Block Diagram -![](../../../../../docs/Figures/stabilizer_ieeest_diagram.png) +![](../../../../../docs/Figures/PhasorDynamics/IEEEST/diagram.png) -Figure 1: Stabilizer IEEEST model. Figure courtesy of [PowerWorld](https://www.powerworld.com/WebHelp/) +Figure 1: IEEEST block diagram. Figure courtesy of [PowerWorld](https://www.powerworld.com/WebHelp/) ## Model Parameters -Symbol | Units | Description | Typical Value -------------|--------|--------------------------------------|-------------- -$A_1$ | [s] | Notch denominator coefficient | 1.013 -$A_2$ | [s²] | Notch denominator coefficient | 0.013 -$A_3$ | [s] | Notch denominator coefficient | 0.0 -$A_4$ | [s²] | Notch denominator coefficient | 0.0 -$A_5$ | [s] | Notch numerator coefficient | 1.013 -$A_6$ | [s²] | Notch 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 - -The IEEE 421.5 IEEEST also defines a cutout window ($V_{cl}$, $V_{cu}$) and an -input delay ($T_{delay}$). These parameters are accepted for input-format -compatibility but are not modeled here. - -### Derived Parameters +Symbol | Units | JSON | Description | Typical Value | Note +------------------------|----------|----------|--------------------------------------|---------------|------ +$A_1$ | [sec] | `A1` | Notch denominator coefficient | 1.013 | +$A_2$ | [sec²] | `A2` | Notch denominator coefficient | 0.013 | +$A_3$ | [sec] | `A3` | Notch denominator coefficient | 0.0 | +$A_4$ | [sec²] | `A4` | Notch denominator coefficient | 0.0 | +$A_5$ | [sec] | `A5` | Notch numerator coefficient | 1.013 | +$A_6$ | [sec²] | `A6` | Notch numerator coefficient | 0.113 | +$T_1$ | [sec] | `T1` | Lead-lag 1 numerator time constant | 0.0 | +$T_2$ | [sec] | `T2` | Lead-lag 1 denominator time constant | 0.02 | +$T_3$ | [sec] | `T3` | Lead-lag 2 numerator time constant | 0.0 | +$T_4$ | [sec] | `T4` | Lead-lag 2 denominator time constant | 0.0 | +$T_5$ | [sec] | `T5` | Washout numerator time constant | 1.65 | +$T_6$ | [sec] | `T6` | Washout denominator time constant | 1.65 | +$K_s$ | [p.u.] | `Ks` | Stabilizer gain | 3.0 | +$L_s^{\min}$ | [p.u.] | `Lsmin` | Minimum stabilizer output limit | -0.1 | +$L_s^{\max}$ | [p.u.] | `Lsmax` | Maximum stabilizer output limit | 0.1 | +$V_{\mathrm{cl}}$ | [p.u.] | `Vcl` | Lower input cutout threshold | 0.0 | Nonzero values warn and are ignored +$V_{\mathrm{cu}}$ | [p.u.] | `Vcu` | Upper input cutout threshold | 0.0 | Nonzero values warn and are ignored +$T_{\mathrm{delay}}$ | [sec] | `Tdelay` | Input delay | 0.0 | Nonzero values warn and are ignored + +### Parameter Validation + +Let $m$ be the numerator order. Exact zero selects the notch-filter topology: + +```math +\begin{aligned} + m &= + \begin{cases} + 2 & A_6 \ne 0 \\ + 1 & A_6 = 0,\ A_5 \ne 0 \\ + 0 & A_5 = A_6 = 0 + \end{cases} \\ + m &\le n \\ + T_2,T_4,T_6 &\ge 0 \\ + L_s^{\min} &< L_s^{\max} \\ + a_n &\ne 0 \quad n>0 +\end{aligned} +``` + +The raw and derived coefficients must be finite. The final condition applies to +the leading coefficient selected by the derived order. Negative denominator +time constants are rejected. + +### Model Derived Parameters + +For one denominator factor, define its exact order by + +```math +d(p,q)= +\begin{cases} + 2 & q \ne 0 \\ + 1 & q = 0,\ p \ne 0 \\ + 0 & p = q = 0. +\end{cases} +``` + +The model derives the notch-filter order and expanded denominator coefficients +once during construction: ```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 + n &= d(A_1,A_2) + d(A_3,A_4) \\ + 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} ``` +Let $\epsilon_T=10^{-3}\ \mathrm{s}$. A denominator time constant below +$\epsilon_T$ is raised to that floor and logs a warning: + +```math +T \leftarrow \max(T,\epsilon_T), +\qquad T\in\{T_2,T_4,T_6\}. +``` + +## Model Ports + +Name | Port | Init | Description +---------|--------|-------|------ +`input` | Input | Known | Stabilizer input signal +`output` | Output | Known | Stabilizer output signal + +The input must be attached to a linked signal. Assigning the output signal is +optional. + ## Model Variables ### Internal Variables #### Differential -Symbol | Units | Description -----------------------|--------|------------ -$x_1, x_2, x_3, x_4$ | [-] | Notch filter states -$x_5$ | [-] | Lead–lag 1 state -$x_6$ | [-] | Lead–lag 2 state -$x_7$ | [-] | Washout state +All seven differential variables are always present. Inactive notch-filter +states remain differential and are frozen by $0=-\dot{x}_i$. + +Symbol | Units | Description | Note +-------|-------------|---------------------------------------|------ +$x_1$ | [p.u.] | Notch-filter signal state | Active for $n\ge1$; frozen for $n=0$ +$x_2$ | [p.u./sec] | First derivative of filtered signal | Active for $n\ge2$; otherwise frozen +$x_3$ | [p.u./sec²] | Second derivative of filtered signal | Active for $n\ge3$; otherwise frozen +$x_4$ | [p.u./sec³] | Third derivative of filtered signal | Active for $n=4$; otherwise frozen +$x_5$ | [p.u.] | Lead-lag 1 state | State 5 in Fig. 1 +$x_6$ | [p.u.] | Lead-lag 2 state | State 6 in Fig. 1 +$x_7$ | [p.u.] | Washout state | State 7 in Fig. 1 #### Algebraic -Symbol | Units | Description ------------|--------|------------ -$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 (model output) +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_{\mathrm{ss}}$ | [p.u.] | Limited stabilizer signal (model output) | ### External Variables +#### Differential + +None. + #### Algebraic -Symbol | Units | Description --------|--------|------------ -$u$ | [p.u.] | Stabilizer input signal +Symbol | Units | Type | Description | Note +-------|--------|-------|-------------------------|------ +$u$ | [p.u.] | Known | Stabilizer input signal | Required signal port `input` ## Model Equations ### Differential Equations +A runtime switch selects one notch-filter realization using the derived order +$n$. The state layout does not change with $n$. + ```math \begin{aligned} -0 &= -\dot{x}_1 + x_2 \\ -0 &= -\dot{x}_2 + x_3 \\ -0 &= -\dot{x}_3 + x_4 \\ -0 &= -\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 \\ -0 &= -T_2 \dot{x}_5 - x_5 + v_4 \\ -0 &= -T_4 \dot{x}_6 - x_6 + v_5 \\ -0 &= -T_6 \dot{x}_7 - x_7 + v_6 + 0 &= -\dot{x}_1 + + \begin{cases} + 0 & n=0 \\ + \dfrac{-x_1+u}{a_1} & n=1 \\ + x_2 & n\in\{2,3,4\} + \end{cases} \\ + 0 &= -\dot{x}_2 + + \begin{cases} + 0 & n\in\{0,1\} \\ + \dfrac{-x_1-a_1x_2+u}{a_2} & n=2 \\ + x_3 & n\in\{3,4\} + \end{cases} \\ + 0 &= -\dot{x}_3 + + \begin{cases} + 0 & n\in\{0,1,2\} \\ + \dfrac{-x_1-a_1x_2-a_2x_3+u}{a_3} & n=3 \\ + x_4 & n=4 + \end{cases} \\ + 0 &= -\dot{x}_4 + + \begin{cases} + 0 & n\in\{0,1,2,3\} \\ + \dfrac{-x_1-a_1x_2-a_2x_3-a_3x_4+u}{a_4} & n=4 + \end{cases} \\ + 0 &= -\dot{x}_5 + \dfrac{v_4-x_5}{T_2} \\ + 0 &= -\dot{x}_6 + \dfrac{v_5-x_6}{T_4} \\ + 0 &= -\dot{x}_7 + \dfrac{v_6-x_7}{T_6}. \end{aligned} ``` @@ -96,19 +194,58 @@ $u$ | [p.u.] | Stabilizer input signal ```math \begin{aligned} -0 &= -v_4 + x_1 + A_5 x_2 + A_6 x_3 \\ -0 &= -T_2(v_5 - x_5) + T_1(v_4 - x_5) \\ -0 &= -T_4(v_6 - x_6) + T_3(v_5 - x_6) \\ -0 &= -T_6 v_7 + K_s T_5(v_6 - x_7) \\ -0 &= -V_{ss} + \text{clamp}(v_7, L_s^{\min}, L_s^{\max}) + 0 &= -v_4 + + \begin{cases} + u & n=0 \\ + x_1 + \dfrac{A_5}{a_1}(-x_1+u) & n=1 \\ + x_1 + A_5x_2 + + \dfrac{A_6}{a_2}(-x_1-a_1x_2+u) & n=2 \\ + x_1 + A_5x_2 + A_6x_3 & n\in\{3,4\} + \end{cases} \\ + 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_{\mathrm{ss}} + + \operatorname{clamp}(v_7,L_s^{\min},L_s^{\max}). \end{aligned} ``` -The output limiter uses GridKit's smooth -[Clamp](../../../../CommonMath.md#derived-functions). +CommonMath defines the smooth +[clamp](../../../../CommonMath.md#clamp). ## Initialization -All states and their derivatives initialize to zero. The stabilizer comes -online at rest and produces signal only in response to deviations in the input -$u$. +### Input Initialization + +```math +u \leftarrow \text{stabilizer input signal}. +``` + +### Internal Initialization + +```math +\begin{aligned} + x_1 &\leftarrow + \begin{cases} + 0 & n=0 \\ + u & n\in\{1,2,3,4\} + \end{cases} \\ + x_2,x_3,x_4 &\leftarrow 0 \\ + v_4,x_5,v_5,x_6,v_6,x_7 &\leftarrow u \\ + v_7 &\leftarrow 0 \\ + V_{\mathrm{ss}} &\leftarrow + \operatorname{clamp}(v_7,L_s^{\min},L_s^{\max}) \\ + \dot{x}_i &\leftarrow 0, + \qquad i\in\{1,\ldots,7\}. +\end{aligned} +``` + +### Output Initialization + +None. + +## Monitorable Outputs + +Output | Units | Description | Note +-------|--------|---------------------------------------------|------ +`vss` | [p.u.] | Limited stabilizer signal $V_{\mathrm{ss}}$ | Model output diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index 1c60b16b0..3a51f1d46 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -439,6 +439,29 @@ namespace GridKit addComponent(hygov); } + // Add IEEEST stabilizers before exciters that read their output during + // initialization. + for (const auto& stabdata : data.stabilizer) + { + auto* stabilizer = new Ieeest(stabdata); + + if (stabdata.signal_inputs.contains(IeeestSignalInputs::input)) + { + IdxT input = stabdata.signal_inputs.at(IeeestSignalInputs::input); + constexpr auto U = IeeestExternalVariables::U; + stabilizer->getSignals().template attachSignalNode(getSignal(input)); + } + + if (stabdata.signal_outputs.contains(IeeestSignalOutputs::output)) + { + IdxT output = stabdata.signal_outputs.at(IeeestSignalOutputs::output); + constexpr auto VSS = IeeestInternalVariables::VSS; + stabilizer->getSignals().template assignSignalNode(getSignal(output)); + } + + addComponent(stabilizer); + } + for (const auto& excitedata : data.exciter) { IdxT bus_index = 0; @@ -548,28 +571,6 @@ namespace GridKit addComponent(exciter); } - // Add IEEEST stabilizers - for (const auto& stabdata : data.stabilizer) - { - auto* stabilizer = new Ieeest(stabdata); - - if (stabdata.signal_inputs.contains(IeeestSignalInputs::input)) - { - IdxT input = stabdata.signal_inputs.at(IeeestSignalInputs::input); - constexpr auto U = IeeestExternalVariables::U; - stabilizer->getSignals().template attachSignalNode(getSignal(input)); - } - - if (stabdata.signal_outputs.contains(IeeestSignalOutputs::output)) - { - IdxT output = stabdata.signal_outputs.at(IeeestSignalOutputs::output); - constexpr auto VSS = IeeestInternalVariables::VSS; - stabilizer->getSignals().template assignSignalNode(getSignal(output)); - } - - addComponent(stabilizer); - } - // Add REPCA plant controllers after the signal producers they read at // initialization for (const auto& repcadata : data.repca) diff --git a/docs/Figures/stabilizer_ieeest_diagram.png b/docs/Figures/PhasorDynamics/IEEEST/diagram.png similarity index 100% rename from docs/Figures/stabilizer_ieeest_diagram.png rename to docs/Figures/PhasorDynamics/IEEEST/diagram.png diff --git a/tests/UnitTests/PhasorDynamics/ComponentConnectionTests.hpp b/tests/UnitTests/PhasorDynamics/ComponentConnectionTests.hpp index bcb193ed7..6efc266d4 100644 --- a/tests/UnitTests/PhasorDynamics/ComponentConnectionTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ComponentConnectionTests.hpp @@ -12,10 +12,13 @@ #include #include #include +#include #include #include #include #include +#include +#include #include #include #include @@ -329,6 +332,68 @@ namespace GridKit return success.report(__func__); } + + /// IEEEST reinitializes before IEEET1 so the exciter does not latch a + /// stale stabilizer output. + TestOutcome ieeestIeeet1() + { + using namespace PhasorDynamics; + using namespace PhasorDynamics::Exciter; + using namespace PhasorDynamics::Stabilizer; + + TestStatus success = true; + + constexpr IdxT bus_id = static_cast(1); + constexpr IdxT input_id = static_cast(20); + constexpr IdxT vs_id = static_cast(21); + SystemModelData data; + + auto& bus = data.bus.emplace_back(); + bus.bus_id = bus_id; + bus.bus_type = BusData::BusType::SLACK; + bus.Vr0 = ONE; + bus.Vi0 = ZERO; + + data.signal = {{"IEEEST input", input_id}, {"IEEEST output", vs_id}}; + + auto& source = data.constant_source.emplace_back(); + source.signal_outputs[ConstantSignalSourceSignalOutputs::sr] = input_id; + source.parameters[ConstantSignalSourceParameters::Sr] = ZERO; + + auto& stabilizer = data.stabilizer.emplace_back(); + stabilizer.signal_inputs[IeeestSignalInputs::input] = input_id; + stabilizer.signal_outputs[IeeestSignalOutputs::output] = vs_id; + + auto& exciter = data.exciter.emplace_back(); + exciter.buses[Ieeet1Buses::bus] = bus_id; + exciter.signal_inputs[Ieeet1SignalInputs::vs] = vs_id; + exciter.parameters[Ieeet1Parameters::Tr] = static_cast(0.02); + + SystemModel system(data); + success *= system.allocate() == 0; + + auto* vs = system.getSignal(vs_id); + if (vs == nullptr || !vs->linked()) + { + success = false; + return success.report(__func__); + } + + // Reinitialization must not let the exciter latch this stale value. + vs->init(static_cast(0.05)); + success *= system.initialize() == 0; + success *= system.evaluateResidual() == 0; + success *= isEqual(vs->read(), ZERO, kTol); + + const auto& residual = system.getResidual(); + const auto* residual_data = residual.getData(); + for (IdxT row = 0; row < residual.getSize(); ++row) + { + success *= isEqual(residual_data[row], ZERO, kTol); + } + + return success.report(__func__); + } }; /// Production SystemModel wiring tests for GASTPTI and synchronous machines. diff --git a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp index 14cd9cd88..4d0e73ead 100644 --- a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp +++ b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp @@ -1,433 +1,918 @@ #pragma once +#include +#include +#include +#include #include #include #include +#include +#include +#include +#include #include #include #include #include #include -#include +#include #include +#include +#include #include namespace GridKit { namespace Testing { - template + using Log = ::GridKit::Utilities::Logger; + + template class StabilizerIeeestTests { public: - using RealT = typename PhasorDynamics::Component::RealT; + using ScalarT = scalar_type; + using IdxT = index_type; + using RealT = typename PhasorDynamics::Component::RealT; StabilizerIeeestTests() = default; ~StabilizerIeeestTests() = default; - TestOutcome constructor() + TestOutcome validation() { TestStatus success = true; - auto data = makeTestData(); - auto* stab = new PhasorDynamics::Stabilizer::Ieeest(data); + IeeestT empty; + success *= (empty.size() == static_cast(Internal::MAXIMUM)); + success *= (empty.getMonitor() == nullptr); + success *= (static_cast(Internal::MAXIMUM) == 12); + success *= (static_cast(Internal::VSS) == 11); + + Fixture configured(makeOrderData(4)); + success *= configured.prepare(); + success *= (configured.model.getMonitor() != nullptr); + + for (size_t row = 0; row < configured.model.tag().size(); ++row) + { + success *= !configured.model.tag()[row]; + } + success *= (configured.model.tagDifferentiable() == 0); + for (size_t row = 0; row < configured.model.tag().size(); ++row) + { + const bool expected = row <= static_cast(Internal::X7); + if (configured.model.tag()[row] != expected) + { + std::cout << "IEEEST differentiability tag " << row << " mismatch\n"; + success = false; + } + } + + noteExpectedLogs("Testing IEEEST invalid configurations. Logged errors " + "and time-constant or unsupported-feature warnings are expected."); + + for (size_t order = 0; order <= 4; ++order) + { + success *= verifies(makeOrderData(order)); + } + success *= verifies(makeOrderThreeReverseData()); + success *= verifies(makeSymmetricOrderFourData()); + + // The required input must be both attached and linked. + { + IeeestT missing(makeOrderData(4)); + success *= (missing.allocate() == 0); + success *= (missing.verify() > 0); + } + { + PhasorDynamics::SignalNode unlinked_node; + IeeestT unlinked(makeOrderData(4)); + unlinked.getSignals().template attachSignalNode(&unlinked_node); + success *= (unlinked.allocate() == 0); + success *= (unlinked.verify() > 0); + } + + // The VSS output assignment remains optional. + Fixture no_output(makeOrderData(4), false); + success *= no_output.prepare(); + + // Integer JSON values are accepted for real parameters; booleans are not. + { + auto data = makeOrderData(4); + data.parameters[Params::Ks] = static_cast(3); + success *= verifies(data); + data.parameters[Params::Ks] = true; + success *= !verifies(data); + } + + const RealT nan = std::numeric_limits::quiet_NaN(); + const RealT infinity = std::numeric_limits::infinity(); + const std::array nonfinite_values{{nan, infinity, -infinity}}; + const std::array real_parameters{{ + Params::A1, + Params::A2, + Params::A3, + Params::A4, + Params::A5, + Params::A6, + Params::T1, + Params::T2, + Params::T3, + Params::T4, + Params::T5, + Params::T6, + Params::Ks, + Params::Lsmin, + Params::Lsmax, + Params::Vcl, + Params::Vcu, + Params::Tdelay, + }}; + + for (const Params parameter : real_parameters) + { + for (const RealT value : nonfinite_values) + { + auto data = makeOrderData(4); + data.parameters[parameter] = value; + success *= !verifies(data); + } + } + + success *= invalidParameterLeavesDefault(Params::Ks, + true, + "boolean parameter fallback"); + success *= invalidParameterLeavesDefault(Params::T1, + nan, + "nonfinite parameter fallback"); + + // Expanded coefficients must remain finite. + { + auto data = makeOrderData(2); + data.parameters[Params::A1] = std::numeric_limits::max(); + data.parameters[Params::A3] = std::numeric_limits::max(); + success *= !verifies(data); + } + + // Finite quadratic factors select order four even when their leading + // product overflows. + { + auto data = makeOrderData(0); + data.parameters[Params::A2] = std::numeric_limits::max(); + data.parameters[Params::A4] = static_cast(2.0); + success *= !verifies(data); + } + + // Exact nonzero factors still select order four; an underflowed a4 is invalid. + { + auto data = makeOrderData(0); + data.parameters[Params::A2] = std::numeric_limits::min(); + data.parameters[Params::A4] = std::numeric_limits::min(); + success *= !verifies(data); + } + + // The numerator may not have a higher order than the denominator. + { + auto data = makeOrderData(0); + data.parameters[Params::A5] = static_cast(0.1); + success *= !verifies(data); + data = makeOrderData(0); + data.parameters[Params::A6] = static_cast(0.1); + success *= !verifies(data); + data = makeOrderData(1); + data.parameters[Params::A6] = static_cast(0.1); + success *= !verifies(data); + success *= verifies(makeOrderData(2)); + } + + for (const Params parameter : {Params::T2, Params::T4, Params::T6}) + { + auto data = makeOrderData(4); + data.parameters[parameter] = static_cast(-0.1); + success *= !verifies(data); + } - success *= (stab != nullptr); - success *= (stab->getMonitor() != nullptr); + { + auto data = makeOrderData(4); + data.parameters[Params::Lsmin] = static_cast(0.1); + data.parameters[Params::Lsmax] = static_cast(0.1); + success *= !verifies(data); + data.parameters[Params::Lsmin] = static_cast(0.2); + success *= !verifies(data); + } + + // Finite unsupported inputs warn but do not invalidate the model. + { + auto data = makeOrderData(4); + data.parameters[Params::Vcl] = static_cast(-0.5); + data.parameters[Params::Vcu] = static_cast(0.5); + data.parameters[Params::Tdelay] = static_cast(0.1); + success *= verifies(data); + } - delete stab; + // Topology uses exact coefficient zeros, including signed zero. + success *= checkTopology(makeOrderData(0), 0, "exact zero"); + { + auto data = makeOrderData(0); + data.parameters[Params::A1] = -ZERO; + data.parameters[Params::A2] = -ZERO; + data.parameters[Params::A3] = -ZERO; + data.parameters[Params::A4] = -ZERO; + success *= checkTopology(data, 0, "signed zero"); + } + { + auto data = makeOrderData(1); + data.parameters[Params::A1] = static_cast(-0.4); + success *= checkTopology(data, 1, "negative coefficient"); + } + { + auto data = makeOrderData(1); + data.parameters[Params::A1] = static_cast(1.0e-12); + success *= checkTopology(data, 1, "small nonzero coefficient"); + } return success.report(__func__); } - /** - * @brief All states initialize to zero (stabilizer at rest). - * With u = 0, all residuals should be zero. - */ - TestOutcome zeroInitialResidual() + TestOutcome initializationAndSignals() { TestStatus success = true; - // Create signal nodes for input (u) and output (Vss) - PhasorDynamics::SignalNode u_node; - PhasorDynamics::SignalNode vss_node; - ScalarT u_value{0.0}; - IdxT u_index = 12; // beyond internal variables - ScalarT vss_value{0.0}; - IdxT vss_index = INVALID_INDEX; + for (size_t order = 0; order <= 4; ++order) + { + Fixture fixture(makeOrderData(order)); + fixture.u() = static_cast(0.25); + success *= fixture.initialize(); + success *= fixture.outputLinked(); + success *= (fixture.outputIndex() == static_cast(Internal::VSS)); + + const auto* y = fixture.model.y().getData(); + const auto* yp = fixture.model.yp().getData(); + for (size_t row = 0; row < static_cast(Internal::MAXIMUM); ++row) + { + RealT expected = 0.0; + if ((order > 0 && row == static_cast(Internal::X1)) + || (row >= static_cast(Internal::X5) + && row <= static_cast(Internal::V6))) + { + expected = 0.25; + } + success *= rowMatches(static_cast(y[row]), expected, "state", row, "initialization"); + success *= rowMatches(static_cast(yp[row]), 0.0, "derivative", row, "initialization"); + } - // Link signal nodes to backing storage - u_node.set(&u_value, &u_index); - vss_node.set(&vss_value, &vss_index); + success *= (fixture.evaluate() == 0); + success *= allResidualsZero(fixture.model); + success *= scalarMatches(static_cast(fixture.output()), + static_cast(y[static_cast(Internal::VSS)]), + "assigned VSS output"); + } - auto data = makeTestData(); - PhasorDynamics::Stabilizer::Ieeest stab(data); + // Zero denominator time constants use the documented floor and retain + // a consistent initial condition. + { + auto data = makeOrderData(4); + data.parameters[Params::T2] = ZERO; + data.parameters[Params::T4] = ZERO; + data.parameters[Params::T6] = ZERO; + Fixture floors(data); + floors.u() = static_cast(0.25); + success *= floors.initialize(); + success *= (floors.evaluate() == 0); + success *= allResidualsZero(floors.model); + } + { + auto data = makeOrderData(4); + data.parameters[Params::T2] = static_cast(5.0e-4); + data.parameters[Params::T4] = static_cast(5.0e-4); + data.parameters[Params::T6] = static_cast(5.0e-4); + Fixture floors(data); + floors.u() = static_cast(0.25); + success *= floors.initialize(); + success *= (floors.evaluate() == 0); + success *= allResidualsZero(floors.model); + } - // Wire: stabilizer reads u_node as input, writes vss_node as output - stab.getSignals().template attachSignalNode(&u_node); - stab.getSignals().template assignSignalNode(&vss_node); + // Initialization applies the output limit to the zero unlimited signal. + { + auto data = makeOrderData(4); + data.parameters[Params::Lsmin] = static_cast(0.2); + data.parameters[Params::Lsmax] = static_cast(0.6); + Fixture limited(data); + limited.u() = static_cast(0.25); + success *= limited.initialize(); + success *= scalarMatches(static_cast(limited.output()), + static_cast(0.2), + "limited initial VSS", + kClampTol); + } - stab.allocate(); - success *= (stab.verify() == 0); - stab.initialize(); - stab.evaluateResidual(); + // A nonfinite input rejects initialization without changing state. + for (const RealT input : {std::numeric_limits::quiet_NaN(), + std::numeric_limits::infinity(), + -std::numeric_limits::infinity()}) + { + Fixture fixture(makeOrderData(4)); + success *= fixture.prepare(); + poisonState(fixture.model); + const auto y_before = copyVector(fixture.model.y()); + const auto yp_before = copyVector(fixture.model.yp()); + fixture.u() = static_cast(input); + success *= (fixture.model.initialize() != 0); + success *= vectorMatches(fixture.model.y(), y_before, "rejected state"); + success *= vectorMatches(fixture.model.yp(), yp_before, "rejected derivative"); + } - auto tol = 10 * std::numeric_limits::epsilon(); - const auto& f = stab.getResidual(); - const auto* f_data = f.getData(); - for (size_t i = 0; i < f.getSize(); ++i) + // Verification failure is also detected before any state mutation. { - if (!isEqual(f_data[i], 0.0, tol)) - { - std::cout << "Non-zero residual at index " << i << ": " << f_data[i] << "\n"; - success = false; - } + IeeestT model(makeOrderData(4)); + success *= (model.allocate() == 0); + poisonState(model); + const auto y_before = copyVector(model.y()); + const auto yp_before = copyVector(model.yp()); + success *= (model.initialize() != 0); + success *= vectorMatches(model.y(), y_before, "unverified state"); + success *= vectorMatches(model.yp(), yp_before, "unverified derivative"); } - // Verify output signal is linked and reads the correct value - success *= vss_node.linked(); - success *= (vss_node.getVariableIndex() == 11); - success *= isEqual(vss_node.read(), static_cast(0.0), tol); + return success.report(__func__); + } + + TestOutcome residualEquations() + { + TestStatus success = true; + + const std::array order_zero{{-0.01, -0.02, -0.03, -0.04, 0.25, 0.24, -0.01, -0.30, -0.25, -0.31, 1.15, 0.0}}; + const std::array order_one{{0.99, -0.02, -0.03, -0.04, 0.25, 0.24, -0.01, -0.20, -0.25, -0.31, 1.15, 0.0}}; + const std::array order_two{{0.19, 1.88, -0.03, -0.04, 0.25, 0.24, -0.01, 0.54, -0.25, -0.31, 1.15, 0.0}}; + const std::array order_three{{0.19, 0.28, 4.153333333333333, -0.04, 0.25, 0.24, -0.01, -0.42, -0.25, -0.31, 1.15, 0.0}}; + const std::array order_three_reverse{{0.19, 0.28, 4.745, -0.04, 0.25, 0.24, -0.01, -0.42, -0.25, -0.31, 1.15, 0.0}}; + const std::array order_four{{0.19, 0.28, 0.37, 1.0975, 0.25, 0.24, -0.01, -0.42, -0.25, -0.31, 1.15, 0.0}}; + const std::array symmetric_four{{0.19, 0.28, 0.37, 2.71, 0.25, 0.24, -0.01, -0.42, -0.25, -0.31, 1.15, 0.0}}; + + success *= checkResidual(makeOrderData(0), order_zero, "order 0"); + success *= checkResidual(makeOrderData(1), order_one, "order 1"); + success *= checkResidual(makeOrderData(2), order_two, "order 2"); + success *= checkResidual(makeOrderData(3), order_three, "order 3 first quadratic"); + success *= checkResidual(makeOrderThreeReverseData(), order_three_reverse, "order 3 second quadratic"); + success *= checkResidual(makeOrderData(4), order_four, "order 4"); + success *= checkResidual(makeSymmetricOrderFourData(), symmetric_four, "symmetric order 4"); 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() + TestOutcome monitor() { TestStatus success = true; - PhasorDynamics::SignalNode u_node; - PhasorDynamics::SignalNode vss_node; - ScalarT u_value{0.5}; - IdxT u_index = 12; - ScalarT vss_value{0.0}; - IdxT vss_index = INVALID_INDEX; - - u_node.set(&u_value, &u_index); - vss_node.set(&vss_value, &vss_index); - - auto data = makeTestData(); - PhasorDynamics::Stabilizer::Ieeest stab(data); - - stab.getSignals().template attachSignalNode(&u_node); - stab.getSignals().template assignSignalNode(&vss_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]: -T2*x5_dot - x5 + v4 - 0.24, // f[5]: -T4*x6_dot - x6 + v5 - -0.05, // f[6]: -T6*x7_dot - x7 + v6 - -0.42, // f[7]: -v4 + x1 + A5*x2 + A6*x3 - -0.25, // f[8]: -T2*(v5 - x5) + T1*(v4 - x5) - -0.31, // f[9]: -T4*(v6 - x6) + T3*(v5 - x6) - 5.75, // f[10]: -T6*v7 + Ks*T5*(v6 - x7) - 0.0, // f[11]: limiter (v7=0.05 within [-0.1, 0.1]) - }; - - // Looser tolerance for f[11] — Math::clamp is a smooth ramp approximation. - const auto loose_tol = static_cast(1.0e-4); - auto& residual = stab.getResidual(); - const auto* residual_data = residual.getData(); - - 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_data[i], res_answer[i], test_tol)) - { - std::cout << "Incorrect result for residual " << i << ": " - << std::setprecision(15) << residual_data[i] - << " != " << res_answer[i] << "\n"; - success = false; - } + Fixture fixture(makeOrderData(4)); + success *= fixture.prepare(); + + auto* y = fixture.model.y().getData(); + y[static_cast(Internal::VSS)] = static_cast(0.075); + fixture.model.y().setDataUpdated(); + + RealT time = 0.0; + Model::VariableMonitorController controller(time); + controller.addMonitor(fixture.model.getMonitor()); + + std::stringstream os; + controller.addSink({Model::VariableMonitorFormat::CSV}, os); + controller.start(); + controller.print(); + controller.stop(); + + std::string header; + std::string values; + std::getline(os, header); + std::getline(os, values); + success *= (header == "t,Ieeest_ieeest_test_vss"); + const auto monitored = Tokenizer(values, ',')(); + success *= (monitored.size() == 2); + if (monitored.size() == 2) + { + success *= scalarMatches(monitored[1], static_cast(0.075), "monitored VSS"); } - // Verify output signal reads the stabilizer output - success *= isEqual(vss_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> cases; + for (size_t order = 0; order <= 4; ++order) + { + cases.emplace_back(makeOrderData(order), order); + } + cases.emplace_back(makeOrderThreeReverseData(), 3); + cases.emplace_back(makeSymmetricOrderFourData(), 4); - std::vector - dependency_tracking_jacobian = DependencyTrackingJacobian(data); + for (const auto& [data, order] : cases) + { + const auto dependency_jacobian = dependencyTrackingJacobian(data, success); + const auto enzyme_jacobian = enzymeJacobian(data, success); - std::vector - enzyme_jacobian = EnzymeJacobian(data); + success *= (dependency_jacobian.size() == enzyme_jacobian.size()); + const auto rows = std::min(dependency_jacobian.size(), enzyme_jacobian.size()); + for (size_t row = 0; row < rows; ++row) + { + if (!isEqual(dependency_jacobian[row], enzyme_jacobian[row], kTol)) + { + std::cout << "IEEEST Jacobian row " << row << " for order " << order + << " differs between dependency tracking and Enzyme\n"; + printDependencyMap(dependency_jacobian[row], "dependency tracking"); + printDependencyMap(enzyme_jacobian[row], "Enzyme"); + success = false; + } + } - // 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)); + for (size_t row = order; row < 4; ++row) + { + const DependencyMap expected{{row, static_cast(-1.0)}}; + if (!isEqual(dependency_jacobian[row], expected, kTol) + || !isEqual(enzyme_jacobian[row], expected, kTol)) + { + std::cout << "IEEEST inactive notch row " << row << " for order " + << order << " is not the frozen derivative diagonal\n"; + success = false; + } + } } return success.report(__func__); } +#endif private: - std::vector DependencyTrackingJacobian( - PhasorDynamics::Stabilizer::IeeestData ieeestdata) + using Params = PhasorDynamics::Stabilizer::IeeestParameters; + using Internal = PhasorDynamics::Stabilizer::IeeestInternalVariables; + using External = PhasorDynamics::Stabilizer::IeeestExternalVariables; + using Mon = PhasorDynamics::Stabilizer::IeeestMonitorableVariables; + using Data = PhasorDynamics::Stabilizer::IeeestData; + using IeeestT = PhasorDynamics::Stabilizer::Ieeest; + using DependencyMap = DependencyTracking::Variable::DependencyMap; + + static constexpr RealT kTol = + static_cast(100.0) * std::numeric_limits::epsilon(); + static constexpr RealT kClampTol = static_cast(1.0e-4); + + static constexpr std::array kRowNames{{"X1", "X2", "X3", "X4", "X5", "X6", "X7", "V4", "V5", "V6", "V7", "VSS"}}; + + template + class Fixture { - using DepVar = DependencyTracking::Variable; + private: + T u_value_{0}; + IdxT u_index_{static_cast(Internal::MAXIMUM)}; + PhasorDynamics::SignalNode u_node_; + PhasorDynamics::SignalNode vss_node_; + + public: + explicit Fixture(const Data& data, bool assign_output = true) + : model(data) + { + u_node_.set(&u_value_, &u_index_); + model.getSignals().template attachSignalNode(&u_node_); + if (assign_output) + { + model.getSignals().template assignSignalNode(&vss_node_); + } + } + + Fixture(const Fixture&) = delete; + Fixture& operator=(const Fixture&) = delete; + + bool prepare() + { + return model.allocate() == 0 && model.verify() == 0; + } + + bool initialize() + { + return prepare() && model.initialize() == 0; + } + + int evaluate() + { + return model.evaluateResidual(); + } - // Set up signal nodes with DependencyTracking scalar type - PhasorDynamics::SignalNode u_node; - PhasorDynamics::SignalNode vss_node; - DepVar u_value{0.5}; - IdxT u_index = 12; - DepVar vss_value{0.0}; - IdxT vss_index = INVALID_INDEX; + T& u() + { + return u_value_; + } - u_node.set(&u_value, &u_index); - vss_node.set(&vss_value, &vss_index); + IdxT uIndex() const + { + return u_index_; + } - PhasorDynamics::Stabilizer::Ieeest stab(ieeestdata); - stab.getSignals().template attachSignalNode(&u_node); - stab.getSignals().template assignSignalNode(&vss_node); + T output() const + { + return vss_node_.read(); + } - stab.allocate(); - stab.initialize(); + bool outputLinked() const + { + return vss_node_.linked(); + } - // --- d/dy: tag internal variables as independent --- - auto* y = stab.y().getData(); - for (size_t i = 0; i < stab.size(); ++i) + IdxT outputIndex() const { - y[i].setVariableNumber(i); + return vss_node_.getVariableIndex(); } - // Tag external signal u as an additional independent variable - u_value.setVariableNumber(stab.size()); - u_value.setValue(0.5); - setStatePointDep(stab); + PhasorDynamics::Stabilizer::Ieeest model; + }; - stab.evaluateResidual(); - auto& residual_y_view = stab.getResidual(); - std::vector residual_y(residual_y_view.getData(), residual_y_view.getData() + residual_y_view.getSize()); + Data makeBaseData() const + { + Data data; + data.device_class = "Ieeest"; + data.disambiguation_string = "ieeest_test"; + data.monitored_variables.insert(Mon::vss); + + data.parameters[Params::A1] = ZERO; + data.parameters[Params::A2] = ZERO; + data.parameters[Params::A3] = ZERO; + data.parameters[Params::A4] = ZERO; + data.parameters[Params::A5] = ZERO; + data.parameters[Params::A6] = ZERO; + data.parameters[Params::T1] = static_cast(0.5); + data.parameters[Params::T2] = static_cast(1.0); + data.parameters[Params::T3] = static_cast(0.3); + data.parameters[Params::T4] = static_cast(1.0); + data.parameters[Params::T5] = static_cast(2.0); + data.parameters[Params::T6] = static_cast(5.0); + data.parameters[Params::Ks] = static_cast(10.0); + data.parameters[Params::Lsmin] = static_cast(-0.1); + data.parameters[Params::Lsmax] = static_cast(0.1); + data.parameters[Params::Vcl] = ZERO; + data.parameters[Params::Vcu] = ZERO; + data.parameters[Params::Tdelay] = ZERO; + return data; + } - // --- d/dy': tag derivatives as independent --- - stab.initialize(); - auto* yp = stab.yp().getData(); - for (size_t i = 0; i < stab.size(); ++i) + Data makeOrderData(size_t order) const + { + auto data = makeBaseData(); + switch (order) { - yp[i].setVariableNumber(i); + case 0: + break; + case 1: + data.parameters[Params::A1] = static_cast(0.4); + data.parameters[Params::A5] = static_cast(0.5); + break; + case 2: + data.parameters[Params::A1] = static_cast(0.1); + data.parameters[Params::A2] = static_cast(0.2); + data.parameters[Params::A5] = static_cast(0.5); + data.parameters[Params::A6] = static_cast(0.6); + break; + case 3: + data.parameters[Params::A1] = static_cast(0.1); + data.parameters[Params::A2] = static_cast(0.2); + data.parameters[Params::A3] = static_cast(0.3); + data.parameters[Params::A5] = static_cast(0.5); + data.parameters[Params::A6] = static_cast(0.6); + break; + case 4: + data.parameters[Params::A1] = static_cast(0.1); + data.parameters[Params::A2] = static_cast(0.2); + data.parameters[Params::A3] = static_cast(0.3); + data.parameters[Params::A4] = static_cast(0.4); + data.parameters[Params::A5] = static_cast(0.5); + data.parameters[Params::A6] = static_cast(0.6); + break; + default: + break; } + return data; + } + + Data makeOrderThreeReverseData() const + { + auto data = makeBaseData(); + data.parameters[Params::A1] = static_cast(0.1); + data.parameters[Params::A3] = static_cast(0.3); + data.parameters[Params::A4] = static_cast(0.4); + data.parameters[Params::A5] = static_cast(0.5); + data.parameters[Params::A6] = static_cast(0.6); + return data; + } + + Data makeSymmetricOrderFourData() const + { + auto data = makeOrderData(4); + data.parameters[Params::A1] = ZERO; + data.parameters[Params::A3] = ZERO; + return data; + } - u_value = 0.5; - setStatePointDep(stab); + bool verifies(const Data& data) const + { + Fixture fixture(data, false); + return fixture.prepare(); + } - stab.evaluateResidual(); - auto& residual_yp_view = stab.getResidual(); - std::vector residual_yp(residual_yp_view.getData(), residual_yp_view.getData() + residual_yp_view.getSize()); + template + bool invalidParameterLeavesDefault(Params parameter, + ValueT invalid_value, + const char* label) const + { + auto expected_data = makeOrderData(4); + expected_data.parameters.erase(parameter); + auto invalid_data = expected_data; + invalid_data.parameters[parameter] = invalid_value; + + Fixture expected(expected_data); + Fixture invalid(invalid_data); + expected.u() = static_cast(0.5); + invalid.u() = static_cast(0.5); + + bool success = expected.prepare(); + success = success && invalid.model.allocate() == 0; + success = success && invalid.model.verify() > 0; + setAnswerKey(expected.model); + setAnswerKey(invalid.model); + success = success && expected.evaluate() == 0; + success = success && invalid.evaluate() == 0; + success = success && vectorMatches(invalid.model.getResidual(), copyVector(expected.model.getResidual()), label); + return success; + } - // Print dependencies for debugging - for (size_t i = 0; i < residual_y.size(); ++i) + bool checkTopology(const Data& data, size_t order, const char* label) const + { + Fixture fixture(data); + fixture.u() = static_cast(0.5); + if (!fixture.initialize()) { - 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"; + std::cout << "IEEEST " << label << " topology failed to initialize\n"; + return false; } - // 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) + setAnswerKey(fixture.model); + if (fixture.evaluate() != 0) { - auto dependency_y = (residual_y[i]).getDependencies(); - auto dependency_yp = (residual_yp[i]).getDependencies(); + return false; + } - for (const auto& pair_y : dependency_y) + bool success = true; + const auto* f = fixture.model.getResidual().getData(); + const auto* yp = fixture.model.yp().getData(); + for (size_t row = 0; row < 4; ++row) + { + const bool frozen = isEqual(static_cast(f[row]), + -static_cast(yp[row]), + kTol); + if (frozen != (row >= order)) { - 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)); - } + std::cout << "IEEEST " << label << " topology row " << row << " mismatch\n"; + success = false; } + } + return success; + } + + static constexpr std::array answerState() + { + return {{0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 0.05, 0.05}}; + } + + static constexpr std::array answerDerivative() + { + return {{0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.0, 0.0, 0.0, 0.0, 0.0}}; + } + + template + void setAnswerKey(PhasorDynamics::Stabilizer::Ieeest& model) const + { + const auto state = answerState(); + const auto derivative = answerDerivative(); + auto* y = model.y().getData(); + auto* yp = model.yp().getData(); + for (size_t row = 0; row < state.size(); ++row) + { + y[row] = static_cast(state[row]); + yp[row] = static_cast(derivative[row]); + } + model.y().setDataUpdated(); + model.yp().setDataUpdated(); + } + + bool checkResidual(const Data& data, + const std::array& expected, + const char* label) const + { + Fixture fixture(data); + fixture.u() = static_cast(0.5); + if (!fixture.initialize()) + { + return false; + } + setAnswerKey(fixture.model); + if (fixture.evaluate() != 0) + { + return false; + } - // Insert yp dependencies that did not exist in the y dependencies - for (const auto& pair_yp : dependency_yp) + bool success = true; + const auto* f = fixture.model.getResidual().getData(); + for (size_t row = 0; row < expected.size(); ++row) + { + const RealT tolerance = row == static_cast(Internal::VSS) ? kClampTol : kTol; + if (!rowMatches(static_cast(f[row]), expected[row], "residual", row, label, tolerance)) { - if (!dependency_y.contains(pair_yp.first)) - { - dependencies[i].insert(std::make_pair(pair_yp.first, pair_yp.second)); - } + success = false; } } + return success; + } - return dependencies; + template + bool allResidualsZero(const ModelT& model) const + { + bool success = true; + const auto* f = model.getResidual().getData(); + for (size_t row = 0; row < static_cast(model.getResidual().getSize()); ++row) + { + if (!rowMatches(static_cast(f[row]), 0.0, "residual", row, "at rest", kClampTol)) + { + success = false; + } + } + return success; } - std::vector EnzymeJacobian( - PhasorDynamics::Stabilizer::IeeestData ieeestdata) + bool rowMatches(RealT actual, + RealT expected, + const char* what, + size_t row, + const char* context, + RealT tolerance = kTol) const { - PhasorDynamics::SignalNode u_node; - PhasorDynamics::SignalNode vss_node; - ScalarT u_value{0.5}; - IdxT u_index = 12; - ScalarT vss_value{0.0}; - IdxT vss_index = INVALID_INDEX; - - u_node.set(&u_value, &u_index); - vss_node.set(&vss_value, &vss_index); - - PhasorDynamics::Stabilizer::Ieeest stab(ieeestdata); - stab.getSignals().template attachSignalNode(&u_node); - stab.getSignals().template assignSignalNode(&vss_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(); - stab.constructCsr(); - auto model_jacobian = stab.getCsrJacobian(); - std::cout << "Sparse Csr Matrix: Ieeest Jacobian\n"; - model_jacobian->print(); - - return GridKit::Testing::MapFromCsr(model_jacobian); + if (isEqual(actual, expected, tolerance)) + { + return true; + } + std::cout << "IEEEST " << what << " row " << kRowNames[row] << ' ' + << context << " mismatch: " + << std::setprecision(std::numeric_limits::max_digits10) + << actual << " != " << expected << '\n'; + return false; } -#endif - private: - static constexpr ScalarT tol_ = 10 * std::numeric_limits::epsilon(); + bool scalarMatches(RealT actual, + RealT expected, + const char* label, + RealT tolerance = kTol) const + { + if (isEqual(actual, expected, tolerance)) + { + return true; + } + std::cout << "IEEEST " << label << " mismatch: " + << std::setprecision(std::numeric_limits::max_digits10) + << actual << " != " << expected << '\n'; + return false; + } - /** - * @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 + template + std::vector copyVector(const VectorT& vector) const { - using Params = PhasorDynamics::Stabilizer::IeeestParameters; + const auto* values = vector.getData(); + std::vector copy(static_cast(vector.getSize())); + for (size_t row = 0; row < copy.size(); ++row) + { + copy[row] = static_cast(values[row]); + } + return copy; + } - PhasorDynamics::Stabilizer::IeeestData data; - data.device_class = "stabilizer"; - data.disambiguation_string = "ieeest_test"; - data.monitored_variables.insert(PhasorDynamics::Stabilizer::IeeestMonitorableVariables::vss); - - 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] = 0.0; - data.parameters[Params::Tdelay] = 0.0; + template + bool vectorMatches(const VectorT& vector, + const std::vector& expected, + const char* label) const + { + bool success = true; + const auto* values = vector.getData(); + for (size_t row = 0; row < expected.size(); ++row) + { + success *= rowMatches(static_cast(values[row]), + expected[row], + label, + row, + "atomic rejection"); + } + return success; + } - return data; + template + void poisonState(ModelT& model) const + { + auto* y = model.y().getData(); + auto* yp = model.yp().getData(); + for (size_t row = 0; row < static_cast(model.size()); ++row) + { + const RealT offset = static_cast(row); + y[row] = static_cast(0.125 + 0.01 * offset); + yp[row] = static_cast(-0.25 - 0.01 * offset); + } + model.y().setDataUpdated(); + model.yp().setDataUpdated(); } - /** - * @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) + void noteExpectedLogs(const char* message) const { - auto* y = stab.y().getData(); - auto* yp = stab.yp().getData(); - - y[0] = 0.1; // x1 - y[1] = 0.2; // x2 - y[2] = 0.3; // x3 - y[3] = 0.4; // x4 - y[4] = 0.5; // x5 - y[5] = 0.6; // x6 - y[6] = 0.7; // x7 - y[7] = 0.8; // v4 - y[8] = 0.9; // v5 - y[9] = 1.0; // v6 - y[10] = 0.05; // v7 (within limiter range) - y[11] = 0.05; // Vss (model output) - - yp[0] = 0.01; // x1_dot - yp[1] = 0.02; // x2_dot - yp[2] = 0.03; // x3_dot - yp[3] = 0.04; // x4_dot - yp[4] = 0.05; // x5_dot - yp[5] = 0.06; // x6_dot - yp[6] = 0.07; // x7_dot - - stab.y().setDataUpdated(); - stab.yp().setDataUpdated(); + const auto previous_verbosity = Log::verbosity(); + Log::setVerbosity(Log::Verbosity::EVERYTHING); + Log::misc() << message << '\n'; + Log::setVerbosity(previous_verbosity); } - /** - * @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) +#ifdef GRIDKIT_ENABLE_ENZYME + void printDependencyMap(const DependencyMap& row, const char* label) const { - auto* y = stab.y().getData(); - auto* yp = stab.yp().getData(); - - y[0].setValue(0.1); - y[1].setValue(0.2); - y[2].setValue(0.3); - y[3].setValue(0.4); - y[4].setValue(0.5); - y[5].setValue(0.6); - y[6].setValue(0.7); - y[7].setValue(0.8); - y[8].setValue(0.9); - y[9].setValue(1.0); - y[10].setValue(0.05); - y[11].setValue(0.05); - - yp[0].setValue(0.01); - yp[1].setValue(0.02); - yp[2].setValue(0.03); - yp[3].setValue(0.04); - yp[4].setValue(0.05); - yp[5].setValue(0.06); - yp[6].setValue(0.07); - - stab.y().setDataUpdated(); - stab.yp().setDataUpdated(); + std::cout << " " << label << ':'; + for (const auto& [column, value] : row) + { + std::cout << " (" << column << ", " + << std::setprecision(std::numeric_limits::max_digits10) + << value << ')'; + } + std::cout << '\n'; } - }; // class StabilizerIeeestTests + + std::vector dependencyTrackingJacobian(const Data& data, + TestStatus& success) const + { + using DepVar = DependencyTracking::Variable; + + Fixture fixture(data); + fixture.u() = static_cast(0.5); + success *= fixture.initialize(); + setAnswerKey(fixture.model); + + auto* y = fixture.model.y().getData(); + auto* yp = fixture.model.yp().getData(); + for (size_t row = 0; row < static_cast(fixture.model.size()); ++row) + { + y[row].setVariableNumber(row); + } + fixture.u().setVariableNumber(fixture.uIndex()); + fixture.model.y().setDataUpdated(); + success *= (fixture.evaluate() == 0); + + std::vector rows(static_cast(fixture.model.size())); + const auto* f = fixture.model.getResidual().getData(); + for (size_t row = 0; row < rows.size(); ++row) + { + rows[row] = f[row].getDependencies(); + } + + // Dependency tracking cannot safely alias y and yp to the same + // variable number in one pass. Evaluate df/dyp separately and add it + // to df/dy, matching Enzyme's df/dy + alpha*df/dyp at alpha = 1. + setAnswerKey(fixture.model); + fixture.u() = 0.5; + for (size_t row = 0; row < static_cast(fixture.model.size()); ++row) + { + yp[row].setVariableNumber(row); + } + fixture.model.yp().setDataUpdated(); + success *= (fixture.evaluate() == 0); + + f = fixture.model.getResidual().getData(); + for (size_t row = 0; row < rows.size(); ++row) + { + for (const auto& [column, value] : f[row].getDependencies()) + { + rows[row][column] += value; + } + } + return rows; + } + + std::vector enzymeJacobian(const Data& data, + TestStatus& success) const + { + Fixture fixture(data); + fixture.u() = static_cast(0.5); + success *= fixture.initialize(); + setAnswerKey(fixture.model); + fixture.model.updateTime(0.0, 1.0); + success *= (fixture.evaluate() == 0); + success *= (fixture.model.evaluateJacobian() == 0); + success *= (fixture.model.constructCsr() == 0); + return MapFromCsr(fixture.model.getCsrJacobian()); + } +#endif + }; } // namespace Testing } // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/runComponentConnectionTests.cpp b/tests/UnitTests/PhasorDynamics/runComponentConnectionTests.cpp index 9a56a471e..586e6b000 100644 --- a/tests/UnitTests/PhasorDynamics/runComponentConnectionTests.cpp +++ b/tests/UnitTests/PhasorDynamics/runComponentConnectionTests.cpp @@ -11,6 +11,7 @@ int main() result += test.genrouEsdc1a(); result += test.genrouHygov(); result += test.regcaRepca(); + result += test.ieeestIeeet1(); result += gastpti.genrouGastPti(); result += gastpti.gensalGastPti(); result += test.regcaReecb(); diff --git a/tests/UnitTests/PhasorDynamics/runStabilizerIeeestTests.cpp b/tests/UnitTests/PhasorDynamics/runStabilizerIeeestTests.cpp index c1999c24a..2b4022e59 100644 --- a/tests/UnitTests/PhasorDynamics/runStabilizerIeeestTests.cpp +++ b/tests/UnitTests/PhasorDynamics/runStabilizerIeeestTests.cpp @@ -6,9 +6,11 @@ int main() GridKit::Testing::StabilizerIeeestTests test; - result += test.constructor(); - result += test.zeroInitialResidual(); - result += test.residual(); + result += test.validation(); + result += test.initializationAndSignals(); + result += test.residualEquations(); + result += test.monitor(); + #ifdef GRIDKIT_ENABLE_ENZYME result += test.jacobian(); #endif From 5676d7834940744dc2183eb167466118bdb3fdd5 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Mon, 10 Aug 2026 18:16:24 -0500 Subject: [PATCH 2/6] remove param names from signature --- .../Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp index 14e4494da..4751bdd3e 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp @@ -105,11 +105,11 @@ namespace GridKit const Model::VariableMonitorBase* getMonitor() const override; __attribute__((always_inline)) inline int evaluateInternalResidual( - const ScalarT* y, - const ScalarT* yp, - const ScalarT* wb, - const ScalarT* ws, - ScalarT* f); + const ScalarT*, + const ScalarT*, + const ScalarT*, + const ScalarT*, + ScalarT*); private: void initializeParameters(const ModelDataT& data); From ebeff30502f895052548b4ea73960b46f614a9c1 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Mon, 10 Aug 2026 18:21:59 -0500 Subject: [PATCH 3/6] no ternary operators --- .../Stabilizer/IEEEST/IeeestImpl.hpp | 14 +++++++++----- .../PhasorDynamics/StabilizerIeeestTests.hpp | 6 +++++- 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp index 1c2b3fc6f..6f1b78a78 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp @@ -118,11 +118,15 @@ namespace GridKit && std::isfinite(a3_) && std::isfinite(a4_), "expanded denominator coefficients must be finite"); - const IdxT numerator_order = A6_ != ZERO - ? static_cast(2) - : A5_ != ZERO - ? static_cast(1) - : static_cast(0); + IdxT numerator_order = static_cast(0); + if (A6_ != ZERO) + { + numerator_order = static_cast(2); + } + else if (A5_ != ZERO) + { + numerator_order = static_cast(1); + } check(numerator_order <= order_, "numerator order must not exceed denominator order"); check(Lsmin_ < Lsmax_, "Lsmin must be less than Lsmax"); diff --git a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp index 4d0e73ead..28b077f11 100644 --- a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp +++ b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp @@ -725,7 +725,11 @@ namespace GridKit const auto* f = fixture.model.getResidual().getData(); for (size_t row = 0; row < expected.size(); ++row) { - const RealT tolerance = row == static_cast(Internal::VSS) ? kClampTol : kTol; + RealT tolerance = kTol; + if (row == static_cast(Internal::VSS)) + { + tolerance = kClampTol; + } if (!rowMatches(static_cast(f[row]), expected[row], "residual", row, label, tolerance)) { success = false; From 8920d539d4998e80bd7a0069032a4bb0ae67cb95 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Mon, 10 Aug 2026 18:27:11 -0500 Subject: [PATCH 4/6] remove operatorname maybe fix gh render --- GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md index d5280ead6..25886126f 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md @@ -206,7 +206,7 @@ $n$. The state layout does not change with $n$. 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_{\mathrm{ss}} - + \operatorname{clamp}(v_7,L_s^{\min},L_s^{\max}). + + \text{clamp}(v_7,L_s^{\min},L_s^{\max}). \end{aligned} ``` @@ -234,7 +234,7 @@ u \leftarrow \text{stabilizer input signal}. v_4,x_5,v_5,x_6,v_6,x_7 &\leftarrow u \\ v_7 &\leftarrow 0 \\ V_{\mathrm{ss}} &\leftarrow - \operatorname{clamp}(v_7,L_s^{\min},L_s^{\max}) \\ + \text{clamp}(v_7,L_s^{\min},L_s^{\max}) \\ \dot{x}_i &\leftarrow 0, \qquad i\in\{1,\ldots,7\}. \end{aligned} From e8652613a3188aa693bd180b050e51847001fc12 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Mon, 10 Aug 2026 18:40:51 -0500 Subject: [PATCH 5/6] test verbosity --- tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp index 28b077f11..577e9dc18 100644 --- a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp +++ b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp @@ -308,6 +308,10 @@ namespace GridKit kClampTol); } + noteExpectedLogs( + "Testing IEEEST rejected initialization paths. " + "The logged errors below are expected."); + // A nonfinite input rejects initialization without changing state. for (const RealT input : {std::numeric_limits::quiet_NaN(), std::numeric_limits::infinity(), From aa89aaaafe58200942e0a143441ac41c1d53b354 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Mon, 10 Aug 2026 18:52:31 -0500 Subject: [PATCH 6/6] remove single use method --- .../PhasorDynamics/StabilizerIeeestTests.hpp | 53 +++++++++---------- 1 file changed, 24 insertions(+), 29 deletions(-) diff --git a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp index 577e9dc18..802b5dce9 100644 --- a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp +++ b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp @@ -424,12 +424,21 @@ namespace GridKit const auto rows = std::min(dependency_jacobian.size(), enzyme_jacobian.size()); for (size_t row = 0; row < rows; ++row) { - if (!isEqual(dependency_jacobian[row], enzyme_jacobian[row], kTol)) + const auto& dependency_map = dependency_jacobian[row].getDependencies(); + if (!isEqual(dependency_map, enzyme_jacobian[row], kTol)) { std::cout << "IEEEST Jacobian row " << row << " for order " << order - << " differs between dependency tracking and Enzyme\n"; - printDependencyMap(dependency_jacobian[row], "dependency tracking"); - printDependencyMap(enzyme_jacobian[row], "Enzyme"); + << " differs between dependency tracking and Enzyme\n" + << " dependency tracking: " + << std::setprecision(std::numeric_limits::max_digits10); + dependency_jacobian[row].print(std::cout); + + std::cout << "\n Enzyme:"; + for (const auto& [column, value] : enzyme_jacobian[row]) + { + std::cout << " (" << column << ", " << value << ')'; + } + std::cout << '\n'; success = false; } } @@ -437,7 +446,8 @@ namespace GridKit for (size_t row = order; row < 4; ++row) { const DependencyMap expected{{row, static_cast(-1.0)}}; - if (!isEqual(dependency_jacobian[row], expected, kTol) + const auto& dependency_map = dependency_jacobian[row].getDependencies(); + if (!isEqual(dependency_map, expected, kTol) || !isEqual(enzyme_jacobian[row], expected, kTol)) { std::cout << "IEEEST inactive notch row " << row << " for order " @@ -844,20 +854,9 @@ namespace GridKit } #ifdef GRIDKIT_ENABLE_ENZYME - void printDependencyMap(const DependencyMap& row, const char* label) const - { - std::cout << " " << label << ':'; - for (const auto& [column, value] : row) - { - std::cout << " (" << column << ", " - << std::setprecision(std::numeric_limits::max_digits10) - << value << ')'; - } - std::cout << '\n'; - } - - std::vector dependencyTrackingJacobian(const Data& data, - TestStatus& success) const + std::vector dependencyTrackingJacobian( + const Data& data, + TestStatus& success) const { using DepVar = DependencyTracking::Variable; @@ -876,12 +875,8 @@ namespace GridKit fixture.model.y().setDataUpdated(); success *= (fixture.evaluate() == 0); - std::vector rows(static_cast(fixture.model.size())); - const auto* f = fixture.model.getResidual().getData(); - for (size_t row = 0; row < rows.size(); ++row) - { - rows[row] = f[row].getDependencies(); - } + const auto* f = fixture.model.getResidual().getData(); + std::vector rows(f, f + static_cast(fixture.model.size())); // Dependency tracking cannot safely alias y and yp to the same // variable number in one pass. Evaluate df/dyp separately and add it @@ -898,10 +893,10 @@ namespace GridKit f = fixture.model.getResidual().getData(); for (size_t row = 0; row < rows.size(); ++row) { - for (const auto& [column, value] : f[row].getDependencies()) - { - rows[row][column] += value; - } + // Add df/dyp dependencies without adding a second residual value. + DepVar yp_row = f[row]; + yp_row.setValue(0.0); + rows[row] += yp_row; } return rows; }