diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp index 4ef2dcb37..1162eee26 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp @@ -6,6 +6,10 @@ #pragma once +#include +#include +#include + #include #include #include @@ -119,6 +123,8 @@ namespace GridKit ScalarT*); private: + static constexpr RealT TIME_CONSTANT_MINIMUM = static_cast(1.0e-3); + RealT A1_{0}; RealT A2_{0}; RealT A3_{0}; @@ -131,6 +137,9 @@ namespace GridKit RealT T4_{1}; RealT T5_{0}; RealT T6_{1}; + RealT T2_inv_{1}; + RealT T4_inv_{1}; + RealT T6_inv_{1}; RealT Ks_{1}; RealT Lsmin_{-0.1}; RealT Lsmax_{0.1}; @@ -138,27 +147,23 @@ namespace GridKit RealT Vcu_{0}; RealT Tdelay_{0}; - RealT a0_{1}; RealT a1_{0}; RealT a2_{0}; RealT a3_{0}; RealT a4_{0}; - // Precomputed masks and safe inverse coefficients for branch-free degenerate paths. - RealT use_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 order_{0}; + RealT s0_{1}; + RealT s1_{0}; + RealT s2_{0}; + RealT s3_{0}; + RealT s4_{0}; + RealT a1_inv_{0}; + RealT a2_inv_{0}; + RealT a3_inv_{0}; + RealT a4_inv_{0}; + + IdxT parameter_error_count_{0}; ComponentSignals signals_; @@ -166,6 +171,7 @@ namespace GridKit void initializeParameters(const ModelDataT& data); void initializeMonitor(); + void setDerivedParameters(); std::vector ws_; std::vector ws_indices_; diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp index 97372cf06..785af0826 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp @@ -6,8 +6,9 @@ * @brief Definition of the IEEEST Power System Stabilizer. */ -#include +#include #include +#include #include #include @@ -26,7 +27,8 @@ namespace GridKit template Ieeest::Ieeest() { - size_ = 12; + size_ = static_cast(IeeestInternalVariables::MAXIMUM); + setDerivedParameters(); } template @@ -35,7 +37,8 @@ namespace GridKit { initializeParameters(data); initializeMonitor(); - size_ = 12; + size_ = static_cast(IeeestInternalVariables::MAXIMUM); + setDerivedParameters(); } template @@ -44,106 +47,98 @@ namespace GridKit } template - void Ieeest::initializeParameters(const ModelDataT& data) + void Ieeest::setDerivedParameters() { - if (data.parameters.contains(ModelDataT::Parameters::A1)) - { - A1_ = std::get(data.parameters.at(ModelDataT::Parameters::A1)); - } - if (data.parameters.contains(ModelDataT::Parameters::A2)) - { - A2_ = std::get(data.parameters.at(ModelDataT::Parameters::A2)); - } - if (data.parameters.contains(ModelDataT::Parameters::A3)) - { - A3_ = std::get(data.parameters.at(ModelDataT::Parameters::A3)); - } - if (data.parameters.contains(ModelDataT::Parameters::A4)) - { - A4_ = std::get(data.parameters.at(ModelDataT::Parameters::A4)); - } - if (data.parameters.contains(ModelDataT::Parameters::A5)) - { - A5_ = std::get(data.parameters.at(ModelDataT::Parameters::A5)); - } - if (data.parameters.contains(ModelDataT::Parameters::A6)) - { - A6_ = std::get(data.parameters.at(ModelDataT::Parameters::A6)); - } - if (data.parameters.contains(ModelDataT::Parameters::T1)) - { - T1_ = std::get(data.parameters.at(ModelDataT::Parameters::T1)); - } - if (data.parameters.contains(ModelDataT::Parameters::T2)) - { - T2_ = std::get(data.parameters.at(ModelDataT::Parameters::T2)); - } - if (data.parameters.contains(ModelDataT::Parameters::T3)) - { - T3_ = std::get(data.parameters.at(ModelDataT::Parameters::T3)); - } - if (data.parameters.contains(ModelDataT::Parameters::T4)) - { - T4_ = std::get(data.parameters.at(ModelDataT::Parameters::T4)); - } - if (data.parameters.contains(ModelDataT::Parameters::T5)) - { - T5_ = std::get(data.parameters.at(ModelDataT::Parameters::T5)); - } - if (data.parameters.contains(ModelDataT::Parameters::T6)) - { - T6_ = std::get(data.parameters.at(ModelDataT::Parameters::T6)); - } - if (data.parameters.contains(ModelDataT::Parameters::Ks)) - { - Ks_ = std::get(data.parameters.at(ModelDataT::Parameters::Ks)); - } - if (data.parameters.contains(ModelDataT::Parameters::Lsmin)) - { - Lsmin_ = std::get(data.parameters.at(ModelDataT::Parameters::Lsmin)); - } - if (data.parameters.contains(ModelDataT::Parameters::Lsmax)) + T2_ = std::max(T2_, TIME_CONSTANT_MINIMUM); + T4_ = std::max(T4_, TIME_CONSTANT_MINIMUM); + T6_ = std::max(T6_, TIME_CONSTANT_MINIMUM); + T2_inv_ = ONE / T2_; + T4_inv_ = ONE / T4_; + T6_inv_ = ONE / T6_; + + a1_ = A1_ + A3_; + a2_ = A2_ + A4_ + A1_ * A3_; + a3_ = A1_ * A4_ + A2_ * A3_; + a4_ = A2_ * A4_; + + order_ = 0; + if (a1_ != ZERO) { - Lsmax_ = std::get(data.parameters.at(ModelDataT::Parameters::Lsmax)); + order_ = 1; } - if (data.parameters.contains(ModelDataT::Parameters::Vcl)) + if (a2_ != ZERO) { - Vcl_ = std::get(data.parameters.at(ModelDataT::Parameters::Vcl)); + order_ = 2; } - if (data.parameters.contains(ModelDataT::Parameters::Vcu)) + if (a3_ != ZERO) { - Vcu_ = std::get(data.parameters.at(ModelDataT::Parameters::Vcu)); + order_ = 3; } - if (data.parameters.contains(ModelDataT::Parameters::Tdelay)) + if (a4_ != ZERO) { - Tdelay_ = std::get(data.parameters.at(ModelDataT::Parameters::Tdelay)); + order_ = 4; } - a0_ = 1; - a1_ = A1_ + A3_; - a2_ = A2_ + A4_ + A1_ * A3_; - a3_ = A1_ * A4_ + A2_ * A3_; - a4_ = A2_ * A4_; + s0_ = (order_ == 0) ? ONE : ZERO; + s1_ = (order_ == 1) ? ONE : ZERO; + s2_ = (order_ == 2) ? ONE : ZERO; + s3_ = (order_ == 3) ? ONE : ZERO; + s4_ = (order_ == 4) ? ONE : ZERO; - // 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_; + a1_inv_ = (order_ == 1) ? ONE / a1_ : ZERO; + a2_inv_ = (order_ == 2) ? ONE / a2_ : ZERO; + a3_inv_ = (order_ == 3) ? ONE / a3_ : ZERO; + a4_inv_ = (order_ == 4) ? ONE / a4_ : ZERO; + } - 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_)); + template + void Ieeest::initializeParameters(const ModelDataT& data) + { + using Params = typename ModelDataT::Parameters; - use_T2_block_ = static_cast(T2_ != 0.0); - bypass_T2_block_ = 1.0 - use_T2_block_; + parameter_error_count_ = 0; - use_T4_block_ = static_cast(T4_ != 0.0); - bypass_T4_block_ = 1.0 - use_T4_block_; + auto load_real = [&](auto key, RealT& target, const char* name) + { + if (!data.parameters.contains(key)) + { + return; + } - use_T6_block_ = static_cast(T6_ != 0.0); - bypass_T6_block_ = 1.0 - use_T6_block_; + const auto& value = data.parameters.at(key); + if (const auto* real_value = std::get_if(&value)) + { + target = *real_value; + } + else if (const auto* index_value = std::get_if(&value)) + { + target = static_cast(*index_value); + } + else + { + Log::error() << "Ieeest: parameter '" << name << "' must be numeric\n"; + ++parameter_error_count_; + } + }; + + 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"); + load_real(Params::Vcl, Vcl_, "Vcl"); + load_real(Params::Vcu, Vcu_, "Vcu"); + load_real(Params::Tdelay, Tdelay_, "Tdelay"); } template @@ -156,19 +151,20 @@ namespace GridKit template int Ieeest::allocate() { + size_ = static_cast(IeeestInternalVariables::MAXIMUM); auto size = static_cast(size_); - f_.resize(size); - y_.resize(size); - yp_.resize(size); - tag_.resize(size); - abs_tol_.resize(size); + + f_.assign(size, ScalarT{0}); + y_.assign(size, ScalarT{0}); + yp_.assign(size, ScalarT{0}); + tag_.assign(size, false); + abs_tol_.assign(size, ScalarT{0}); variable_indices_.resize(size); residual_indices_.resize(size); - ws_.resize(1); - ws_indices_.resize(1); - ws_[0] = 0.0; - ws_indices_[0] = INVALID_INDEX; + 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) { @@ -179,16 +175,19 @@ namespace GridKit if (signals_.template isAssigned()) { signals_.template getSignalNode()->set( - &y_[11], &(this->getVariableIndex(11))); + &y_[static_cast(IeeestInternalVariables::VSS)], + &(this->getVariableIndex(static_cast(IeeestInternalVariables::VSS)))); } + tagDifferentiable(); + return 0; } template int Ieeest::verify() const { - int ret = 0; + int ret = static_cast(parameter_error_count_); if (signals_.template isAttached()) { @@ -204,9 +203,9 @@ namespace GridKit ret += 1; } - if (a4_ == 0 && a3_ == 0 && a2_ == 0 && a1_ != 0) + if (order_ == 1 && A6_ != ZERO) { - Log::error() << "Ieeest: a2, a3, and a4 are all zero - no valid notch filter\n"; + Log::error() << "Ieeest: unsupported first-order notch filter with second-order numerator\n"; ret += 1; } @@ -216,30 +215,65 @@ namespace GridKit template int Ieeest::initialize() { - for (IdxT i = 0; i < size_; ++i) + if (verify() > 0) { - y_[static_cast(i)] = 0.0; - yp_[static_cast(i)] = 0.0; + Log::error() << "Ieeest: cannot initialize with invalid configuration\n"; + return 1; } + 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); + + std::fill(y_.begin(), y_.end(), ZERO); + std::fill(yp_.begin(), yp_.end(), ZERO); + + const ScalarT u = signals_.template readExternalVariable(); + ws_[U] = u; + ws_indices_[U] = signals_.template readExternalVariableIndex(); + + y_[X1] = u; + y_[X2] = ZERO; + y_[X3] = ZERO; + y_[X4] = ZERO; + 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_); + return 0; } 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; + auto index = [](IeeestInternalVariables variable) + { + return static_cast(variable); + }; + + std::fill(tag_.begin(), tag_.end(), false); + tag_[index(IeeestInternalVariables::X1)] = true; + tag_[index(IeeestInternalVariables::X2)] = true; + tag_[index(IeeestInternalVariables::X3)] = true; + tag_[index(IeeestInternalVariables::X4)] = true; + tag_[index(IeeestInternalVariables::X5)] = true; + tag_[index(IeeestInternalVariables::X6)] = true; + tag_[index(IeeestInternalVariables::X7)] = true; return 0; } @@ -271,43 +305,73 @@ 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]; + + const RealT s0 = s0_; + const RealT s1 = s1_; + const RealT s2 = s2_; + const RealT s3 = s3_; + const RealT s4 = s4_; + + const ScalarT x1_rhs = (-x1 + u) * a1_inv_; + const ScalarT x2_rhs = (-x1 - a1_ * x2 + u) * a2_inv_; + const ScalarT x3_rhs = (-x1 - a1_ * x2 - a2_ * x3 + u) * a3_inv_; + const ScalarT x4_rhs = (-x1 - a1_ * x2 - a2_ * x3 - a3_ * x4 + u) * a4_inv_; + const ScalarT x5_rhs = (v4 - x5) * T2_inv_; + const ScalarT x6_rhs = (v5 - x6) * T4_inv_; + const ScalarT x7_rhs = (v6 - x7) * T6_inv_; + + f[X1] = -x1_dot + s1 * x1_rhs + (s2 + s3 + s4) * x2; + f[X2] = -x2_dot + s2 * x2_rhs + (s3 + s4) * x3; + f[X3] = -x3_dot + s3 * x3_rhs + s4 * x4; + f[X4] = -x4_dot + s4 * x4_rhs; + f[X5] = -x5_dot + x5_rhs; + f[X6] = -x6_dot + x6_rhs; + f[X7] = -x7_dot + x7_rhs; + + f[V4] = -v4 + s0 * u + + s1 * (x1 + A5_ * x1_rhs) + + s2 * (x1 + A5_ * x2 + A6_ * x2_rhs) + + (s3 + s4) * (x1 + A5_ * x2 + A6_ * x3); + 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; } @@ -315,10 +379,15 @@ namespace GridKit template int Ieeest::evaluateResidual() { + const auto U = static_cast(IeeestExternalVariables::U); + + std::fill(ws_.begin(), ws_.end(), ZERO); + std::fill(ws_indices_.begin(), ws_indices_.end(), INVALID_INDEX); + if (signals_.template isAttached()) { - ws_[0] = signals_.template readExternalVariable(); - ws_indices_[0] = signals_.template readExternalVariableIndex(); + ws_[U] = signals_.template readExternalVariable(); + ws_indices_[U] = signals_.template readExternalVariableIndex(); } evaluateInternalResidual(y_.data(), yp_.data(), wb_.data(), ws_.data(), f_.data()); @@ -336,8 +405,13 @@ namespace GridKit void Ieeest::initializeMonitor() { using Variable = typename ModelDataT::MonitorableVariables; - monitor_->set(Variable::vss, [this] - { return y_[11]; }); + auto index = [](IeeestInternalVariables variable) + { + return static_cast(variable); + }; + + monitor_->set(Variable::vss, [this, index] + { return y_[index(IeeestInternalVariables::VSS)]; }); } } // namespace Stabilizer diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md index e5765ba1a..cdae661d3 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md @@ -3,6 +3,10 @@ Standard IEEE power system stabilizer: 4th-order notch filter, two lead–lag blocks, washout, and output limiter. +Notes: +- $V_{\mathrm{cl}}$, $V_{\mathrm{cu}}$, and $T_{\mathrm{delay}}$ are accepted for input-format + compatibility but are not modeled. + ## Block Diagram ![](../../../../../docs/Figures/stabilizer_ieeest_diagram.png) @@ -11,37 +15,57 @@ Figure 1: Stabilizer IEEEST model. Figure courtesy of [PowerWorld](https://www.p ## 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 | Accepted but not modeled +$V_{\mathrm{cu}}$ | [p.u.] | `Vcu` | Upper input cutout threshold | 0.0 | Accepted but not modeled +$T_{\mathrm{delay}}$ | [sec] | `Tdelay` | Input delay | 0.0 | Accepted but not modeled + +### Parameter Validation + +IEEEST denominator time constants are conditioned and unsupported notch forms are rejected by the following checks. Let $\epsilon_T=10^{-3}$. + +```math +\begin{aligned} + T &\leftarrow \max\!\left(T, \epsilon_T\right) + \quad T\in\{T_2,T_4,T_6\} \\ + A_6 &= 0 + \quad\text{when}\quad + n = 1 +\end{aligned} +``` + +### Model Derived Parameters ```math \begin{aligned} -a_0 &= 1 \\ -a_1 &= A_1 + A_3 \\ -a_2 &= A_2 + A_4 + A_1 A_3 \\ -a_3 &= A_1 A_4 + A_2 A_3 \\ -a_4 &= A_2 A_4 + 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 &= + \begin{cases} + 4 & a_4 \ne 0 \\ + 3 & a_4 = 0,\ a_3 \ne 0 \\ + 2 & a_4 = a_3 = 0,\ a_2 \ne 0 \\ + 1 & a_4 = a_3 = a_2 = 0,\ a_1 \ne 0 \\ + 0 & a_4 = a_3 = a_2 = a_1 = 0 + \end{cases} \end{aligned} ``` @@ -51,30 +75,34 @@ a_4 &= A_2 A_4 #### 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 +Symbol | Units | Description | Note +----------------------|--------|-----------------------|------ +$x_1, x_2, x_3, x_4$ | [-] | Notch filter states | States 1–4 in Fig. 1 +$x_5$ | [-] | Lead–lag 1 state | State 5 in Fig. 1 +$x_6$ | [-] | Lead–lag 2 state | State 6 in Fig. 1 +$x_7$ | [-] | 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 | ## Model Equations @@ -82,13 +110,50 @@ $u$ | [p.u.] | Stabilizer input signal ```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{1}{a_1}\left(-x_1 + u\right), & n = 1 \\ + x_2, & n = 2,3,4 + \end{cases} \\ + 0 &= + -\dot{x}_2 + + \begin{cases} + 0, + & n = 0,1 \\ + \dfrac{1}{a_2}\left(-x_1 - a_1x_2 + u\right), + & n = 2 \\ + x_3, + & n = 3,4 + \end{cases} \\ + 0 &= + -\dot{x}_3 + + \begin{cases} + 0, + & n = 0,1,2 \\ + \dfrac{1}{a_3}\left(-x_1 - a_1x_2 - a_2x_3 + u\right), + & n = 3 \\ + x_4, + & n = 4 + \end{cases} \\ + 0 &= + -\dot{x}_4 + + \begin{cases} + 0, + & n = 0,1,2,3 \\ + \dfrac{1}{a_4}\left(-x_1 - a_1x_2 - a_2x_3 - a_3x_4 + u\right), + & n = 4 + \end{cases} \\ + 0 &= + -\dot{x}_5 + + \dfrac{1}{T_2}\left(v_4 - x_5\right) \\ + 0 &= + -\dot{x}_6 + + \dfrac{1}{T_4}\left(v_5 - x_6\right) \\ + 0 &= + -\dot{x}_7 + + \dfrac{1}{T_6}\left(v_6 - x_7\right) \end{aligned} ``` @@ -96,19 +161,59 @@ $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}\left(-x_1 + u\right), + & n = 1 \\ + x_1 + A_5x_2 + + \dfrac{A_6}{a_2}\left(-x_1 - a_1x_2 + u\right), + & n = 2 \\ + x_1 + A_5x_2 + A_6x_3, + & n = 3,4 + \end{cases} \\ + 0 &= -v_5 + x_5 + \dfrac{T_1}{T_2}\left(v_4 - x_5\right) \\ + 0 &= -v_6 + x_6 + \dfrac{T_3}{T_4}\left(v_5 - x_6\right) \\ + 0 &= -v_7 + K_s\dfrac{T_5}{T_6}\left(v_6 - x_7\right) \\ + 0 &= + -V_{\mathrm{ss}} + + \text{clamp}(v_7, L_s^{\min}, L_s^{\max}) \end{aligned} ``` The output limiter uses GridKit's smooth -[Clamp](../../../../CommonMath.md#derived-functions). +[clamp](../../../../CommonMath.md#derived-functions). ## 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$. +### External Priors + +```math +\begin{aligned} + u_0 &\leftarrow \text{stabilizer input signal} +\end{aligned} +``` + +### Internal Initialization + +```math +\begin{aligned} + x_{1,0} &= v_{4,0} = x_{5,0} = v_{5,0} = x_{6,0} = v_{6,0} = x_{7,0} = u_0 \\ + x_{2,0} &= x_{3,0} = x_{4,0} = 0 \\ + v_{7,0} &= 0 \\ + V_{\mathrm{ss},0} &= \text{clamp}(v_{7,0}, L_s^{\min}, L_s^{\max}) \\ + \dot{x}_{i,0} &= 0 \quad i\in\{1,\ldots,7\} +\end{aligned} +``` + +### External Solved + +None. + +## Monitorable Outputs + +Output | Units | Description | Note +-------|--------|---------------------------|------ +`vss` | [p.u.] | Limited stabilizer signal | Exported through `output` when assigned diff --git a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp index 15b64876f..82fcf27db 100644 --- a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp +++ b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp @@ -1,8 +1,10 @@ #pragma once +#include #include #include #include +#include #include #include @@ -22,312 +24,505 @@ namespace GridKit { public: using RealT = typename PhasorDynamics::Component::RealT; + using DataT = PhasorDynamics::Stabilizer::IeeestData; StabilizerIeeestTests() = default; ~StabilizerIeeestTests() = default; - TestOutcome constructor() + TestOutcome init() { TestStatus success = true; - auto data = makeTestData(); - auto* stab = new PhasorDynamics::Stabilizer::Ieeest(data); + using Params = PhasorDynamics::Stabilizer::IeeestParameters; + + struct InitCase + { + RealT T6; + RealT Ks; + ScalarT u; + RealT Lsmin; + RealT Lsmax; + ScalarT expected_v7; + ScalarT expected_vss; + }; - success *= (stab != nullptr); - success *= (stab->getMonitor() != nullptr); + const auto loose_tol = static_cast(1.0e-4); + const std::vector cases = { + {0.0, 0.0, 0.25, -1.0, 1.0, 0.0, 0.0}, + {0.0, 4.0, 0.25, 0.2, 0.6, 0.0, 0.2}, + {5.0, 3.0, 0.25, -1.0, 1.0, 0.0, 0.0}, + }; - delete stab; + for (const auto& test : cases) + { + PhasorDynamics::SignalNode u_node; + PhasorDynamics::SignalNode vss_node; + ScalarT u_value{test.u}; + 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 = makeData(); + data.parameters[Params::T6] = test.T6; + data.parameters[Params::Ks] = test.Ks; + data.parameters[Params::Lsmin] = test.Lsmin; + data.parameters[Params::Lsmax] = test.Lsmax; + + PhasorDynamics::Stabilizer::Ieeest model(data); + model.getSignals().template attachSignalNode(&u_node); + model.getSignals().template assignSignalNode(&vss_node); + + model.allocate(); + success *= (model.verify() == 0); + model.initialize(); + + success *= vss_node.linked(); + success *= (vss_node.getVariableIndex() == 11); + success *= isEqual(model.y()[10], test.expected_v7, tol_); + success *= isEqual(model.y()[11], test.expected_vss, loose_tol); + success *= isEqual(vss_node.read(), test.expected_vss, loose_tol); + } return success.report(__func__); } - /** - * @brief All states initialize to zero (stabilizer at rest). - * With u = 0, all residuals should be zero. - */ - TestOutcome zeroInitialResidual() + TestOutcome residual() { 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; - - // Link signal nodes to backing storage - u_node.set(&u_value, &u_index); - vss_node.set(&vss_value, &vss_index); + using Params = PhasorDynamics::Stabilizer::IeeestParameters; - auto data = makeTestData(); - PhasorDynamics::Stabilizer::Ieeest stab(data); + struct ResidualCase + { + const char* name; + std::function edit; + const std::vector expected; + }; - // Wire: stabilizer reads u_node as input, writes vss_node as output - stab.getSignals().template attachSignalNode(&u_node); - stab.getSignals().template assignSignalNode(&vss_node); + const std::vector cases = { + {"baseline", + [](DataT&) {}, + {0.19, 0.28, 0.37, 1.0975, 0.25, 0.24, -0.01, -0.42, -0.25, -0.31, 1.15, 0.0}}, + {"a4_zero", + [](DataT& data) + { + data.parameters[Params::A4] = static_cast(0.0); + }, + {0.19, 0.28, 4.153333333333333, -0.04, 0.25, 0.24, -0.01, -0.42, -0.25, -0.31, 1.15, 0.0}}, + {"a3_a4_zero", + [](DataT& data) + { + data.parameters[Params::A3] = static_cast(0.0); + data.parameters[Params::A4] = static_cast(0.0); + }, + {0.19, 1.88, -0.03, -0.04, 0.25, 0.24, -0.01, 0.54, -0.25, -0.31, 1.15, 0.0}}, + {"notch_bypass", + [](DataT& data) + { + data.parameters[Params::A1] = static_cast(0.0); + data.parameters[Params::A2] = static_cast(0.0); + data.parameters[Params::A3] = static_cast(0.0); + data.parameters[Params::A4] = static_cast(0.0); + data.parameters[Params::A5] = static_cast(0.0); + data.parameters[Params::A6] = static_cast(0.0); + }, + {-0.01, -0.02, -0.03, -0.04, 0.25, 0.24, -0.01, -0.3, -0.25, -0.31, 1.15, 0.0}}, + {"first_order", + [](DataT& data) + { + data.parameters[Params::A2] = static_cast(0.0); + data.parameters[Params::A3] = static_cast(0.0); + data.parameters[Params::A4] = static_cast(0.0); + data.parameters[Params::A6] = static_cast(0.0); + }, + {3.99, -0.02, -0.03, -0.04, 0.25, 0.24, -0.01, 1.3, -0.25, -0.31, 1.15, 0.0}}, + }; - stab.allocate(); - success *= (stab.verify() == 0); - stab.initialize(); - stab.evaluateResidual(); + const auto loose_tol = static_cast(1.0e-4); - auto tol = 10 * std::numeric_limits::epsilon(); - const auto& f = stab.getResidual(); - for (size_t i = 0; i < f.size(); ++i) + for (const auto& test : cases) { - if (!isEqual(f[i], 0.0, tol)) + 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 = makeData(); + test.edit(data); + + PhasorDynamics::Stabilizer::Ieeest model(data); + model.getSignals().template attachSignalNode(&u_node); + model.getSignals().template assignSignalNode(&vss_node); + + model.allocate(); + model.initialize(); + + model.y()[0] = 0.1; + model.y()[1] = 0.2; + model.y()[2] = 0.3; + model.y()[3] = 0.4; + model.y()[4] = 0.5; + model.y()[5] = 0.6; + model.y()[6] = 0.7; + model.y()[7] = 0.8; + model.y()[8] = 0.9; + model.y()[9] = 1.0; + model.y()[10] = 0.05; + model.y()[11] = 0.05; + + model.yp()[0] = 0.01; + model.yp()[1] = 0.02; + model.yp()[2] = 0.03; + model.yp()[3] = 0.04; + model.yp()[4] = 0.05; + model.yp()[5] = 0.06; + model.yp()[6] = 0.07; + + model.evaluateResidual(); + + for (size_t i = 0; i < test.expected.size(); ++i) { - std::cout << "Non-zero residual at index " << i << ": " << f[i] << "\n"; - success = false; + auto test_tol = (i == 11) ? loose_tol : tol_; + if (!isEqual(model.getResidual()[i], test.expected[i], test_tol)) + { + std::cout << "Incorrect residual for " << test.name + << " row " << i << ": " + << std::setprecision(15) << model.getResidual()[i] + << " != " << test.expected[i] << "\n"; + success = false; + } } } - // 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); + { + 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 = makeData(); + data.parameters[Params::A2] = static_cast(0.0); + data.parameters[Params::A4] = static_cast(0.0); + data.parameters[Params::T2] = static_cast(0.0); + data.parameters[Params::T4] = static_cast(0.0); + data.parameters[Params::T6] = static_cast(0.0); + + PhasorDynamics::Stabilizer::Ieeest model(data); + model.getSignals().template attachSignalNode(&u_node); + model.getSignals().template assignSignalNode(&vss_node); + + model.allocate(); + success *= (model.verify() == 0); + success *= (model.initialize() == 0); + success *= (model.evaluateResidual() == 0); + + const auto loose_tol = static_cast(1.0e-4); + for (size_t i = 0; i < model.getResidual().size(); ++i) + { + if (!isEqual(model.getResidual()[i], static_cast(0.0), loose_tol)) + { + std::cout << "Conditioned-constant residual row " << i << " is " + << std::setprecision(15) << model.getResidual()[i] << "\n"; + success = false; + } + } + } 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 tags() { 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; + using Params = PhasorDynamics::Stabilizer::IeeestParameters; + PhasorDynamics::SignalNode u_node; + ScalarT u_value{0.0}; + IdxT u_index{12}; 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(); + auto data = makeData(); + data.parameters[Params::A2] = static_cast(0.0); + data.parameters[Params::A4] = static_cast(0.0); + data.parameters[Params::T2] = static_cast(0.0); + data.parameters[Params::T4] = static_cast(0.0); + data.parameters[Params::T6] = static_cast(0.0); + + PhasorDynamics::Stabilizer::Ieeest model(data); + model.getSignals().template attachSignalNode(&u_node); + + model.allocate(); - for (size_t i = 0; i < res_answer.size(); ++i) + for (size_t i = 0; i < model.tag().size(); ++i) { - auto test_tol = (i == 11) ? loose_tol : static_cast(10 * std::numeric_limits::epsilon()); - if (!isEqual(residual[i], res_answer[i], test_tol)) + const bool expected = (i <= static_cast(PhasorDynamics::Stabilizer::IeeestInternalVariables::X7)); + if (model.tag()[i] != expected) { - std::cout << "Incorrect result for residual " << i << ": " - << std::setprecision(15) << residual[i] - << " != " << res_answer[i] << "\n"; + std::cout << "Incorrect differential tag at row " << i << "\n"; success = false; } } - // 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() + TestOutcome verify() { TestStatus success = true; + using Params = PhasorDynamics::Stabilizer::IeeestParameters; - auto data = makeTestData(); + { + PhasorDynamics::Stabilizer::Ieeest model(makeData()); + model.allocate(); + success *= (model.verify() != 0); + } - std::vector - dependency_tracking_jacobian = DependencyTrackingJacobian(data); + { + PhasorDynamics::SignalNode u_node; + PhasorDynamics::Stabilizer::Ieeest model(makeData()); + model.getSignals().template attachSignalNode(&u_node); + model.allocate(); + success *= (model.verify() != 0); + } - std::vector - enzyme_jacobian = EnzymeJacobian(data); + { + PhasorDynamics::SignalNode u_node; + ScalarT u_value{0.0}; + IdxT u_index{12}; + u_node.set(&u_value, &u_index); + + auto data = makeData(); + data.parameters[Params::A1] = static_cast(1.0); + data.parameters[Params::A2] = static_cast(0.0); + data.parameters[Params::A3] = static_cast(0.0); + data.parameters[Params::A4] = static_cast(0.0); + + PhasorDynamics::Stabilizer::Ieeest model(data); + model.getSignals().template attachSignalNode(&u_node); + model.allocate(); + success *= (model.verify() != 0); + } - // 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)); + PhasorDynamics::SignalNode u_node; + ScalarT u_value{0.0}; + IdxT u_index{12}; + u_node.set(&u_value, &u_index); + + auto data = makeData(); + data.parameters[Params::A1] = static_cast(1.0); + data.parameters[Params::A2] = static_cast(0.0); + data.parameters[Params::A3] = static_cast(0.0); + data.parameters[Params::A4] = static_cast(0.0); + data.parameters[Params::A6] = static_cast(0.0); + + PhasorDynamics::Stabilizer::Ieeest model(data); + model.getSignals().template attachSignalNode(&u_node); + model.allocate(); + success *= (model.verify() == 0); } return success.report(__func__); } - private: - std::vector DependencyTrackingJacobian( - PhasorDynamics::Stabilizer::IeeestData ieeestdata) +#ifdef GRIDKIT_ENABLE_ENZYME + TestOutcome jacobian() { - using DepVar = DependencyTracking::Variable; - - // 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; - - 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); + TestStatus success = true; + using DepVar = DependencyTracking::Variable; - stab.allocate(); - stab.initialize(); + std::vector dependency_tracking_jacobian; - // --- d/dy: tag internal variables as independent --- - for (size_t i = 0; i < stab.size(); ++i) { - stab.y()[i].setVariableNumber(i); - } - // Tag external signal u as an additional independent variable - u_value.setVariableNumber(stab.size()); - u_value.setValue(0.5); - - setStatePointDep(stab); + 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}; - stab.evaluateResidual(); - std::vector residual_y = stab.getResidual(); - - // --- d/dy': tag derivatives as independent --- - stab.initialize(); - for (size_t i = 0; i < stab.size(); ++i) - { - stab.yp()[i].setVariableNumber(i); - } + u_node.set(&u_value, &u_index); + vss_node.set(&vss_value, &vss_index); - u_value = 0.5; - setStatePointDep(stab); + PhasorDynamics::Stabilizer::Ieeest model(makeData()); + model.getSignals().template attachSignalNode(&u_node); + model.getSignals().template assignSignalNode(&vss_node); - stab.evaluateResidual(); - std::vector residual_yp = stab.getResidual(); + model.allocate(); + model.initialize(); - // Print dependencies for debugging - for (size_t i = 0; i < residual_y.size(); ++i) - { - std::cout << i << "th residual, y: "; - (residual_y[i]).print(std::cout); - std::cout << "\n"; - std::cout << i << "th residual, yp: "; - (residual_yp[i]).print(std::cout); - std::cout << "\n"; - } - - // Merge d/dy and d/dy' into a single dependency map - std::vector dependencies(residual_y.size()); - for (IdxT i = 0; i < residual_y.size(); ++i) - { - auto dependency_y = (residual_y[i]).getDependencies(); - auto dependency_yp = (residual_yp[i]).getDependencies(); - - for (const auto& pair_y : dependency_y) + for (size_t i = 0; i < model.size(); ++i) { - 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 + model.y()[i].setVariableNumber(i); + } + u_value.setVariableNumber(model.size()); + u_value.setValue(0.5); + + model.y()[0].setValue(0.1); + model.y()[1].setValue(0.2); + model.y()[2].setValue(0.3); + model.y()[3].setValue(0.4); + model.y()[4].setValue(0.5); + model.y()[5].setValue(0.6); + model.y()[6].setValue(0.7); + model.y()[7].setValue(0.8); + model.y()[8].setValue(0.9); + model.y()[9].setValue(1.0); + model.y()[10].setValue(0.05); + model.y()[11].setValue(0.05); + + model.yp()[0].setValue(0.01); + model.yp()[1].setValue(0.02); + model.yp()[2].setValue(0.03); + model.yp()[3].setValue(0.04); + model.yp()[4].setValue(0.05); + model.yp()[5].setValue(0.06); + model.yp()[6].setValue(0.07); + + model.evaluateResidual(); + std::vector residual_y = model.getResidual(); + + model.initialize(); + for (size_t i = 0; i < model.size(); ++i) + { + model.y()[i] = model.y()[i].getValue(); + model.yp()[i].setVariableNumber(i); + } + u_value = 0.5; + + model.y()[0].setValue(0.1); + model.y()[1].setValue(0.2); + model.y()[2].setValue(0.3); + model.y()[3].setValue(0.4); + model.y()[4].setValue(0.5); + model.y()[5].setValue(0.6); + model.y()[6].setValue(0.7); + model.y()[7].setValue(0.8); + model.y()[8].setValue(0.9); + model.y()[9].setValue(1.0); + model.y()[10].setValue(0.05); + model.y()[11].setValue(0.05); + + model.yp()[0].setValue(0.01); + model.yp()[1].setValue(0.02); + model.yp()[2].setValue(0.03); + model.yp()[3].setValue(0.04); + model.yp()[4].setValue(0.05); + model.yp()[5].setValue(0.06); + model.yp()[6].setValue(0.07); + + model.evaluateResidual(); + std::vector residual_yp = model.getResidual(); + + dependency_tracking_jacobian.resize(residual_y.size()); + for (IdxT i = 0; i < residual_y.size(); ++i) + { + auto dependency_y = residual_y[i].getDependencies(); + auto dependency_yp = residual_yp[i].getDependencies(); + + for (const auto& pair_y : dependency_y) { - dependencies[i].insert(std::make_pair(pair_y.first, pair_y.second)); + auto it_yp = dependency_yp.find(pair_y.first); + if (it_yp != dependency_yp.end()) + { + dependency_tracking_jacobian[i].insert(std::make_pair(pair_y.first, pair_y.second + it_yp->second)); + } + else + { + dependency_tracking_jacobian[i].insert(std::make_pair(pair_y.first, pair_y.second)); + } } - } - // Insert yp dependencies that did not exist in the y dependencies - for (const auto& pair_yp : dependency_yp) - { - if (!dependency_y.contains(pair_yp.first)) + for (const auto& pair_yp : dependency_yp) { - dependencies[i].insert(std::make_pair(pair_yp.first, pair_yp.second)); + if (!dependency_y.contains(pair_yp.first)) + { + dependency_tracking_jacobian[i].insert(std::make_pair(pair_yp.first, pair_yp.second)); + } } } } - return dependencies; - } - - std::vector EnzymeJacobian( - PhasorDynamics::Stabilizer::IeeestData ieeestdata) - { - PhasorDynamics::SignalNode u_node; - PhasorDynamics::SignalNode 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); + std::vector enzyme_jacobian; - stab.updateTime(0.0, 1.0); // alpha = 1.0 to verify d/dy' term + { + 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 model(makeData()); + model.getSignals().template attachSignalNode(&u_node); + model.getSignals().template assignSignalNode(&vss_node); + + model.allocate(); + model.initialize(); + + model.y()[0] = 0.1; + model.y()[1] = 0.2; + model.y()[2] = 0.3; + model.y()[3] = 0.4; + model.y()[4] = 0.5; + model.y()[5] = 0.6; + model.y()[6] = 0.7; + model.y()[7] = 0.8; + model.y()[8] = 0.9; + model.y()[9] = 1.0; + model.y()[10] = 0.05; + model.y()[11] = 0.05; + + model.yp()[0] = 0.01; + model.yp()[1] = 0.02; + model.yp()[2] = 0.03; + model.yp()[3] = 0.04; + model.yp()[4] = 0.05; + model.yp()[5] = 0.06; + model.yp()[6] = 0.07; + + model.updateTime(0.0, 1.0); + model.evaluateResidual(); + model.evaluateJacobian(); + model.constructCsr(); + auto model_jacobian = model.getCsrJacobian(); + enzyme_jacobian = GridKit::Testing::MapFromCsr(model_jacobian); + } - stab.evaluateResidual(); - stab.evaluateJacobian(); - stab.constructCsr(); - auto model_jacobian = stab.getCsrJacobian(); - std::cout << "Sparse Csr Matrix: Ieeest Jacobian\n"; - model_jacobian->print(); + for (size_t i = 0; i < dependency_tracking_jacobian.size(); ++i) + { + success *= GridKit::Testing::isEqual(dependency_tracking_jacobian[i], enzyme_jacobian[i], tol_); + } - return GridKit::Testing::MapFromCsr(model_jacobian); + return success.report(__func__); } #endif private: static constexpr ScalarT tol_ = 10 * std::numeric_limits::epsilon(); - /** - * @brief Standard IEEEST parameter set for all tests. - * Derived: a0=1, a1=0.4, a2=0.63, a3=0.1, a4=0.08 - */ - auto makeTestData() -> PhasorDynamics::Stabilizer::IeeestData + auto makeData() -> DataT { using Params = PhasorDynamics::Stabilizer::IeeestParameters; - PhasorDynamics::Stabilizer::IeeestData data; + DataT data; data.device_class = "stabilizer"; data.disambiguation_string = "ieeest_test"; data.monitored_variables.insert(PhasorDynamics::Stabilizer::IeeestMonitorableVariables::vss); @@ -353,62 +548,6 @@ namespace GridKit return data; } - - /** - * @brief Set a non-trivial operating point for residual/Jacobian tests. - * Avoids zeros and ones to catch coefficient errors. - */ - void setStatePoint(PhasorDynamics::Stabilizer::Ieeest& stab) - { - stab.y()[0] = 0.1; // x1 - stab.y()[1] = 0.2; // x2 - stab.y()[2] = 0.3; // x3 - stab.y()[3] = 0.4; // x4 - stab.y()[4] = 0.5; // x5 - stab.y()[5] = 0.6; // x6 - stab.y()[6] = 0.7; // x7 - stab.y()[7] = 0.8; // v4 - stab.y()[8] = 0.9; // v5 - stab.y()[9] = 1.0; // v6 - stab.y()[10] = 0.05; // v7 (within limiter range) - stab.y()[11] = 0.05; // Vss (model output) - - stab.yp()[0] = 0.01; // x1_dot - stab.yp()[1] = 0.02; // x2_dot - stab.yp()[2] = 0.03; // x3_dot - stab.yp()[3] = 0.04; // x4_dot - stab.yp()[4] = 0.05; // x5_dot - stab.yp()[5] = 0.06; // x6_dot - stab.yp()[6] = 0.07; // x7_dot - } - - /** - * @brief Set the same operating point for DependencyTracking variables. - * Uses setValue() to set the numeric value while preserving dependency info. - */ - void setStatePointDep(PhasorDynamics::Stabilizer::Ieeest& stab) - { - stab.y()[0].setValue(0.1); - stab.y()[1].setValue(0.2); - stab.y()[2].setValue(0.3); - stab.y()[3].setValue(0.4); - stab.y()[4].setValue(0.5); - stab.y()[5].setValue(0.6); - stab.y()[6].setValue(0.7); - stab.y()[7].setValue(0.8); - stab.y()[8].setValue(0.9); - stab.y()[9].setValue(1.0); - stab.y()[10].setValue(0.05); - stab.y()[11].setValue(0.05); - - stab.yp()[0].setValue(0.01); - stab.yp()[1].setValue(0.02); - stab.yp()[2].setValue(0.03); - stab.yp()[3].setValue(0.04); - stab.yp()[4].setValue(0.05); - stab.yp()[5].setValue(0.06); - stab.yp()[6].setValue(0.07); - } }; // class StabilizerIeeestTests } // namespace Testing diff --git a/tests/UnitTests/PhasorDynamics/runStabilizerIeeestTests.cpp b/tests/UnitTests/PhasorDynamics/runStabilizerIeeestTests.cpp index c1999c24a..2a8f60f89 100644 --- a/tests/UnitTests/PhasorDynamics/runStabilizerIeeestTests.cpp +++ b/tests/UnitTests/PhasorDynamics/runStabilizerIeeestTests.cpp @@ -6,9 +6,10 @@ int main() GridKit::Testing::StabilizerIeeestTests test; - result += test.constructor(); - result += test.zeroInitialResidual(); + result += test.init(); result += test.residual(); + result += test.tags(); + result += test.verify(); #ifdef GRIDKIT_ENABLE_ENZYME result += test.jacobian(); #endif