From c0878bdf2dd2d4539e235497f88f031a44730907 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Mon, 8 Jun 2026 04:44:18 -0500 Subject: [PATCH 1/4] Fix bad parameter initialization and polish docs --- .../Stabilizer/IEEEST/Ieeest.hpp | 18 +- .../Stabilizer/IEEEST/IeeestImpl.hpp | 105 +-- .../Stabilizer/IEEEST/README.md | 140 +++- .../PhasorDynamics/StabilizerIeeestTests.hpp | 620 +++++++++--------- .../runStabilizerIeeestTests.cpp | 4 +- 5 files changed, 497 insertions(+), 390 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp index 4ef2dcb37..399398ae8 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp @@ -138,34 +138,18 @@ 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}; - ComponentSignals signals_; std::unique_ptr monitor_; 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..ec2649aa5 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp @@ -43,6 +43,15 @@ namespace GridKit { } + template + void Ieeest::setDerivedParameters() + { + a1_ = A1_ + A3_; + a2_ = A2_ + A4_ + A1_ * A3_; + a3_ = A1_ * A4_ + A2_ * A3_; + a4_ = A2_ * A4_; + } + template void Ieeest::initializeParameters(const ModelDataT& data) { @@ -119,31 +128,7 @@ namespace GridKit Tdelay_ = std::get(data.parameters.at(ModelDataT::Parameters::Tdelay)); } - a0_ = 1; - a1_ = A1_ + A3_; - a2_ = A2_ + A4_ + A1_ * A3_; - a3_ = A1_ * A4_ + A2_ * A3_; - a4_ = A2_ * A4_; - - // Precompute masks and safe inverse coefficients so the residual stays branch-free. - use_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_; + setDerivedParameters(); } template @@ -182,6 +167,8 @@ namespace GridKit &y_[11], &(this->getVariableIndex(11))); } + tagDifferentiable(); + return 0; } @@ -204,7 +191,7 @@ namespace GridKit ret += 1; } - if (a4_ == 0 && a3_ == 0 && a2_ == 0 && a1_ != 0) + if (a2_ == ZERO && a3_ == ZERO && a4_ == ZERO && a1_ != ZERO) { Log::error() << "Ieeest: a2, a3, and a4 are all zero - no valid notch filter\n"; ret += 1; @@ -222,19 +209,53 @@ namespace GridKit yp_[static_cast(i)] = 0.0; } + ScalarT u{0.0}; + if (signals_.template isAttached()) + { + u = signals_.template readExternalVariable(); + ws_[0] = u; + ws_indices_[0] = signals_.template readExternalVariableIndex(); + } + + const ScalarT zero{0.0}; + const ScalarT x1 = u; + const ScalarT x2 = zero; + const ScalarT x3 = zero; + const ScalarT x4 = zero; + const ScalarT v4 = x1 + A5_ * x2 + A6_ * x3; + const ScalarT x5 = v4; + const ScalarT v5 = v4; + const ScalarT x6 = v5; + const ScalarT v6 = v5; + const ScalarT x7 = v6; + const ScalarT v7 = (T6_ == ZERO) ? Ks_ * v6 : zero; + + y_[0] = x1; + y_[1] = x2; + y_[2] = x3; + y_[3] = x4; + y_[4] = x5; + y_[5] = x6; + y_[6] = x7; + y_[7] = v4; + y_[8] = v5; + y_[9] = v6; + y_[10] = v7; + y_[11] = Math::clamp(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_[0] = (a2_ != ZERO || a3_ != ZERO || a4_ != ZERO); + tag_[1] = tag_[0]; + tag_[2] = (a3_ != ZERO || a4_ != ZERO); + tag_[3] = (a4_ != ZERO); + tag_[4] = (T2_ != ZERO); + tag_[5] = (T4_ != ZERO); + tag_[6] = (T6_ != ZERO); tag_[7] = false; tag_[8] = false; tag_[9] = false; @@ -294,19 +315,17 @@ namespace GridKit 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[0] = -tag_[0] * x1_dot + x2; + f[1] = -tag_[1] * x2_dot + x3; + f[2] = -tag_[2] * x3_dot + x4; + f[3] = -a4_ * x4_dot - x1 - a1_ * x2 - a2_ * x3 - a3_ * x4 + u; 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[7] = -v4 + x1 + A5_ * x2 + A6_ * x3; + f[8] = tag_[4] * (-T2_ * (v5 - x5) + T1_ * (v4 - x5)) + (1 - tag_[4]) * (v4 - v5); + f[9] = tag_[5] * (-T4_ * (v6 - x6) + T3_ * (v5 - x6)) + (1 - tag_[5]) * (v5 - v6); + f[10] = tag_[6] * (-T6_ * v7 + Ks_ * T5_ * (v6 - x7)) + (1 - tag_[6]) * (Ks_ * v6 - v7); f[11] = -vss + Math::clamp(v7, Lsmin_, Lsmax_); return 0; diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md index e5765ba1a..d1d5a8870 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md @@ -3,6 +3,12 @@ Standard IEEE power system stabilizer: 4th-order notch filter, two lead–lag blocks, washout, and output limiter. +Notes: +- $V_{cl}$, $V_{cu}$, and $T_{delay}$ are accepted for input-format + compatibility but are not modeled. +- A zero denominator time constant bypasses its corresponding lead–lag or + washout block. + ## Block Diagram ![](../../../../../docs/Figures/stabilizer_ieeest_diagram.png) @@ -11,33 +17,43 @@ 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$ | [s] | `A1` | Notch denominator coefficient | 1.013 | +$A_2$ | [s²] | `A2` | Notch denominator coefficient | 0.013 | +$A_3$ | [s] | `A3` | Notch denominator coefficient | 0.0 | +$A_4$ | [s²] | `A4` | Notch denominator coefficient | 0.0 | +$A_5$ | [s] | `A5` | Notch numerator coefficient | 1.013 | +$A_6$ | [s²] | `A6` | Notch numerator coefficient | 0.113 | +$T_1$ | [s] | `T1` | Lead–lag 1 numerator time constant | 0.0 | +$T_2$ | [s] | `T2` | Lead–lag 1 denominator time constant | 0.02 | +$T_3$ | [s] | `T3` | Lead–lag 2 numerator time constant | 0.0 | +$T_4$ | [s] | `T4` | Lead–lag 2 denominator time constant | 0.0 | +$T_5$ | [s] | `T5` | Washout numerator time constant | 1.65 | +$T_6$ | [s] | `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_{cl}$ | [p.u.] | `Vcl` | Lower input cutout threshold | 0.0 | Accepted but not modeled +$V_{cu}$ | [p.u.] | `Vcu` | Upper input cutout threshold | 0.0 | Accepted but not modeled +$T_\text{delay}$ | [s] | `Tdelay` | Input delay | 0.0 | Accepted but not modeled + +### Parameter Validation + +The fixed realization rejects a first-order-only notch denominator: + +```math +\begin{aligned} + a_2 \ne 0 \lor a_3 \ne 0 \lor a_4 \ne 0 \lor a_1 = 0 +\end{aligned} +``` + +### Model Derived Parameters + +The notch-filter denominator expands to: ```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 \\ @@ -45,6 +61,24 @@ a_4 &= A_2 A_4 \end{aligned} ``` +The binary DAE selectors choose the active notch-filter order: + +```math +\begin{aligned} +\delta_1 &= \delta_2 = +\begin{cases} +1 & a_2 \ne 0 \lor a_3 \ne 0 \lor a_4 \ne 0 \\ +0 & \text{otherwise} +\end{cases} +\\ +\delta_3 &= +\begin{cases} +1 & a_3 \ne 0 \lor a_4 \ne 0 \\ +0 & \text{otherwise} +\end{cases} +\end{aligned} +``` + ## Model Variables ### Internal Variables @@ -58,6 +92,9 @@ $x_5$ | [-] | Lead–lag 1 state $x_6$ | [-] | Lead–lag 2 state $x_7$ | [-] | Washout state +For reduced-order notch filters, unused notch states remain in the fixed +component state vector and are pinned by algebraic residuals. + #### Algebraic Symbol | Units | Description @@ -82,10 +119,10 @@ $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 &= -\delta_1\dot{x}_1 + x_2 \\ +0 &= -\delta_2\dot{x}_2 + x_3 \\ +0 &= -\delta_3\dot{x}_3 + x_4 \\ +0 &= -a_4\dot{x}_4 - x_1 - a_1x_2 - a_2x_3 - a_3x_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 @@ -96,10 +133,22 @@ $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_4 + x_1 + A_5x_2 + A_6x_3 \\ +0 &= +\begin{cases} +-T_2(v_5 - x_5) + T_1(v_4 - x_5) & T_2 \ne 0 \\ +v_4 - v_5 & T_2 = 0 +\end{cases} \\ +0 &= +\begin{cases} +-T_4(v_6 - x_6) + T_3(v_5 - x_6) & T_4 \ne 0 \\ +v_5 - v_6 & T_4 = 0 +\end{cases} \\ +0 &= +\begin{cases} +-T_6 v_7 + K_s T_5(v_6 - x_7) & T_6 \ne 0 \\ +K_s v_6 - v_7 & T_6 = 0 +\end{cases} \\ 0 &= -V_{ss} + \text{clamp}(v_7, L_s^{\min}, L_s^{\max}) \end{aligned} ``` @@ -109,6 +158,27 @@ The output limiter uses GridKit's smooth ## 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$. +States and derivatives initialize to the steady state implied by the attached +input $u$: + +```math +\begin{aligned} +x_1 &= v_4 = x_5 = v_5 = x_6 = v_6 = x_7 = u \\ +x_2 &= x_3 = x_4 = 0 \\ +v_7 &= +\begin{cases} +K_su & T_6 = 0 \\ +0 & \text{otherwise} +\end{cases} +& +V_{ss} &= \text{clamp}(v_7, L_s^{\min}, L_s^{\max}) +\end{aligned} +``` + +All internal derivatives initialize to zero. + +## Model 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..e960ae308 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,400 @@ 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); - - success *= (stab != nullptr); - success *= (stab->getMonitor() != nullptr); - - delete stab; + using Params = PhasorDynamics::Stabilizer::IeeestParameters; - return success.report(__func__); - } + struct InitCase + { + RealT T6; + RealT Ks; + ScalarT u; + RealT Lsmin; + RealT Lsmax; + ScalarT raw_v7; + ScalarT expected_vss; + }; - /** - * @brief All states initialize to zero (stabilizer at rest). - * With u = 0, all residuals should be zero. - */ - TestOutcome zeroInitialResidual() - { - TestStatus success = true; + 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, 1.0, 0.25, -1.0, 1.0, 0.25, 0.25}, + {0.0, 2.0, 0.25, -1.0, 1.0, 0.50, 0.50}, + {0.0, 4.0, 0.25, -1.0, 0.6, 1.00, 0.60}, + {5.0, 3.0, 0.25, -1.0, 1.0, 0.00, 0.00}, + }; - // 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); - - auto data = makeTestData(); - PhasorDynamics::Stabilizer::Ieeest stab(data); - - // Wire: stabilizer reads u_node as input, writes vss_node as output - stab.getSignals().template attachSignalNode(&u_node); - stab.getSignals().template assignSignalNode(&vss_node); - - stab.allocate(); - success *= (stab.verify() == 0); - stab.initialize(); - stab.evaluateResidual(); - - auto tol = 10 * std::numeric_limits::epsilon(); - const auto& f = stab.getResidual(); - for (size_t i = 0; i < f.size(); ++i) + for (const auto& test : cases) { - if (!isEqual(f[i], 0.0, tol)) - { - std::cout << "Non-zero residual at index " << i << ": " << f[i] << "\n"; - success = false; - } + 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.raw_v7, tol_); + success *= isEqual(model.y()[11], test.expected_vss, loose_tol); + success *= isEqual(vss_node.read(), test.expected_vss, loose_tol); } - // 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__); } - /** - * @brief Residual evaluation against hand-computed answer key. - * - * Sets specific y/yp values and verifies residuals match - * pre-computed expected values. See plan for derivation. - */ TestOutcome residual() { TestStatus success = true; - PhasorDynamics::SignalNode u_node; - PhasorDynamics::SignalNode 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]) + using Params = PhasorDynamics::Stabilizer::IeeestParameters; + + struct ResidualCase + { + const char* name; + std::function edit; + const std::vector expected; }; - // 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 std::vector cases = { + {"baseline", + [](DataT&) {}, + {0.19, 0.28, 0.37, 0.0878, 0.25, 0.24, -0.05, -0.42, -0.25, -0.31, 5.75, 0.0}}, + {"a4_zero", + [](DataT& data) + { + data.parameters[Params::A4] = static_cast(0.0); + }, + {0.19, 0.28, 0.37, 0.227, 0.25, 0.24, -0.05, -0.42, -0.25, -0.31, 5.75, 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, 0.28, 0.40, 0.32, 0.25, 0.24, -0.05, -0.42, -0.25, -0.31, 5.75, 0.0}}, + {"time_zero", + [](DataT& data) + { + data.parameters[Params::T2] = static_cast(0.0); + data.parameters[Params::T4] = static_cast(0.0); + data.parameters[Params::T6] = static_cast(0.0); + }, + {0.19, 0.28, 0.37, 0.0878, 0.30, 0.30, 0.30, -0.42, -0.10, -0.10, 9.95, 0.0}}, + }; - for (size_t i = 0; i < res_answer.size(); ++i) + const auto loose_tol = static_cast(1.0e-4); + for (const auto& test : cases) { - auto test_tol = (i == 11) ? loose_tol : static_cast(10 * std::numeric_limits::epsilon()); - if (!isEqual(residual[i], res_answer[i], test_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 << "Incorrect result for residual " << i << ": " - << std::setprecision(15) << residual[i] - << " != " << res_answer[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 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(); - - std::vector - dependency_tracking_jacobian = DependencyTrackingJacobian(data); + { + PhasorDynamics::Stabilizer::Ieeest model(makeData()); + model.allocate(); + success *= (model.verify() != 0); + } - std::vector - enzyme_jacobian = EnzymeJacobian(data); + { + PhasorDynamics::SignalNode u_node; + PhasorDynamics::Stabilizer::Ieeest model(makeData()); + 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); + + 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); - - stab.allocate(); - stab.initialize(); - - // --- d/dy: tag internal variables as independent --- - for (size_t i = 0; i < stab.size(); ++i) - { - stab.y()[i].setVariableNumber(i); - } - // Tag external signal u as an additional independent variable - u_value.setVariableNumber(stab.size()); - u_value.setValue(0.5); - - setStatePointDep(stab); + TestStatus success = true; + using DepVar = DependencyTracking::Variable; - stab.evaluateResidual(); - std::vector residual_y = stab.getResidual(); + std::vector dependency_tracking_jacobian; - // --- d/dy': tag derivatives as independent --- - stab.initialize(); - for (size_t i = 0; i < stab.size(); ++i) { - stab.yp()[i].setVariableNumber(i); - } + 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_value = 0.5; - setStatePointDep(stab); + u_node.set(&u_value, &u_index); + vss_node.set(&vss_value, &vss_index); - stab.evaluateResidual(); - std::vector residual_yp = stab.getResidual(); - - // Print dependencies for debugging - for (size_t i = 0; i < residual_y.size(); ++i) - { - std::cout << i << "th residual, y: "; - (residual_y[i]).print(std::cout); - std::cout << "\n"; - std::cout << i << "th residual, yp: "; - (residual_yp[i]).print(std::cout); - std::cout << "\n"; - } + PhasorDynamics::Stabilizer::Ieeest model(makeData()); + model.getSignals().template attachSignalNode(&u_node); + model.getSignals().template assignSignalNode(&vss_node); - // 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(); + model.allocate(); + model.initialize(); - 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 enzyme_jacobian; - 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); - - 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); + { + 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); + } + + for (size_t i = 0; i < dependency_tracking_jacobian.size(); ++i) + { + success *= GridKit::Testing::isEqual(dependency_tracking_jacobian[i], enzyme_jacobian[i], tol_); + } + + return success.report(__func__); } #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 +443,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..2a76a60b0 100644 --- a/tests/UnitTests/PhasorDynamics/runStabilizerIeeestTests.cpp +++ b/tests/UnitTests/PhasorDynamics/runStabilizerIeeestTests.cpp @@ -6,9 +6,9 @@ int main() GridKit::Testing::StabilizerIeeestTests test; - result += test.constructor(); - result += test.zeroInitialResidual(); + result += test.init(); result += test.residual(); + result += test.verify(); #ifdef GRIDKIT_ENABLE_ENZYME result += test.jacobian(); #endif From 39d3a962c2e29bbc64dcdff55dfd6064806460d9 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Mon, 8 Jun 2026 07:54:37 -0500 Subject: [PATCH 2/4] Update Docs --- .../Stabilizer/IEEEST/README.md | 92 ++++++++++--------- 1 file changed, 49 insertions(+), 43 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md index d1d5a8870..a0ea1064d 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md @@ -4,7 +4,7 @@ Standard IEEE power system stabilizer: 4th-order notch filter, two lead–lag blocks, washout, and output limiter. Notes: -- $V_{cl}$, $V_{cu}$, and $T_{delay}$ are accepted for input-format +- $V_{\mathrm{cl}}$, $V_{\mathrm{cu}}$, and $T_{\mathrm{delay}}$ are accepted for input-format compatibility but are not modeled. - A zero denominator time constant bypasses its corresponding lead–lag or washout block. @@ -17,26 +17,26 @@ Figure 1: Stabilizer IEEEST model. Figure courtesy of [PowerWorld](https://www.p ## Model Parameters -Symbol | Units | JSON | Description | Typical Value | Note --------------|--------|----------|--------------------------------------|---------------|------ -$A_1$ | [s] | `A1` | Notch denominator coefficient | 1.013 | -$A_2$ | [s²] | `A2` | Notch denominator coefficient | 0.013 | -$A_3$ | [s] | `A3` | Notch denominator coefficient | 0.0 | -$A_4$ | [s²] | `A4` | Notch denominator coefficient | 0.0 | -$A_5$ | [s] | `A5` | Notch numerator coefficient | 1.013 | -$A_6$ | [s²] | `A6` | Notch numerator coefficient | 0.113 | -$T_1$ | [s] | `T1` | Lead–lag 1 numerator time constant | 0.0 | -$T_2$ | [s] | `T2` | Lead–lag 1 denominator time constant | 0.02 | -$T_3$ | [s] | `T3` | Lead–lag 2 numerator time constant | 0.0 | -$T_4$ | [s] | `T4` | Lead–lag 2 denominator time constant | 0.0 | -$T_5$ | [s] | `T5` | Washout numerator time constant | 1.65 | -$T_6$ | [s] | `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_{cl}$ | [p.u.] | `Vcl` | Lower input cutout threshold | 0.0 | Accepted but not modeled -$V_{cu}$ | [p.u.] | `Vcu` | Upper input cutout threshold | 0.0 | Accepted but not modeled -$T_\text{delay}$ | [s] | `Tdelay` | Input delay | 0.0 | Accepted but not modeled +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 @@ -85,33 +85,37 @@ The binary DAE selectors choose the active notch-filter order: #### 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 For reduced-order notch filters, unused notch states remain in the fixed component state vector and are pinned by algebraic residuals. #### 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 | Description | Note +-------|--------|-------------------------|------ +$u$ | [p.u.] | Stabilizer input signal | ## Model Equations @@ -149,17 +153,19 @@ v_5 - v_6 & T_4 = 0 -T_6 v_7 + K_s T_5(v_6 - x_7) & T_6 \ne 0 \\ K_s v_6 - v_7 & T_6 = 0 \end{cases} \\ -0 &= -V_{ss} + \text{clamp}(v_7, L_s^{\min}, L_s^{\max}) +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 -States and derivatives initialize to the steady state implied by the attached -input $u$: +Initialization is performed by evaluating the steady-state residuals in +dependency order. Let subscript $0$ denote initial values and set all internal +derivatives to zero. States and derivatives initialize to the steady state +implied by the attached input $u$: ```math \begin{aligned} @@ -171,7 +177,7 @@ K_su & T_6 = 0 \\ 0 & \text{otherwise} \end{cases} & -V_{ss} &= \text{clamp}(v_7, L_s^{\min}, L_s^{\max}) +V_{\mathrm{ss}} &= \text{clamp}(v_7, L_s^{\min}, L_s^{\max}) \end{aligned} ``` @@ -180,5 +186,5 @@ All internal derivatives initialize to zero. ## Model Outputs Output | Units | Description | Note --------|--------|---------------------------|----- +-------|--------|---------------------------|------ `vss` | [p.u.] | Limited stabilizer signal | Exported through `output` when assigned From 352e3bcd023c4cb8340c3fb8ae8ebc12610f3054 Mon Sep 17 00:00:00 2001 From: Wyatt Lowery Date: Wed, 24 Jun 2026 19:23:11 -0500 Subject: [PATCH 3/4] isEqual change and tidy --- .../Stabilizer/IEEEST/IeeestImpl.hpp | 39 ++++++++++++------- 1 file changed, 25 insertions(+), 14 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp index ec2649aa5..72b0b401c 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp @@ -13,6 +13,7 @@ #include #include #include +#include #include namespace GridKit @@ -22,6 +23,7 @@ namespace GridKit namespace Stabilizer { using Log = ::GridKit::Utilities::Logger; + using Testing::isEqual; template Ieeest::Ieeest() @@ -191,7 +193,7 @@ namespace GridKit ret += 1; } - if (a2_ == ZERO && a3_ == ZERO && a4_ == ZERO && a1_ != ZERO) + if (isEqual(a2_, ZERO) && isEqual(a3_, ZERO) && isEqual(a4_, ZERO) && !isEqual(a1_, ZERO)) { Log::error() << "Ieeest: a2, a3, and a4 are all zero - no valid notch filter\n"; ret += 1; @@ -228,7 +230,7 @@ namespace GridKit const ScalarT x6 = v5; const ScalarT v6 = v5; const ScalarT x7 = v6; - const ScalarT v7 = (T6_ == ZERO) ? Ks_ * v6 : zero; + const ScalarT v7 = isEqual(T6_, ZERO) ? Ks_ * v6 : zero; y_[0] = x1; y_[1] = x2; @@ -249,13 +251,13 @@ namespace GridKit template int Ieeest::tagDifferentiable() { - tag_[0] = (a2_ != ZERO || a3_ != ZERO || a4_ != ZERO); + tag_[0] = (!isEqual(a2_, ZERO) || !isEqual(a3_, ZERO) || !isEqual(a4_, ZERO)); tag_[1] = tag_[0]; - tag_[2] = (a3_ != ZERO || a4_ != ZERO); - tag_[3] = (a4_ != ZERO); - tag_[4] = (T2_ != ZERO); - tag_[5] = (T4_ != ZERO); - tag_[6] = (T6_ != ZERO); + tag_[2] = (!isEqual(a3_, ZERO) || !isEqual(a4_, ZERO)); + tag_[3] = (!isEqual(a4_, ZERO)); + tag_[4] = (!isEqual(T2_, ZERO)); + tag_[5] = (!isEqual(T4_, ZERO)); + tag_[6] = (!isEqual(T6_, ZERO)); tag_[7] = false; tag_[8] = false; tag_[9] = false; @@ -315,17 +317,26 @@ namespace GridKit ScalarT u = ws[0]; - f[0] = -tag_[0] * x1_dot + x2; - f[1] = -tag_[1] * x2_dot + x3; - f[2] = -tag_[2] * x3_dot + x4; + const ScalarT zero{ZERO}; + const ScalarT one{ONE}; + const ScalarT dx1 = tag_[0] ? one : zero; + const ScalarT dx2 = tag_[1] ? one : zero; + const ScalarT dx3 = tag_[2] ? one : zero; + const ScalarT t2 = tag_[4] ? one : zero; + const ScalarT t4 = tag_[5] ? one : zero; + const ScalarT t6 = tag_[6] ? one : zero; + + f[0] = -dx1 * x1_dot + x2; + f[1] = -dx2 * x2_dot + x3; + f[2] = -dx3 * x3_dot + x4; f[3] = -a4_ * x4_dot - x1 - a1_ * x2 - a2_ * x3 - a3_ * x4 + u; f[4] = -T2_ * x5_dot - x5 + v4; f[5] = -T4_ * x6_dot - x6 + v5; f[6] = -T6_ * x7_dot - x7 + v6; f[7] = -v4 + x1 + A5_ * x2 + A6_ * x3; - f[8] = tag_[4] * (-T2_ * (v5 - x5) + T1_ * (v4 - x5)) + (1 - tag_[4]) * (v4 - v5); - f[9] = tag_[5] * (-T4_ * (v6 - x6) + T3_ * (v5 - x6)) + (1 - tag_[5]) * (v5 - v6); - f[10] = tag_[6] * (-T6_ * v7 + Ks_ * T5_ * (v6 - x7)) + (1 - tag_[6]) * (Ks_ * v6 - v7); + f[8] = t2 * (-T2_ * (v5 - x5) + T1_ * (v4 - x5)) + (one - t2) * (v4 - v5); + f[9] = t4 * (-T4_ * (v6 - x6) + T3_ * (v5 - x6)) + (one - t4) * (v5 - v6); + f[10] = t6 * (-T6_ * v7 + Ks_ * T5_ * (v6 - x7)) + (one - t6) * (Ks_ * v6 - v7); f[11] = -vss + Math::clamp(v7, Lsmin_, Lsmax_); return 0; From 282fc1da68e16ef26a3504c4202800d3db508619 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Fri, 26 Jun 2026 07:43:08 -0500 Subject: [PATCH 4/4] Hessenberg form --- .../Stabilizer/IEEEST/Ieeest.hpp | 22 + .../Stabilizer/IEEEST/IeeestImpl.hpp | 408 ++++++++++-------- .../Stabilizer/IEEEST/README.md | 175 ++++---- .../PhasorDynamics/StabilizerIeeestTests.hpp | 133 +++++- .../runStabilizerIeeestTests.cpp | 1 + 5 files changed, 470 insertions(+), 269 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp index 399398ae8..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}; @@ -143,6 +152,19 @@ namespace GridKit RealT a3_{0}; RealT a4_{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_; std::unique_ptr monitor_; diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp index 72b0b401c..785af0826 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp @@ -6,14 +6,14 @@ * @brief Definition of the IEEEST Power System Stabilizer. */ -#include +#include #include +#include #include #include #include #include -#include #include namespace GridKit @@ -23,12 +23,12 @@ namespace GridKit namespace Stabilizer { using Log = ::GridKit::Utilities::Logger; - using Testing::isEqual; template Ieeest::Ieeest() { - size_ = 12; + size_ = static_cast(IeeestInternalVariables::MAXIMUM); + setDerivedParameters(); } template @@ -37,7 +37,8 @@ namespace GridKit { initializeParameters(data); initializeMonitor(); - size_ = 12; + size_ = static_cast(IeeestInternalVariables::MAXIMUM); + setDerivedParameters(); } template @@ -48,89 +49,96 @@ namespace GridKit template void Ieeest::setDerivedParameters() { + 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_; - } - template - void Ieeest::initializeParameters(const ModelDataT& data) - { - 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)) + order_ = 0; + if (a1_ != ZERO) { - Ks_ = std::get(data.parameters.at(ModelDataT::Parameters::Ks)); + order_ = 1; } - if (data.parameters.contains(ModelDataT::Parameters::Lsmin)) + if (a2_ != ZERO) { - Lsmin_ = std::get(data.parameters.at(ModelDataT::Parameters::Lsmin)); + order_ = 2; } - if (data.parameters.contains(ModelDataT::Parameters::Lsmax)) + if (a3_ != ZERO) { - Lsmax_ = std::get(data.parameters.at(ModelDataT::Parameters::Lsmax)); + order_ = 3; } - if (data.parameters.contains(ModelDataT::Parameters::Vcl)) + if (a4_ != ZERO) { - Vcl_ = std::get(data.parameters.at(ModelDataT::Parameters::Vcl)); + order_ = 4; } - if (data.parameters.contains(ModelDataT::Parameters::Vcu)) - { - Vcu_ = std::get(data.parameters.at(ModelDataT::Parameters::Vcu)); - } - if (data.parameters.contains(ModelDataT::Parameters::Tdelay)) + + 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; + + 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; + } + + 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) { - Tdelay_ = std::get(data.parameters.at(ModelDataT::Parameters::Tdelay)); - } + if (!data.parameters.contains(key)) + { + return; + } - setDerivedParameters(); + 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 @@ -143,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) { @@ -166,7 +175,8 @@ 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(); @@ -177,7 +187,7 @@ namespace GridKit template int Ieeest::verify() const { - int ret = 0; + int ret = static_cast(parameter_error_count_); if (signals_.template isAttached()) { @@ -193,9 +203,9 @@ namespace GridKit ret += 1; } - if (isEqual(a2_, ZERO) && isEqual(a3_, ZERO) && isEqual(a4_, ZERO) && !isEqual(a1_, ZERO)) + 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; } @@ -205,45 +215,45 @@ 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; } - ScalarT u{0.0}; - if (signals_.template isAttached()) - { - u = signals_.template readExternalVariable(); - ws_[0] = u; - ws_indices_[0] = signals_.template readExternalVariableIndex(); - } - - const ScalarT zero{0.0}; - const ScalarT x1 = u; - const ScalarT x2 = zero; - const ScalarT x3 = zero; - const ScalarT x4 = zero; - const ScalarT v4 = x1 + A5_ * x2 + A6_ * x3; - const ScalarT x5 = v4; - const ScalarT v5 = v4; - const ScalarT x6 = v5; - const ScalarT v6 = v5; - const ScalarT x7 = v6; - const ScalarT v7 = isEqual(T6_, ZERO) ? Ks_ * v6 : zero; - - y_[0] = x1; - y_[1] = x2; - y_[2] = x3; - y_[3] = x4; - y_[4] = x5; - y_[5] = x6; - y_[6] = x7; - y_[7] = v4; - y_[8] = v5; - y_[9] = v6; - y_[10] = v7; - y_[11] = 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); + + 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; } @@ -251,18 +261,19 @@ namespace GridKit template int Ieeest::tagDifferentiable() { - tag_[0] = (!isEqual(a2_, ZERO) || !isEqual(a3_, ZERO) || !isEqual(a4_, ZERO)); - tag_[1] = tag_[0]; - tag_[2] = (!isEqual(a3_, ZERO) || !isEqual(a4_, ZERO)); - tag_[3] = (!isEqual(a4_, ZERO)); - tag_[4] = (!isEqual(T2_, ZERO)); - tag_[5] = (!isEqual(T4_, ZERO)); - tag_[6] = (!isEqual(T6_, ZERO)); - 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; } @@ -294,50 +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]; - - const ScalarT zero{ZERO}; - const ScalarT one{ONE}; - const ScalarT dx1 = tag_[0] ? one : zero; - const ScalarT dx2 = tag_[1] ? one : zero; - const ScalarT dx3 = tag_[2] ? one : zero; - const ScalarT t2 = tag_[4] ? one : zero; - const ScalarT t4 = tag_[5] ? one : zero; - const ScalarT t6 = tag_[6] ? one : zero; - - f[0] = -dx1 * x1_dot + x2; - f[1] = -dx2 * x2_dot + x3; - f[2] = -dx3 * x3_dot + x4; - f[3] = -a4_ * x4_dot - x1 - a1_ * x2 - a2_ * x3 - a3_ * x4 + u; - f[4] = -T2_ * x5_dot - x5 + v4; - f[5] = -T4_ * x6_dot - x6 + v5; - f[6] = -T6_ * x7_dot - x7 + v6; - f[7] = -v4 + x1 + A5_ * x2 + A6_ * x3; - f[8] = t2 * (-T2_ * (v5 - x5) + T1_ * (v4 - x5)) + (one - t2) * (v4 - v5); - f[9] = t4 * (-T4_ * (v6 - x6) + T3_ * (v5 - x6)) + (one - t4) * (v5 - v6); - f[10] = t6 * (-T6_ * v7 + Ks_ * T5_ * (v6 - x7)) + (one - t6) * (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; } @@ -345,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()); @@ -366,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 a0ea1064d..cdae661d3 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md @@ -6,8 +6,6 @@ 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. -- A zero denominator time constant bypasses its corresponding lead–lag or - washout block. ## Block Diagram @@ -40,42 +38,34 @@ $T_{\mathrm{delay}}$ | [sec] | `Tdelay` | Input delay ### Parameter Validation -The fixed realization rejects a first-order-only notch denominator: +IEEEST denominator time constants are conditioned and unsupported notch forms are rejected by the following checks. Let $\epsilon_T=10^{-3}$. ```math \begin{aligned} - a_2 \ne 0 \lor a_3 \ne 0 \lor a_4 \ne 0 \lor a_1 = 0 + 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 -The notch-filter denominator expands to: - -```math -\begin{aligned} -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} -``` - -The binary DAE selectors choose the active notch-filter order: - ```math \begin{aligned} -\delta_1 &= \delta_2 = -\begin{cases} -1 & a_2 \ne 0 \lor a_3 \ne 0 \lor a_4 \ne 0 \\ -0 & \text{otherwise} -\end{cases} -\\ -\delta_3 &= -\begin{cases} -1 & a_3 \ne 0 \lor a_4 \ne 0 \\ -0 & \text{otherwise} -\end{cases} + 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} ``` @@ -92,9 +82,6 @@ $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 -For reduced-order notch filters, unused notch states remain in the fixed -component state vector and are pinned by algebraic residuals. - #### Algebraic Symbol | Units | Description | Note @@ -113,9 +100,9 @@ None. #### Algebraic -Symbol | Units | Description | Note --------|--------|-------------------------|------ -$u$ | [p.u.] | Stabilizer input signal | +Symbol | Units | Type | Description | Note +-------|--------|-------|-------------------------|------ +$u$ | [p.u.] | Known | Stabilizer input signal | ## Model Equations @@ -123,13 +110,50 @@ $u$ | [p.u.] | Stabilizer input signal | ```math \begin{aligned} -0 &= -\delta_1\dot{x}_1 + x_2 \\ -0 &= -\delta_2\dot{x}_2 + x_3 \\ -0 &= -\delta_3\dot{x}_3 + x_4 \\ -0 &= -a_4\dot{x}_4 - x_1 - a_1x_2 - a_2x_3 - a_3x_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} ``` @@ -137,23 +161,25 @@ $u$ | [p.u.] | Stabilizer input signal | ```math \begin{aligned} -0 &= -v_4 + x_1 + A_5x_2 + A_6x_3 \\ -0 &= -\begin{cases} --T_2(v_5 - x_5) + T_1(v_4 - x_5) & T_2 \ne 0 \\ -v_4 - v_5 & T_2 = 0 -\end{cases} \\ -0 &= -\begin{cases} --T_4(v_6 - x_6) + T_3(v_5 - x_6) & T_4 \ne 0 \\ -v_5 - v_6 & T_4 = 0 -\end{cases} \\ -0 &= -\begin{cases} --T_6 v_7 + K_s T_5(v_6 - x_7) & T_6 \ne 0 \\ -K_s v_6 - v_7 & T_6 = 0 -\end{cases} \\ -0 &= -V_{\mathrm{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} ``` @@ -162,28 +188,31 @@ The output limiter uses GridKit's smooth ## Initialization -Initialization is performed by evaluating the steady-state residuals in -dependency order. Let subscript $0$ denote initial values and set all internal -derivatives to zero. States and derivatives initialize to the steady state -implied by the attached input $u$: +### External Priors ```math \begin{aligned} -x_1 &= v_4 = x_5 = v_5 = x_6 = v_6 = x_7 = u \\ -x_2 &= x_3 = x_4 = 0 \\ -v_7 &= -\begin{cases} -K_su & T_6 = 0 \\ -0 & \text{otherwise} -\end{cases} -& -V_{\mathrm{ss}} &= \text{clamp}(v_7, L_s^{\min}, L_s^{\max}) + u_0 &\leftarrow \text{stabilizer input signal} \end{aligned} ``` -All internal derivatives initialize to zero. +### 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. -## Model Outputs +## Monitorable Outputs Output | Units | Description | Note -------|--------|---------------------------|------ diff --git a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp index e960ae308..82fcf27db 100644 --- a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp +++ b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp @@ -42,17 +42,15 @@ namespace GridKit ScalarT u; RealT Lsmin; RealT Lsmax; - ScalarT raw_v7; + ScalarT expected_v7; ScalarT expected_vss; }; 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, 1.0, 0.25, -1.0, 1.0, 0.25, 0.25}, - {0.0, 2.0, 0.25, -1.0, 1.0, 0.50, 0.50}, - {0.0, 4.0, 0.25, -1.0, 0.6, 1.00, 0.60}, - {5.0, 3.0, 0.25, -1.0, 1.0, 0.00, 0.00}, + {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}, }; for (const auto& test : cases) @@ -83,7 +81,7 @@ namespace GridKit success *= vss_node.linked(); success *= (vss_node.getVariableIndex() == 11); - success *= isEqual(model.y()[10], test.raw_v7, tol_); + 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); } @@ -107,31 +105,44 @@ namespace GridKit const std::vector cases = { {"baseline", [](DataT&) {}, - {0.19, 0.28, 0.37, 0.0878, 0.25, 0.24, -0.05, -0.42, -0.25, -0.31, 5.75, 0.0}}, + {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, 0.37, 0.227, 0.25, 0.24, -0.05, -0.42, -0.25, -0.31, 5.75, 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, 0.28, 0.40, 0.32, 0.25, 0.24, -0.05, -0.42, -0.25, -0.31, 5.75, 0.0}}, - {"time_zero", + {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::T2] = static_cast(0.0); - data.parameters[Params::T4] = static_cast(0.0); - data.parameters[Params::T6] = static_cast(0.0); + 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.19, 0.28, 0.37, 0.0878, 0.30, 0.30, 0.30, -0.42, -0.10, -0.10, 9.95, 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}}, }; const auto loose_tol = static_cast(1.0e-4); + for (const auto& test : cases) { PhasorDynamics::SignalNode u_node; @@ -191,6 +202,81 @@ namespace GridKit } } + { + 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__); + } + + TestOutcome tags() + { + TestStatus success = true; + + 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); + + 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 < model.tag().size(); ++i) + { + const bool expected = (i <= static_cast(PhasorDynamics::Stabilizer::IeeestInternalVariables::X7)); + if (model.tag()[i] != expected) + { + std::cout << "Incorrect differential tag at row " << i << "\n"; + success = false; + } + } + return success.report(__func__); } @@ -231,6 +317,25 @@ namespace GridKit success *= (model.verify() != 0); } + { + 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__); } diff --git a/tests/UnitTests/PhasorDynamics/runStabilizerIeeestTests.cpp b/tests/UnitTests/PhasorDynamics/runStabilizerIeeestTests.cpp index 2a76a60b0..2a8f60f89 100644 --- a/tests/UnitTests/PhasorDynamics/runStabilizerIeeestTests.cpp +++ b/tests/UnitTests/PhasorDynamics/runStabilizerIeeestTests.cpp @@ -8,6 +8,7 @@ int main() result += test.init(); result += test.residual(); + result += test.tags(); result += test.verify(); #ifdef GRIDKIT_ENABLE_ENZYME result += test.jacobian();