diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/CMakeLists.txt
index abc5bf20cf..462e2d873a 100644
--- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/CMakeLists.txt
+++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/CMakeLists.txt
@@ -19,8 +19,9 @@ if(GRIDKIT_ENABLE_ENZYME)
       ClangEnzymeFlags
     COMPILE_OPTIONS
       PRIVATE
-      -mllvm
-      -enzyme-auto-sparsity=1
+      "SHELL:-mllvm -enzyme-auto-sparsity=1"
+      # Preserve runtime-order stores for Enzyme sparsification.
+      "SHELL:-mllvm -simplifycfg-sink-common=false"
       -fno-math-errno)
 else()
   gridkit_add_library(
diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp
index 448e0144c9..8ef16e95bb 100644
--- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp
+++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp
@@ -6,6 +6,11 @@
 
 #pragma once
 
+#include <cstddef>
+#include <memory>
+#include <vector>
+
+#include <GridKit/Definitions.hpp>
 #include <GridKit/Model/PhasorDynamics/Component.hpp>
 #include <GridKit/Model/PhasorDynamics/SignalNode/SignalNodeSet.hpp>
 #include <GridKit/Model/PhasorDynamics/SignalPorts.hpp>
@@ -16,6 +21,9 @@ namespace GridKit
 {
   namespace PhasorDynamics
   {
+    template <typename scalar_type, typename index_type>
+    class SignalNode;
+
     namespace Stabilizer
     {
       /// Internal variables of `Ieeest`
@@ -44,41 +52,42 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       class Ieeest : public Component<scalar_type, index_type>
       {
-        using Component<scalar_type, index_type>::gridkit_component_id_;
+        using Component<scalar_type, index_type>::abs_tol_;
+        using Component<scalar_type, index_type>::allocated_;
         using Component<scalar_type, index_type>::alpha_;
         using Component<scalar_type, index_type>::f_;
+        using Component<scalar_type, index_type>::gridkit_component_id_;
+        using Component<scalar_type, index_type>::J_cols_buffer_;
+        using Component<scalar_type, index_type>::J_rows_buffer_;
+        using Component<scalar_type, index_type>::J_vals_buffer_;
         using Component<scalar_type, index_type>::nnz_;
+        using Component<scalar_type, index_type>::residual_indices_;
         using Component<scalar_type, index_type>::size_;
         using Component<scalar_type, index_type>::tag_;
-        using Component<scalar_type, index_type>::abs_tol_;
-        using Component<scalar_type, index_type>::time_;
+        using Component<scalar_type, index_type>::variable_indices_;
+        using Component<scalar_type, index_type>::wb_;
         using Component<scalar_type, index_type>::y_;
         using Component<scalar_type, index_type>::yp_;
-        using Component<scalar_type, index_type>::wb_;
         using Component<scalar_type, index_type>::ws_;
         using Component<scalar_type, index_type>::ws_indices_;
-        using Component<scalar_type, index_type>::h_;
-        using Component<scalar_type, index_type>::J_rows_buffer_;
-        using Component<scalar_type, index_type>::J_cols_buffer_;
-        using Component<scalar_type, index_type>::J_vals_buffer_;
-        using Component<scalar_type, index_type>::variable_indices_;
-        using Component<scalar_type, index_type>::residual_indices_;
-        using Component<scalar_type, index_type>::allocated_;
 
       public:
-        using ScalarT        = scalar_type;
-        using IdxT           = index_type;
-        using RealT          = typename Component<ScalarT, IdxT>::RealT;
-        using ModelDataT     = IeeestData<RealT, IdxT>;
-        using SignalNodeSetT = SignalNodeSet<ScalarT, IdxT>;
-        using SignalPortsT   = SignalPorts<ScalarT, ModelDataT>;
-        using MonitorT       = Model::VariableMonitor<Ieeest, IeeestData>;
+        using ScalarT            = scalar_type;
+        using IdxT               = index_type;
+        using RealT              = typename Component<ScalarT, IdxT>::RealT;
+        using SignalT            = SignalNode<ScalarT, IdxT>;
+        using ModelDataT         = IeeestData<RealT, IdxT>;
+        using SignalNodeSetT     = SignalNodeSet<ScalarT, IdxT>;
+        using SignalPortsT       = SignalPorts<ScalarT, ModelDataT>;
+        using MonitorT           = Model::VariableMonitor<Ieeest, IeeestData>;
+        using InternalVariablesT = IeeestInternalVariables;
+        using ExternalVariablesT = IeeestExternalVariables;
 
         Ieeest();
-        Ieeest(const ModelDataT& data);
+        explicit Ieeest(const ModelDataT& data);
         ~Ieeest();
 
-        int setGridKitComponentID(IdxT) override final;
+        int setGridKitComponentID(IdxT component_id) override final;
         int allocate() override final;
         int verify() const override final;
         int initialize() override final;
@@ -102,6 +111,12 @@ namespace GridKit
             ScalarT*);
 
       private:
+        void initializeParameters(const ModelDataT& data);
+        void initializeMonitor();
+        void setDerivedParameters();
+
+        static constexpr RealT TIME_CONSTANT_MINIMUM = static_cast<RealT>(1.0e-3);
+
         RealT A1_{0};
         RealT A2_{0};
         RealT A3_{0};
@@ -121,34 +136,17 @@ namespace GridKit
         RealT Vcu_{0};
         RealT Tdelay_{0};
 
-        RealT a0_{1};
+        IdxT order_{0};
+
         RealT a1_{0};
         RealT a2_{0};
         RealT a3_{0};
         RealT a4_{0};
 
-        // Precomputed masks and safe inverse coefficients for branch-free degenerate paths.
-        RealT use_notch_{0};
-        RealT bypass_notch_{1};
-        RealT use_4th_order_{0};
-        RealT use_3rd_order_{0};
-        RealT use_2nd_order_{0};
-        RealT safe_inv_a4_{0};
-        RealT safe_inv_a3_{0};
-        RealT safe_inv_a2_{0};
-        RealT use_T2_block_{1};
-        RealT bypass_T2_block_{0};
-        RealT use_T4_block_{1};
-        RealT bypass_T4_block_{0};
-        RealT use_T6_block_{1};
-        RealT bypass_T6_block_{0};
-
-        SignalPortsT ports_;
+        IdxT parameter_error_count_{0};
 
+        SignalPortsT              ports_;
         std::unique_ptr<MonitorT> monitor_;
-
-        void initializeParameters(const ModelDataT& data);
-        void initializeMonitor();
       };
 
     } // namespace Stabilizer
diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp
index 4a421a8475..17baad8be2 100644
--- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp
+++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp
@@ -14,21 +14,11 @@ namespace GridKit
   {
     namespace Stabilizer
     {
-      /**
-       * @brief Jacobian evaluation experimental
-       *
-       * @tparam ScalarT - Scalar data type
-       * @tparam IdxT    - Index data type
-       * @return int - error code, 0 = success
-       */
       template <typename scalar_type, typename index_type>
       int Ieeest<scalar_type, index_type>::evaluateJacobian()
       {
         if (J_rows_buffer_ == nullptr)
         {
-          // Reserve space for the dense blocks.
-          // The size of the buffer is the sum of maximum capacities of the blocks.
-          // Enzyme will compute the appropriate nnz from sparsification.
           auto size        = static_cast<size_t>(size_);
           auto signal_size = static_cast<size_t>(ws_.getSize());
           auto buffer_size = 2 * size * size + size * signal_size;
@@ -37,63 +27,61 @@ namespace GridKit
           J_vals_buffer_   = new RealT[buffer_size];
         }
 
+        using ModelT = GridKit::PhasorDynamics::Stabilizer::Ieeest<scalar_type, index_type>;
+        using Fn     = GridKit::Enzyme::Sparse::MemberFunctions;
+
         nnz_ = 0;
 
-        GridKit::Enzyme::Sparse::DfDy<GridKit::PhasorDynamics::Stabilizer::Ieeest<ScalarT, IdxT>,
-                                      GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this,
-                                                                                                                  static_cast<size_t>(f_.getSize()),
-                                                                                                                  static_cast<size_t>(y_.getSize()),
-                                                                                                                  (this->getResidualIndices()).data(),
-                                                                                                                  (this->getVariableIndices()).data(),
-                                                                                                                  y_.getData(),
-                                                                                                                  yp_.getData(),
-                                                                                                                  nullptr,
-                                                                                                                  ws_.getData(),
-                                                                                                                  J_rows_buffer_,
-                                                                                                                  J_cols_buffer_,
-                                                                                                                  J_vals_buffer_,
-                                                                                                                  nnz_);
+        GridKit::Enzyme::Sparse::DfDy<ModelT, Fn::InternalResidualWithSignal>::eval(this,
+                                                                                    static_cast<size_t>(f_.getSize()),
+                                                                                    static_cast<size_t>(y_.getSize()),
+                                                                                    (this->getResidualIndices()).data(),
+                                                                                    (this->getVariableIndices()).data(),
+                                                                                    y_.getData(),
+                                                                                    yp_.getData(),
+                                                                                    nullptr,
+                                                                                    ws_.getData(),
+                                                                                    J_rows_buffer_,
+                                                                                    J_cols_buffer_,
+                                                                                    J_vals_buffer_,
+                                                                                    nnz_);
 
-        GridKit::Enzyme::Sparse::DfDyp<GridKit::PhasorDynamics::Stabilizer::Ieeest<ScalarT, IdxT>,
-                                       GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this,
-                                                                                                                   static_cast<size_t>(f_.getSize()),
-                                                                                                                   static_cast<size_t>(y_.getSize()),
-                                                                                                                   (this->getResidualIndices()).data(),
-                                                                                                                   (this->getVariableIndices()).data(),
-                                                                                                                   y_.getData(),
-                                                                                                                   yp_.getData(),
-                                                                                                                   nullptr,
-                                                                                                                   ws_.getData(),
-                                                                                                                   alpha_,
-                                                                                                                   J_rows_buffer_,
-                                                                                                                   J_cols_buffer_,
-                                                                                                                   J_vals_buffer_,
-                                                                                                                   nnz_);
+        GridKit::Enzyme::Sparse::DfDyp<ModelT, Fn::InternalResidualWithSignal>::eval(this,
+                                                                                     static_cast<size_t>(f_.getSize()),
+                                                                                     static_cast<size_t>(y_.getSize()),
+                                                                                     (this->getResidualIndices()).data(),
+                                                                                     (this->getVariableIndices()).data(),
+                                                                                     y_.getData(),
+                                                                                     yp_.getData(),
+                                                                                     nullptr,
+                                                                                     ws_.getData(),
+                                                                                     alpha_,
+                                                                                     J_rows_buffer_,
+                                                                                     J_cols_buffer_,
+                                                                                     J_vals_buffer_,
+                                                                                     nnz_);
 
-        GridKit::Enzyme::Sparse::DfDws<GridKit::PhasorDynamics::Stabilizer::Ieeest<ScalarT, IdxT>,
-                                       GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::eval(this,
-                                                                                                                   static_cast<size_t>(f_.getSize()),
-                                                                                                                   static_cast<size_t>(ws_.getSize()),
-                                                                                                                   (this->getResidualIndices()).data(),
-                                                                                                                   ws_indices_.data(),
-                                                                                                                   y_.getData(),
-                                                                                                                   yp_.getData(),
-                                                                                                                   nullptr,
-                                                                                                                   ws_.getData(),
-                                                                                                                   J_rows_buffer_,
-                                                                                                                   J_cols_buffer_,
-                                                                                                                   J_vals_buffer_,
-                                                                                                                   nnz_);
+        GridKit::Enzyme::Sparse::DfDws<ModelT, Fn::InternalResidualWithSignal>::eval(this,
+                                                                                     static_cast<size_t>(f_.getSize()),
+                                                                                     static_cast<size_t>(ws_.getSize()),
+                                                                                     (this->getResidualIndices()).data(),
+                                                                                     ws_indices_.data(),
+                                                                                     y_.getData(),
+                                                                                     yp_.getData(),
+                                                                                     nullptr,
+                                                                                     ws_.getData(),
+                                                                                     J_rows_buffer_,
+                                                                                     J_cols_buffer_,
+                                                                                     J_vals_buffer_,
+                                                                                     nnz_);
 
         this->constructCoo();
 
         return 0;
       }
 
-      // Available template instantiations
       template class Ieeest<double, long int>;
       template class Ieeest<double, size_t>;
-
     } // namespace Stabilizer
   } // namespace PhasorDynamics
 } // namespace GridKit
diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp
index 2bc836363e..6052ee6ee0 100644
--- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp
+++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp
@@ -6,10 +6,11 @@
  * @brief Definition of the IEEEST Power System Stabilizer.
  */
 
-#include <iostream>
+#include <algorithm>
+#include <cmath>
+#include <variant>
 
 #include <GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp>
-#include <GridKit/Model/PhasorDynamics/SignalNode/SignalNodeSet.hpp>
 #include <GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp>
 #include <GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestData.hpp>
 #include <GridKit/Model/VariableMonitorImpl.hpp>
@@ -26,7 +27,8 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       Ieeest<scalar_type, index_type>::Ieeest()
       {
-        size_ = 12;
+        size_ = static_cast<IdxT>(Utilities::enum_size<IeeestInternalVariables>());
+        setDerivedParameters();
       }
 
       template <typename scalar_type, typename index_type>
@@ -35,7 +37,7 @@ namespace GridKit
       {
         initializeParameters(data);
         initializeMonitor();
-        size_ = 12;
+        size_ = static_cast<IdxT>(Utilities::enum_size<IeeestInternalVariables>());
       }
 
       template <typename scalar_type, typename index_type>
@@ -43,110 +45,6 @@ namespace GridKit
       {
       }
 
-      template <typename scalar_type, typename index_type>
-      void Ieeest<scalar_type, index_type>::initializeParameters(const ModelDataT& data)
-      {
-        using Parameter = typename ModelDataT::Parameters;
-        if (data.parameters.contains(Parameter::A1))
-        {
-          A1_ = std::get<RealT>(data.parameters.at(Parameter::A1));
-        }
-        if (data.parameters.contains(Parameter::A2))
-        {
-          A2_ = std::get<RealT>(data.parameters.at(Parameter::A2));
-        }
-        if (data.parameters.contains(Parameter::A3))
-        {
-          A3_ = std::get<RealT>(data.parameters.at(Parameter::A3));
-        }
-        if (data.parameters.contains(Parameter::A4))
-        {
-          A4_ = std::get<RealT>(data.parameters.at(Parameter::A4));
-        }
-        if (data.parameters.contains(Parameter::A5))
-        {
-          A5_ = std::get<RealT>(data.parameters.at(Parameter::A5));
-        }
-        if (data.parameters.contains(Parameter::A6))
-        {
-          A6_ = std::get<RealT>(data.parameters.at(Parameter::A6));
-        }
-        if (data.parameters.contains(Parameter::T1))
-        {
-          T1_ = std::get<RealT>(data.parameters.at(Parameter::T1));
-        }
-        if (data.parameters.contains(Parameter::T2))
-        {
-          T2_ = std::get<RealT>(data.parameters.at(Parameter::T2));
-        }
-        if (data.parameters.contains(Parameter::T3))
-        {
-          T3_ = std::get<RealT>(data.parameters.at(Parameter::T3));
-        }
-        if (data.parameters.contains(Parameter::T4))
-        {
-          T4_ = std::get<RealT>(data.parameters.at(Parameter::T4));
-        }
-        if (data.parameters.contains(Parameter::T5))
-        {
-          T5_ = std::get<RealT>(data.parameters.at(Parameter::T5));
-        }
-        if (data.parameters.contains(Parameter::T6))
-        {
-          T6_ = std::get<RealT>(data.parameters.at(Parameter::T6));
-        }
-        if (data.parameters.contains(Parameter::Ks))
-        {
-          Ks_ = std::get<RealT>(data.parameters.at(Parameter::Ks));
-        }
-        if (data.parameters.contains(Parameter::Lsmin))
-        {
-          Lsmin_ = std::get<RealT>(data.parameters.at(Parameter::Lsmin));
-        }
-        if (data.parameters.contains(Parameter::Lsmax))
-        {
-          Lsmax_ = std::get<RealT>(data.parameters.at(Parameter::Lsmax));
-        }
-        if (data.parameters.contains(Parameter::Vcl))
-        {
-          Vcl_ = std::get<RealT>(data.parameters.at(Parameter::Vcl));
-        }
-        if (data.parameters.contains(Parameter::Vcu))
-        {
-          Vcu_ = std::get<RealT>(data.parameters.at(Parameter::Vcu));
-        }
-        if (data.parameters.contains(Parameter::Tdelay))
-        {
-          Tdelay_ = std::get<RealT>(data.parameters.at(Parameter::Tdelay));
-        }
-
-        a0_ = 1;
-        a1_ = A1_ + A3_;
-        a2_ = A2_ + A4_ + A1_ * A3_;
-        a3_ = A1_ * A4_ + A2_ * A3_;
-        a4_ = A2_ * A4_;
-
-        // Precompute masks and safe inverse coefficients so the residual stays branch-free.
-        use_notch_    = static_cast<RealT>(a2_ != 0.0 || a3_ != 0.0 || a4_ != 0.0);
-        bypass_notch_ = 1.0 - use_notch_;
-
-        use_4th_order_ = static_cast<RealT>(a4_ != 0.0);
-        use_3rd_order_ = static_cast<RealT>(a4_ == 0.0 && a3_ != 0.0);
-        use_2nd_order_ = static_cast<RealT>(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<RealT>(T2_ != 0.0);
-        bypass_T2_block_ = 1.0 - use_T2_block_;
-
-        use_T4_block_    = static_cast<RealT>(T4_ != 0.0);
-        bypass_T4_block_ = 1.0 - use_T4_block_;
-
-        use_T6_block_    = static_cast<RealT>(T6_ != 0.0);
-        bypass_T6_block_ = 1.0 - use_T6_block_;
-      }
-
       template <typename scalar_type, typename index_type>
       int Ieeest<scalar_type, index_type>::setGridKitComponentID(IdxT component_id)
       {
@@ -161,26 +59,29 @@ namespace GridKit
         {
           this->allocateVectors(size_);
         }
-        auto size = static_cast<size_t>(size_);
-
-        tag_.resize(size);
+        const auto size = static_cast<size_t>(size_);
 
+        tag_.assign(size, false);
         variable_indices_.resize(size);
         residual_indices_.resize(size);
+
+        const auto signal_size = Utilities::enum_size<IeeestExternalVariables>();
+        ws_.resize(static_cast<IdxT>(signal_size));
+        ws_.setToZero();
+        ws_indices_.assign(signal_size, INVALID_INDEX<IdxT>);
+
         for (IdxT j = 0; j < size_; ++j)
         {
           this->setVariableIndex(j, j);
           this->setResidualIndex(j, j);
         }
 
-        ws_.resize(1);
-        ws_.setToZero();
-        ws_indices_.resize(1);
-        ws_indices_[0] = INVALID_INDEX<IdxT>;
-
         if (auto output_port = ports_.out.template port<IeeestSignalOutputs::output>())
         {
-          output_port.link(&y_.getData()[11], &(this->getVariableIndex(11)));
+          auto* y = y_.getData();
+          output_port.link(
+              &y[static_cast<size_t>(IeeestInternalVariables::VSS)],
+              &(this->getVariableIndex(static_cast<IdxT>(IeeestInternalVariables::VSS))));
         }
 
         allocated_ = true;
@@ -190,23 +91,73 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       int Ieeest<scalar_type, index_type>::verify() const
       {
-        int ret = 0;
+        int ret = static_cast<int>(parameter_error_count_);
 
-        auto input_port = ports_.in.template port<IeeestSignalInputs::input>();
-        if (!input_port.connected())
+        auto check = [&](bool condition, const char* message)
+        {
+          if (!condition)
+          {
+            Log::error() << "Ieeest: " << message << '\n';
+            ret += 1;
+          }
+        };
+
+        if (!ports_.in.template port<IeeestSignalInputs::input>().connected())
         {
           Log::error() << "Ieeest: required input signal U is not attached\n";
           ret += 1;
         }
-        if (input_port.connected() && !input_port.linked())
+        else if (!ports_.in.template port<IeeestSignalInputs::input>().linked())
         {
           Log::error() << "Ieeest: input signal U attached with no linked source\n";
           ret += 1;
         }
 
-        if (a4_ == 0 && a3_ == 0 && a2_ == 0 && a1_ != 0)
+        check(std::isfinite(a1_) && std::isfinite(a2_)
+                  && std::isfinite(a3_) && std::isfinite(a4_),
+              "expanded denominator coefficients must be finite");
+
+        IdxT numerator_order = static_cast<IdxT>(0);
+        if (A6_ != ZERO<RealT>)
         {
-          Log::error() << "Ieeest: a2, a3, and a4 are all zero - no valid notch filter\n";
+          numerator_order = static_cast<IdxT>(2);
+        }
+        else if (A5_ != ZERO<RealT>)
+        {
+          numerator_order = static_cast<IdxT>(1);
+        }
+        check(numerator_order <= order_,
+              "numerator order must not exceed denominator order");
+        check(Lsmin_ < Lsmax_, "Lsmin must be less than Lsmax");
+
+        RealT       leading_coefficient = ONE<RealT>;
+        const char* leading_name        = nullptr;
+        switch (order_)
+        {
+        case 0:
+          break;
+        case 1:
+          leading_coefficient = a1_;
+          leading_name        = "a1";
+          break;
+        case 2:
+          leading_coefficient = a2_;
+          leading_name        = "a2";
+          break;
+        case 3:
+          leading_coefficient = a3_;
+          leading_name        = "a3";
+          break;
+        case 4:
+          leading_coefficient = a4_;
+          leading_name        = "a4";
+          break;
+        }
+
+        if (leading_name != nullptr && leading_coefficient == ZERO<RealT>)
+        {
+          Log::error() << "Ieeest: active leading denominator coefficient '"
+                       << leading_name << "' must be nonzero\n";
           ret += 1;
         }
 
@@ -216,38 +167,59 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       int Ieeest<scalar_type, index_type>::initialize()
       {
+        const auto X1  = static_cast<size_t>(IeeestInternalVariables::X1);
+        const auto X5  = static_cast<size_t>(IeeestInternalVariables::X5);
+        const auto X6  = static_cast<size_t>(IeeestInternalVariables::X6);
+        const auto X7  = static_cast<size_t>(IeeestInternalVariables::X7);
+        const auto V4  = static_cast<size_t>(IeeestInternalVariables::V4);
+        const auto V5  = static_cast<size_t>(IeeestInternalVariables::V5);
+        const auto V6  = static_cast<size_t>(IeeestInternalVariables::V6);
+        const auto V7  = static_cast<size_t>(IeeestInternalVariables::V7);
+        const auto VSS = static_cast<size_t>(IeeestInternalVariables::VSS);
+        const auto U   = static_cast<size_t>(IeeestExternalVariables::U);
+
         if (verify() != 0)
         {
           Log::error() << "Ieeest: cannot initialize with invalid configuration\n";
           return 1;
         }
 
-        auto          input_port = ports_.in.template port<IeeestSignalInputs::input>();
-        const ScalarT u          = input_port.readSignal();
+        const ScalarT u = ports_.in.template port<IeeestSignalInputs::input>().readSignal();
+        if (!std::isfinite(static_cast<RealT>(u)))
+        {
+          Log::error() << "Ieeest: initial input signal U must be finite\n";
+          return 1;
+        }
 
         auto* y  = y_.getData();
         auto* yp = yp_.getData();
+        auto* ws = ws_.getData();
+        std::fill_n(y, static_cast<size_t>(size_), ScalarT{ZERO<RealT>});
+        std::fill_n(yp, static_cast<size_t>(size_), ScalarT{ZERO<RealT>});
 
-        for (IdxT i = 0; i < size_; ++i)
+        ws[U]          = u;
+        ws_indices_[U] = ports_.in.template port<IeeestSignalInputs::input>().signalVariableIndex();
+
+        switch (order_)
         {
-          y[static_cast<size_t>(i)]  = 0.0;
-          yp[static_cast<size_t>(i)] = 0.0;
+        case 0:
+          break;
+        case 1:
+        case 2:
+        case 3:
+        case 4:
+          y[X1] = u;
+          break;
         }
 
-        ws_.getData()[0] = u;
-        ws_indices_[0]   = input_port.signalVariableIndex();
-
-        y[0] = use_notch_ * u;
-        y[4] = u;
-        y[5] = u;
-        y[6] = u;
-        y[7] = u;
-        y[8] = u;
-        y[9] = u;
-
-        // Preserve the current T6 = 0 bypass behavior.
-        y[10] = bypass_T6_block_ * Ks_ * u;
-        y[11] = Math::clamp(y[10], Lsmin_, Lsmax_);
+        y[X5]  = u;
+        y[X6]  = u;
+        y[X7]  = u;
+        y[V4]  = u;
+        y[V5]  = u;
+        y[V6]  = u;
+        y[V7]  = ZERO<RealT>;
+        y[VSS] = Math::clamp(y[V7], Lsmin_, Lsmax_);
 
         // For DependencyTracking::Variable, set variable numbers
         if constexpr (std::is_same_v<scalar_type, DependencyTracking::Variable>)
@@ -264,18 +236,9 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       int Ieeest<scalar_type, index_type>::tagDifferentiable()
       {
-        tag_[0]  = true;
-        tag_[1]  = true;
-        tag_[2]  = true;
-        tag_[3]  = true;
-        tag_[4]  = (T2_ != 0.0);
-        tag_[5]  = (T4_ != 0.0);
-        tag_[6]  = (T6_ != 0.0);
-        tag_[7]  = false;
-        tag_[8]  = false;
-        tag_[9]  = false;
-        tag_[10] = false;
-        tag_[11] = false;
+        std::fill(tag_.begin(), tag_.end(), false);
+        constexpr auto X7 = static_cast<size_t>(IeeestInternalVariables::X7);
+        std::fill_n(tag_.begin(), X7 + 1, true);
 
         return 0;
       }
@@ -307,45 +270,106 @@ namespace GridKit
           const ScalarT*                  ws,
           ScalarT*                        f)
       {
-        ScalarT x1  = y[0];
-        ScalarT x2  = y[1];
-        ScalarT x3  = y[2];
-        ScalarT x4  = y[3];
-        ScalarT x5  = y[4];
-        ScalarT x6  = y[5];
-        ScalarT x7  = y[6];
-        ScalarT v4  = y[7];
-        ScalarT v5  = y[8];
-        ScalarT v6  = y[9];
-        ScalarT v7  = y[10];
-        ScalarT vss = y[11];
-
-        ScalarT x1_dot = yp[0];
-        ScalarT x2_dot = yp[1];
-        ScalarT x3_dot = yp[2];
-        ScalarT x4_dot = yp[3];
-        ScalarT x5_dot = yp[4];
-        ScalarT x6_dot = yp[5];
-        ScalarT x7_dot = yp[6];
-
-        ScalarT u = ws[0];
-
-        const ScalarT x2_rhs = (use_4th_order_ + use_3rd_order_) * x3
-                               + use_2nd_order_ * (-a0_ * x1 - a1_ * x2 + u) * safe_inv_a2_;
-
-        f[0] = -x1_dot + use_notch_ * x2;
-        f[1] = -x2_dot + x2_rhs;
-        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 + A6_ * x2_rhs);
-        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<size_t>(IeeestInternalVariables::X1);
+        const auto X2  = static_cast<size_t>(IeeestInternalVariables::X2);
+        const auto X3  = static_cast<size_t>(IeeestInternalVariables::X3);
+        const auto X4  = static_cast<size_t>(IeeestInternalVariables::X4);
+        const auto X5  = static_cast<size_t>(IeeestInternalVariables::X5);
+        const auto X6  = static_cast<size_t>(IeeestInternalVariables::X6);
+        const auto X7  = static_cast<size_t>(IeeestInternalVariables::X7);
+        const auto V4  = static_cast<size_t>(IeeestInternalVariables::V4);
+        const auto V5  = static_cast<size_t>(IeeestInternalVariables::V5);
+        const auto V6  = static_cast<size_t>(IeeestInternalVariables::V6);
+        const auto V7  = static_cast<size_t>(IeeestInternalVariables::V7);
+        const auto VSS = static_cast<size_t>(IeeestInternalVariables::VSS);
+        const auto U   = static_cast<size_t>(IeeestExternalVariables::U);
+
+        const ScalarT x1  = y[X1];
+        const ScalarT x2  = y[X2];
+        const ScalarT x3  = y[X3];
+        const ScalarT x4  = y[X4];
+        const ScalarT x5  = y[X5];
+        const ScalarT x6  = y[X6];
+        const ScalarT x7  = y[X7];
+        const ScalarT v4  = y[V4];
+        const ScalarT v5  = y[V5];
+        const ScalarT v6  = y[V6];
+        const ScalarT v7  = y[V7];
+        const ScalarT vss = y[VSS];
+
+        const ScalarT x1_dot = yp[X1];
+        const ScalarT x2_dot = yp[X2];
+        const ScalarT x3_dot = yp[X3];
+        const ScalarT x4_dot = yp[X4];
+        const ScalarT x5_dot = yp[X5];
+        const ScalarT x6_dot = yp[X6];
+        const ScalarT x7_dot = yp[X7];
+
+        const ScalarT u = ws[U];
+
+        switch (order_)
+        {
+        case 0:
+          f[X1] = -x1_dot;
+          f[X2] = -x2_dot;
+          f[X3] = -x3_dot;
+          f[X4] = -x4_dot;
+          f[V4] = -v4 + u;
+          break;
+        case 1:
+        {
+          const ScalarT x1_rhs = (u - x1) / a1_;
+          f[X1]                = -x1_dot + x1_rhs;
+          f[X2]                = -x2_dot;
+          f[X3]                = -x3_dot;
+          f[X4]                = -x4_dot;
+          f[V4]                = -v4 + x1 + A5_ * x1_rhs;
+          break;
+        }
+        case 2:
+        {
+          const ScalarT x2_rhs = (u - x1 - a1_ * x2) / a2_;
+          f[X1]                = -x1_dot + x2;
+          f[X2]                = -x2_dot + x2_rhs;
+          f[X3]                = -x3_dot;
+          f[X4]                = -x4_dot;
+          f[V4]                = -v4 + x1 + A5_ * x2 + A6_ * x2_rhs;
+          break;
+        }
+        case 3:
+        {
+          const ScalarT x3_rhs = (u - x1 - a1_ * x2 - a2_ * x3) / a3_;
+          f[X1]                = -x1_dot + x2;
+          f[X2]                = -x2_dot + x3;
+          f[X3]                = -x3_dot + x3_rhs;
+          f[X4]                = -x4_dot;
+          f[V4]                = -v4 + x1 + A5_ * x2 + A6_ * x3;
+          break;
+        }
+        case 4:
+        {
+          const ScalarT x4_rhs =
+              (u - x1 - a1_ * x2 - a2_ * x3 - a3_ * x4) / a4_;
+          f[X1] = -x1_dot + x2;
+          f[X2] = -x2_dot + x3;
+          f[X3] = -x3_dot + x4;
+          f[X4] = -x4_dot + x4_rhs;
+          f[V4] = -v4 + x1 + A5_ * x2 + A6_ * x3;
+          break;
+        }
+        }
+
+        const ScalarT x5_rhs = (v4 - x5) / T2_;
+        const ScalarT x6_rhs = (v5 - x6) / T4_;
+        const ScalarT x7_rhs = (v6 - x7) / T6_;
+
+        f[X5]  = -x5_dot + x5_rhs;
+        f[X6]  = -x6_dot + x6_rhs;
+        f[X7]  = -x7_dot + x7_rhs;
+        f[V5]  = -v5 + x5 + T1_ * x5_rhs;
+        f[V6]  = -v6 + x6 + T3_ * x6_rhs;
+        f[V7]  = -v7 + Ks_ * T5_ * x7_rhs;
+        f[VSS] = -vss + Math::clamp(v7, Lsmin_, Lsmax_);
 
         return 0;
       }
@@ -353,13 +377,11 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       int Ieeest<scalar_type, index_type>::evaluateResidual()
       {
-        auto* ws = ws_.getData();
+        auto*      ws = ws_.getData();
+        const auto U  = static_cast<size_t>(IeeestExternalVariables::U);
 
-        if (auto input_port = ports_.in.template port<IeeestSignalInputs::input>())
-        {
-          ws[0]          = input_port.readSignal();
-          ws_indices_[0] = input_port.signalVariableIndex();
-        }
+        ws[U]          = ports_.in.template port<IeeestSignalInputs::input>().readSignal();
+        ws_indices_[U] = ports_.in.template port<IeeestSignalInputs::input>().signalVariableIndex();
 
         const auto* y  = y_.getData();
         const auto* yp = yp_.getData();
@@ -376,12 +398,145 @@ namespace GridKit
         return monitor_.get();
       }
 
+      //
+      //  Private methods
+      //
+
+      template <typename scalar_type, typename index_type>
+      void Ieeest<scalar_type, index_type>::initializeParameters(const ModelDataT& data)
+      {
+        using Params = typename ModelDataT::Parameters;
+
+        parameter_error_count_ = 0;
+
+        auto load_real = [&](auto key, RealT& target, const char* name) -> bool
+        {
+          if (!data.parameters.contains(key))
+          {
+            return false;
+          }
+
+          const auto& value = data.parameters.at(key);
+          RealT       parsed_value{};
+          if (const auto* real_value = std::get_if<RealT>(&value))
+          {
+            parsed_value = *real_value;
+          }
+          else if (const auto* index_value = std::get_if<IdxT>(&value))
+          {
+            parsed_value = static_cast<RealT>(*index_value);
+          }
+          else
+          {
+            Log::error() << "Ieeest: parameter '" << name << "' must be numeric\n";
+            ++parameter_error_count_;
+            return false;
+          }
+
+          if (!std::isfinite(parsed_value))
+          {
+            Log::error() << "Ieeest: parameter '" << name << "' must be finite\n";
+            ++parameter_error_count_;
+            return false;
+          }
+
+          target = parsed_value;
+          return true;
+        };
+
+        load_real(Params::A1, A1_, "A1");
+        load_real(Params::A2, A2_, "A2");
+        load_real(Params::A3, A3_, "A3");
+        load_real(Params::A4, A4_, "A4");
+        load_real(Params::A5, A5_, "A5");
+        load_real(Params::A6, A6_, "A6");
+        load_real(Params::T1, T1_, "T1");
+        load_real(Params::T2, T2_, "T2");
+        load_real(Params::T3, T3_, "T3");
+        load_real(Params::T4, T4_, "T4");
+        load_real(Params::T5, T5_, "T5");
+        load_real(Params::T6, T6_, "T6");
+        load_real(Params::Ks, Ks_, "Ks");
+        load_real(Params::Lsmin, Lsmin_, "Lsmin");
+        load_real(Params::Lsmax, Lsmax_, "Lsmax");
+        if (load_real(Params::Vcl, Vcl_, "Vcl") && Vcl_ != ZERO<RealT>)
+        {
+          Log::warning() << "Ieeest: nonzero Vcl requests a lower input cutout, "
+                            "but input cutout is not implemented and Vcl is ignored\n";
+        }
+        if (load_real(Params::Vcu, Vcu_, "Vcu") && Vcu_ != ZERO<RealT>)
+        {
+          Log::warning() << "Ieeest: nonzero Vcu requests an upper input cutout, "
+                            "but input cutout is not implemented and Vcu is ignored\n";
+        }
+        if (load_real(Params::Tdelay, Tdelay_, "Tdelay") && Tdelay_ != ZERO<RealT>)
+        {
+          Log::warning() << "Ieeest: nonzero Tdelay requests an input delay, "
+                            "but input delay is not implemented and Tdelay is ignored\n";
+        }
+
+        setDerivedParameters();
+      }
+
       template <typename scalar_type, typename index_type>
       void Ieeest<scalar_type, index_type>::initializeMonitor()
       {
         using Variable = typename ModelDataT::MonitorableVariables;
+
+        constexpr auto VSS = static_cast<size_t>(IeeestInternalVariables::VSS);
+
         monitor_->set(Variable::vss, [this]
-                      { return y_.getData()[11]; });
+                      { return y_.getData()[VSS]; });
+      }
+
+      template <typename scalar_type, typename index_type>
+      void Ieeest<scalar_type, index_type>::setDerivedParameters()
+      {
+        auto check_non_negative = [&](RealT value, const char* name)
+        {
+          if (value < ZERO<RealT>)
+          {
+            Log::error() << "Ieeest: " << name << " must be non-negative\n";
+            ++parameter_error_count_;
+          }
+        };
+
+        check_non_negative(T2_, "T2");
+        check_non_negative(T4_, "T4");
+        check_non_negative(T6_, "T6");
+
+        if (T2_ < TIME_CONSTANT_MINIMUM
+            || T4_ < TIME_CONSTANT_MINIMUM
+            || T6_ < TIME_CONSTANT_MINIMUM)
+        {
+          Log::warning() << "Ieeest: T2, T4, and T6 below "
+                         << TIME_CONSTANT_MINIMUM
+                         << " s are raised to preserve Hessenberg form\n";
+        }
+
+        T2_ = std::max(T2_, TIME_CONSTANT_MINIMUM);
+        T4_ = std::max(T4_, TIME_CONSTANT_MINIMUM);
+        T6_ = std::max(T6_, TIME_CONSTANT_MINIMUM);
+
+        a1_ = A1_ + A3_;
+        a2_ = A2_ + A4_ + A1_ * A3_;
+        a3_ = A1_ * A4_ + A2_ * A3_;
+        a4_ = A2_ * A4_;
+
+        auto factor_order = [](RealT linear, RealT quadratic) -> IdxT
+        {
+          if (quadratic != ZERO<RealT>)
+          {
+            return static_cast<IdxT>(2);
+          }
+          if (linear != ZERO<RealT>)
+          {
+            return static_cast<IdxT>(1);
+          }
+          return static_cast<IdxT>(0);
+        };
+
+        order_ = factor_order(A1_, A2_) + factor_order(A3_, A4_);
       }
 
     } // namespace Stabilizer
diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md
index c023a495df..25886126f5 100644
--- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md
+++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md
@@ -1,75 +1,115 @@
-# IEEEST
+# **IEEE Stabilizer Model (IEEEST)**
 
-Standard IEEE power system stabilizer: notch filter of order up to four, two lead–lag
-blocks, washout, and output limiter.
+IEEEST is a standard IEEE power system stabilizer with a derived-order notch
+filter, two lead-lag blocks, washout, and an output limiter.
+
+## Notes
+
+> [!NOTE]
+> $V_{\mathrm{cl}}$, $V_{\mathrm{cu}}$, and $T_{\mathrm{delay}}$ are accepted for
+> input-format compatibility but are not modeled. A finite nonzero value logs a
+> warning and is ignored.
+
+- The notch-filter order is derived from the two configured denominator
+  factors. It is not an input parameter.
+- The numerator must be proper: its order cannot exceed the derived denominator
+  order.
 
 ## Block Diagram
 
-![](../../../../../docs/Figures/stabilizer_ieeest_diagram.png)
+![](../../../../../docs/Figures/PhasorDynamics/IEEEST/diagram.png)
 
-Figure 1: Stabilizer IEEEST model. Figure courtesy of [PowerWorld](https://www.powerworld.com/WebHelp/)
+Figure 1: IEEEST block diagram. Figure courtesy of [PowerWorld](https://www.powerworld.com/WebHelp/)
 
 ## Model Parameters
 
-Symbol       | Units  | 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_\mathrm{s}$        | [p.u.] | `Ks`    | Stabilizer gain                      | 3.0           |
-$L_\mathrm{s}^{\min}$ | [p.u.] | `Lsmin` | Minimum stabilizer output limit      | -0.1          |
-$L_\mathrm{s}^{\max}$ | [p.u.] | `Lsmax` | Maximum stabilizer output limit      | 0.1           |
-
-The IEEE 421.5 IEEEST also defines a cutout window ($V_\mathrm{cl}$, $V_\mathrm{cu}$) and an
-input delay ($T_\mathrm{delay}$). These parameters are accepted for input-format
-compatibility but are not modeled here.
+Symbol                  | Units    | JSON     | Description                          | Typical Value | Note
+------------------------|----------|----------|--------------------------------------|---------------|------
+$A_1$                   | [sec]    | `A1`     | Notch denominator coefficient        | 1.013         |
+$A_2$                   | [sec²]   | `A2`     | Notch denominator coefficient        | 0.013         |
+$A_3$                   | [sec]    | `A3`     | Notch denominator coefficient        | 0.0           |
+$A_4$                   | [sec²]   | `A4`     | Notch denominator coefficient        | 0.0           |
+$A_5$                   | [sec]    | `A5`     | Notch numerator coefficient          | 1.013         |
+$A_6$                   | [sec²]   | `A6`     | Notch numerator coefficient          | 0.113         |
+$T_1$                   | [sec]    | `T1`     | Lead-lag 1 numerator time constant   | 0.0           |
+$T_2$                   | [sec]    | `T2`     | Lead-lag 1 denominator time constant | 0.02          |
+$T_3$                   | [sec]    | `T3`     | Lead-lag 2 numerator time constant   | 0.0           |
+$T_4$                   | [sec]    | `T4`     | Lead-lag 2 denominator time constant | 0.0           |
+$T_5$                   | [sec]    | `T5`     | Washout numerator time constant      | 1.65          |
+$T_6$                   | [sec]    | `T6`     | Washout denominator time constant    | 1.65          |
+$K_s$                   | [p.u.]   | `Ks`     | Stabilizer gain                      | 3.0           |
+$L_s^{\min}$            | [p.u.]   | `Lsmin`  | Minimum stabilizer output limit      | -0.1          |
+$L_s^{\max}$            | [p.u.]   | `Lsmax`  | Maximum stabilizer output limit      | 0.1           |
+$V_{\mathrm{cl}}$       | [p.u.]   | `Vcl`    | Lower input cutout threshold         | 0.0           | Nonzero values warn and are ignored
+$V_{\mathrm{cu}}$       | [p.u.]   | `Vcu`    | Upper input cutout threshold         | 0.0           | Nonzero values warn and are ignored
+$T_{\mathrm{delay}}$    | [sec]    | `Tdelay` | Input delay                          | 0.0           | Nonzero values warn and are ignored
 
 ### Parameter Validation
 
-A valid IEEEST parameter set must satisfy the following conditions:
+Let $m$ be the numerator order. Exact zero selects the notch-filter topology:
 
 ```math
-a_2=a_3=a_4=0 \quad\Longrightarrow\quad a_1=0
+\begin{aligned}
+  m &=
+  \begin{cases}
+    2 & A_6 \ne 0 \\
+    1 & A_6 = 0,\ A_5 \ne 0 \\
+    0 & A_5 = A_6 = 0
+  \end{cases} \\
+  m &\le n \\
+  T_2,T_4,T_6 &\ge 0 \\
+  L_s^{\min} &< L_s^{\max} \\
+  a_n &\ne 0 \quad n>0
+\end{aligned}
 ```
 
+The raw and derived coefficients must be finite. The final condition applies to
+the leading coefficient selected by the derived order. Negative denominator
+time constants are rejected.
+
 ### Model Derived Parameters
 
+For one denominator factor, define its exact order by
+
+```math
+d(p,q)=
+\begin{cases}
+  2 & q \ne 0 \\
+  1 & q = 0,\ p \ne 0 \\
+  0 & p = q = 0.
+\end{cases}
+```
+
+The model derives the notch-filter order and expanded denominator coefficients
+once during construction:
+
 ```math
 \begin{aligned}
-a_0 &= 1 \\
-a_1 &= A_1 + A_3 \\
-a_2 &= A_2 + A_4 + A_1 A_3 \\
-a_3 &= A_1 A_4 + A_2 A_3 \\
-a_4 &= A_2 A_4
+  n   &= d(A_1,A_2) + d(A_3,A_4) \\
+  a_1 &= A_1 + A_3 \\
+  a_2 &= A_2 + A_4 + A_1 A_3 \\
+  a_3 &= A_1 A_4 + A_2 A_3 \\
+  a_4 &= A_2 A_4.
 \end{aligned}
 ```
 
-The notch order is
+Let $\epsilon_T=10^{-3}\ \mathrm{s}$. A denominator time constant below
+$\epsilon_T$ is raised to that floor and logs a warning:
 
 ```math
-n = \begin{cases}
-4 & a_4\ne0 \\
-3 & a_4=0,\ a_3\ne0 \\
-2 & a_4=a_3=0,\ a_2\ne0 \\
-0 & a_4=a_3=a_2=a_1=0
-\end{cases}
+T \leftarrow \max(T,\epsilon_T),
+\qquad T\in\{T_2,T_4,T_6\}.
 ```
 
 ## Model Ports
 
 Name     | Port   | Init  | Description
----------|--------|-------|---------------------------------
-`input`  | Input  | Known | Required stabilizer input signal
-`output` | Output | Known | Limited stabilizer output signal
+---------|--------|-------|------
+`input`  | Input  | Known | Stabilizer input signal
+`output` | Output | Known | Stabilizer output signal
+
+The input must be attached to a linked signal. Assigning the output signal is
+optional.
 
 ## Model Variables
 
@@ -77,23 +117,26 @@ Name     | Port   | Init  | Description
 
 #### Differential
 
-Symbol | Units     | Description        | Note
--------|-----------|--------------------|------------------------
-$x_1$  | [p.u.]    | Notch filter state | Held at zero when $n=0$
-$x_2$  | [p.u./s]  | Notch filter state | Held at zero when $n=0$
-$x_3$  | [p.u./s²] | Notch filter state | Held at zero when $n<3$
-$x_4$  | [p.u./s³] | Notch filter state | Held at zero when $n<4$
-$x_5$  | [p.u.]    | Lead–lag 1 state   | Algebraic when $T_2=0$
-$x_6$  | [p.u.]    | Lead–lag 2 state   | Algebraic when $T_4=0$
-$x_7$  | [p.u.]    | Washout state      | Algebraic when $T_6=0$
+All seven differential variables are always present. Inactive notch-filter
+states remain differential and are frozen by $0=-\dot{x}_i$.
+
+Symbol | Units       | Description                           | Note
+-------|-------------|---------------------------------------|------
+$x_1$  | [p.u.]      | Notch-filter signal state             | Active for $n\ge1$; frozen for $n=0$
+$x_2$  | [p.u./sec]  | First derivative of filtered signal   | Active for $n\ge2$; otherwise frozen
+$x_3$  | [p.u./sec²] | Second derivative of filtered signal  | Active for $n\ge3$; otherwise frozen
+$x_4$  | [p.u./sec³] | Third derivative of filtered signal   | Active for $n=4$; otherwise frozen
+$x_5$  | [p.u.]      | Lead-lag 1 state                      | State 5 in Fig. 1
+$x_6$  | [p.u.]      | Lead-lag 2 state                      | State 6 in Fig. 1
+$x_7$  | [p.u.]      | Washout state                         | State 7 in Fig. 1
 
 #### Algebraic
 
 Symbol            | Units  | Description                              | 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_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) |
 
@@ -105,73 +148,104 @@ None.
 
 #### Algebraic
 
-Symbol | Units  | Description             | Note
--------|--------|-------------------------|-----
-$u$    | [p.u.] | Stabilizer input signal |
+Symbol | Units  | Type  | Description             | Note
+-------|--------|-------|-------------------------|------
+$u$    | [p.u.] | Known | Stabilizer input signal | Required signal port `input`
 
 ## Model Equations
 
-Smooth functions: [`clamp`](../../../../CommonMath.md#clamp).
-
-### Internal Equations
+### Differential Equations
 
-#### Differential
+A runtime switch selects one notch-filter realization using the derived order
+$n$. The state layout does not change with $n$.
 
 ```math
 \begin{aligned}
-0 &= -\dot{x}_k+x_{k+1}, && 1\le k<n \\
-0 &= -a_n\dot{x}_n-\sum_{j=0}^{n-1}a_jx_{j+1}+u, && n>0 \\
-0 &= -\dot{x}_k, && n<k\le4 \\
-0 &= -T_2\dot{x}_5-x_5+v_4 \\
-0 &= -T_4\dot{x}_6-x_6+v_5 \\
-0 &= -T_6\dot{x}_7-x_7+v_6
+  0 &= -\dot{x}_1
+    + \begin{cases}
+      0 & n=0 \\
+      \dfrac{-x_1+u}{a_1} & n=1 \\
+      x_2 & n\in\{2,3,4\}
+    \end{cases} \\
+  0 &= -\dot{x}_2
+    + \begin{cases}
+      0 & n\in\{0,1\} \\
+      \dfrac{-x_1-a_1x_2+u}{a_2} & n=2 \\
+      x_3 & n\in\{3,4\}
+    \end{cases} \\
+  0 &= -\dot{x}_3
+    + \begin{cases}
+      0 & n\in\{0,1,2\} \\
+      \dfrac{-x_1-a_1x_2-a_2x_3+u}{a_3} & n=3 \\
+      x_4 & n=4
+    \end{cases} \\
+  0 &= -\dot{x}_4
+    + \begin{cases}
+      0 & n\in\{0,1,2,3\} \\
+      \dfrac{-x_1-a_1x_2-a_2x_3-a_3x_4+u}{a_4} & n=4
+    \end{cases} \\
+  0 &= -\dot{x}_5 + \dfrac{v_4-x_5}{T_2} \\
+  0 &= -\dot{x}_6 + \dfrac{v_5-x_6}{T_4} \\
+  0 &= -\dot{x}_7 + \dfrac{v_6-x_7}{T_6}.
 \end{aligned}
 ```
 
-#### Algebraic
+### Algebraic Equations
 
 ```math
 \begin{aligned}
-0 &= -v_4+\begin{cases}
-x_1+A_5x_2+A_6\dot{x}_2 & n>0 \\
-u & n=0
-\end{cases} \\
-0 &= \begin{cases}
--T_2(v_5-x_5)+T_1(v_4-x_5) & T_2\ne0 \\
--v_5+v_4 & T_2=0
-\end{cases} \\
-0 &= \begin{cases}
--T_4(v_6-x_6)+T_3(v_5-x_6) & T_4\ne0 \\
--v_6+v_5 & T_4=0
-\end{cases} \\
-0 &= \begin{cases}
--T_6v_7+K_\mathrm{s}T_5(v_6-x_7) & T_6\ne0 \\
--v_7+K_\mathrm{s}v_6 & T_6=0
-\end{cases} \\
-0 &= -V_{\mathrm{ss}}+\text{clamp}(v_7;L_\mathrm{s}^{\min},L_\mathrm{s}^{\max})
+  0 &= -v_4
+    + \begin{cases}
+      u & n=0 \\
+      x_1 + \dfrac{A_5}{a_1}(-x_1+u) & n=1 \\
+      x_1 + A_5x_2
+        + \dfrac{A_6}{a_2}(-x_1-a_1x_2+u) & n=2 \\
+      x_1 + A_5x_2 + A_6x_3 & n\in\{3,4\}
+    \end{cases} \\
+  0 &= -v_5 + x_5 + \dfrac{T_1}{T_2}(v_4-x_5) \\
+  0 &= -v_6 + x_6 + \dfrac{T_3}{T_4}(v_5-x_6) \\
+  0 &= -v_7 + K_s\dfrac{T_5}{T_6}(v_6-x_7) \\
+  0 &= -V_{\mathrm{ss}}
+    + \text{clamp}(v_7,L_s^{\min},L_s^{\max}).
 \end{aligned}
 ```
 
-### External Equations
-
-None.
+CommonMath defines the smooth
+[clamp](../../../../CommonMath.md#clamp).
 
 ## Initialization
 
-The initial input determines the steady state; all derivatives initialize to zero.
+### Input Initialization
+
+```math
+u \leftarrow \text{stabilizer input signal}.
+```
+
+### Internal Initialization
 
 ```math
 \begin{aligned}
-x_1 &\leftarrow \begin{cases}u & n>0 \\ 0 & n=0\end{cases} \\
-x_2,x_3,x_4 &\leftarrow 0 \\
-x_5,x_6,x_7,v_4,v_5,v_6 &\leftarrow u \\
-v_7 &\leftarrow \begin{cases}0 & T_6\ne0 \\ K_\mathrm{s}u & T_6=0\end{cases} \\
-V_{\mathrm{ss}} &\leftarrow \text{clamp}(v_7;L_\mathrm{s}^{\min},L_\mathrm{s}^{\max})
+  x_1 &\leftarrow
+    \begin{cases}
+      0 & n=0 \\
+      u & n\in\{1,2,3,4\}
+    \end{cases} \\
+  x_2,x_3,x_4 &\leftarrow 0 \\
+  v_4,x_5,v_5,x_6,v_6,x_7 &\leftarrow u \\
+  v_7 &\leftarrow 0 \\
+  V_{\mathrm{ss}} &\leftarrow
+    \text{clamp}(v_7,L_s^{\min},L_s^{\max}) \\
+  \dot{x}_i &\leftarrow 0,
+    \qquad i\in\{1,\ldots,7\}.
 \end{aligned}
 ```
 
-## Monitors
+### Output Initialization
+
+None.
+
+## Monitorable Outputs
 
-Monitor | Units  | Description               | Note
---------|--------|---------------------------|--------------------------------
-`vss`   | [p.u.] | Limited stabilizer signal | $V_{\mathrm{ss}}$; model output
+Output | Units  | Description                                 | Note
+-------|--------|---------------------------------------------|------
+`vss`  | [p.u.] | Limited stabilizer signal $V_{\mathrm{ss}}$ | Model output
diff --git a/docs/Figures/stabilizer_ieeest_diagram.png b/docs/Figures/PhasorDynamics/IEEEST/diagram.png
similarity index 100%
rename from docs/Figures/stabilizer_ieeest_diagram.png
rename to docs/Figures/PhasorDynamics/IEEEST/diagram.png
diff --git a/tests/UnitTests/PhasorDynamics/ComponentConnectionTests.hpp b/tests/UnitTests/PhasorDynamics/ComponentConnectionTests.hpp
index 5cdb8db8b3..32b69ad8d2 100644
--- a/tests/UnitTests/PhasorDynamics/ComponentConnectionTests.hpp
+++ b/tests/UnitTests/PhasorDynamics/ComponentConnectionTests.hpp
@@ -12,10 +12,13 @@
 #include <GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp>
 #include <GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp>
 #include <GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aData.hpp>
+#include <GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Data.hpp>
 #include <GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPti.hpp>
 #include <GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp>
 #include <GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovData.hpp>
 #include <GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp>
+#include <GridKit/Model/PhasorDynamics/SignalSource/ConstantSignalSourceData.hpp>
+#include <GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestData.hpp>
 #include <GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/Genrou.hpp>
 #include <GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/Gensal.hpp>
 #include <GridKit/Model/PhasorDynamics/SynchronousMachine/GenClassical/GenClassical.hpp>
@@ -414,6 +417,68 @@ namespace GridKit
 
         return success.report(__func__);
       }
+
+      /// IEEEST reinitializes before IEEET1 so the exciter does not latch a
+      /// stale stabilizer output.
+      TestOutcome ieeestIeeet1()
+      {
+        using namespace PhasorDynamics;
+        using namespace PhasorDynamics::Exciter;
+        using namespace PhasorDynamics::Stabilizer;
+
+        TestStatus success = true;
+
+        constexpr IdxT               bus_id   = static_cast<IdxT>(1);
+        constexpr IdxT               input_id = static_cast<IdxT>(20);
+        constexpr IdxT               vs_id    = static_cast<IdxT>(21);
+        SystemModelData<RealT, IdxT> data;
+
+        auto& bus    = data.bus.emplace_back();
+        bus.bus_id   = bus_id;
+        bus.bus_type = BusData<RealT, IdxT>::BusType::SLACK;
+        bus.Vr0      = ONE<RealT>;
+        bus.Vi0      = ZERO<RealT>;
+
+        data.signal = {{"IEEEST input", input_id}, {"IEEEST output", vs_id}};
+
+        auto& source                                                 = data.constant_source.emplace_back();
+        source.signal_outputs[ConstantSignalSourceSignalOutputs::sr] = input_id;
+        source.parameters[ConstantSignalSourceParameters::Sr]        = ZERO<RealT>;
+
+        auto& stabilizer                                       = data.stabilizer.emplace_back();
+        stabilizer.signal_inputs[IeeestSignalInputs::input]    = input_id;
+        stabilizer.signal_outputs[IeeestSignalOutputs::output] = vs_id;
+
+        auto& exciter                                 = data.exciter.emplace_back();
+        exciter.buses[Ieeet1Buses::bus]               = bus_id;
+        exciter.signal_inputs[Ieeet1SignalInputs::vs] = vs_id;
+        exciter.parameters[Ieeet1Parameters::Tr]      = static_cast<RealT>(0.02);
+
+        SystemModel<ScalarT, IdxT> system(data);
+        success *= system.allocate() == 0;
+
+        auto* vs = system.getSignalNode(vs_id);
+        if (vs == nullptr || !vs->linked())
+        {
+          success = false;
+          return success.report(__func__);
+        }
+
+        // Reinitialization must not let the exciter latch this stale value.
+        vs->init(static_cast<ScalarT>(0.05));
+        success *= system.initialize() == 0;
+        success *= system.evaluateResidual() == 0;
+        success *= isEqual(vs->read(), ZERO<ScalarT>, kTol);
+
+        const auto& residual      = system.getResidual();
+        const auto* residual_data = residual.getData();
+        for (IdxT row = 0; row < residual.getSize(); ++row)
+        {
+          success *= isEqual(residual_data[row], ZERO<ScalarT>, kTol);
+        }
+
+        return success.report(__func__);
+      }
     };
 
     /// Production SystemModel wiring tests for GASTPTI and synchronous machines.
diff --git a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp
index 155a88a164..dad3c0553e 100644
--- a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp
+++ b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp
@@ -1,16 +1,25 @@
 #pragma once
 
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <cstddef>
 #include <iomanip>
 #include <iostream>
 #include <limits>
+#include <sstream>
+#include <string>
+#include <utility>
+#include <vector>
 
 #include <GridKit/AutomaticDifferentiation/DependencyTracking/Variable.hpp>
 #include <GridKit/Definitions.hpp>
 #include <GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp>
 #include <GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp>
 #include <GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestData.hpp>
-#include <GridKit/Testing/TestHelpers.hpp>
+#include <GridKit/Model/VariableMonitorController.hpp>
 #include <GridKit/Testing/Testing.hpp>
+#include <GridKit/Testing/Tokenizer.hpp>
 #include <GridKit/Utilities/Logger/Logger.hpp>
 #include <GridKit/Utilities/MapFromCsr.hpp>
 
@@ -20,367 +29,917 @@ namespace GridKit
   {
     using Log = ::GridKit::Utilities::Logger;
 
-    template <class ScalarT, typename IdxT>
+    template <typename scalar_type, typename index_type>
     class StabilizerIeeestTests
     {
     public:
-      using RealT = typename PhasorDynamics::Component<ScalarT, IdxT>::RealT;
+      using ScalarT = scalar_type;
+      using IdxT    = index_type;
+      using RealT   = typename PhasorDynamics::Component<ScalarT, IdxT>::RealT;
 
       StabilizerIeeestTests()  = default;
       ~StabilizerIeeestTests() = default;
 
-      TestOutcome constructor()
+      TestOutcome validation()
       {
         TestStatus success = true;
 
-        auto  data = makeTestData();
-        auto* stab = new PhasorDynamics::Stabilizer::Ieeest<ScalarT, IdxT>(data);
+        IeeestT empty;
+        success *= (empty.size() == static_cast<IdxT>(Utilities::enum_size<Internal>()));
+        success *= (empty.getMonitor() == nullptr);
+        success *= (Utilities::enum_size<Internal>() == 12);
+        success *= (static_cast<size_t>(Internal::VSS) == 11);
 
-        success *= (stab != nullptr);
-        success *= (stab->getMonitor() != nullptr);
+        Fixture<ScalarT> configured(makeOrderData(4));
+        success *= configured.prepare();
+        success *= (configured.model.getMonitor() != nullptr);
 
-        delete stab;
+        for (size_t row = 0; row < configured.model.tag().size(); ++row)
+        {
+          success *= !configured.model.tag()[row];
+        }
+        success *= (configured.model.tagDifferentiable() == 0);
+        for (size_t row = 0; row < configured.model.tag().size(); ++row)
+        {
+          const bool expected = row <= static_cast<size_t>(Internal::X7);
+          if (configured.model.tag()[row] != expected)
+          {
+            std::cout << "IEEEST differentiability tag " << row << " mismatch\n";
+            success = false;
+          }
+        }
+
+        noteExpectedLogs("Testing IEEEST invalid configurations. Logged errors "
+                         "and time-constant or unsupported-feature warnings are expected.");
+
+        for (size_t order = 0; order <= 4; ++order)
+        {
+          success *= verifies(makeOrderData(order));
+        }
+        success *= verifies(makeOrderThreeReverseData());
+        success *= verifies(makeSymmetricOrderFourData());
+
+        // The required input must be both attached and linked.
+        {
+          IeeestT missing(makeOrderData(4));
+          success *= (missing.allocate() == 0);
+          success *= (missing.verify() > 0);
+        }
+        {
+          PhasorDynamics::SignalNode<ScalarT, IdxT> unlinked_node;
+          IeeestT                                   unlinked(makeOrderData(4));
+          unlinked.getPorts().in.template port<PhasorDynamics::Stabilizer::IeeestSignalInputs::input>().connect(&unlinked_node);
+          success *= (unlinked.allocate() == 0);
+          success *= (unlinked.verify() > 0);
+        }
+
+        // The VSS output assignment remains optional.
+        Fixture<ScalarT> no_output(makeOrderData(4), false);
+        success *= no_output.prepare();
+
+        // Integer JSON values are accepted for real parameters; booleans are not.
+        {
+          auto data                    = makeOrderData(4);
+          data.parameters[Params::Ks]  = static_cast<IdxT>(3);
+          success                     *= verifies(data);
+          data.parameters[Params::Ks]  = true;
+          success                     *= !verifies(data);
+        }
+
+        const RealT                  nan      = std::numeric_limits<RealT>::quiet_NaN();
+        const RealT                  infinity = std::numeric_limits<RealT>::infinity();
+        const std::array<RealT, 3>   nonfinite_values{{nan, infinity, -infinity}};
+        const std::array<Params, 18> real_parameters{{
+            Params::A1,
+            Params::A2,
+            Params::A3,
+            Params::A4,
+            Params::A5,
+            Params::A6,
+            Params::T1,
+            Params::T2,
+            Params::T3,
+            Params::T4,
+            Params::T5,
+            Params::T6,
+            Params::Ks,
+            Params::Lsmin,
+            Params::Lsmax,
+            Params::Vcl,
+            Params::Vcu,
+            Params::Tdelay,
+        }};
+
+        for (const Params parameter : real_parameters)
+        {
+          for (const RealT value : nonfinite_values)
+          {
+            auto data                   = makeOrderData(4);
+            data.parameters[parameter]  = value;
+            success                    *= !verifies(data);
+          }
+        }
+
+        success *= invalidParameterLeavesDefault(Params::Ks,
+                                                 true,
+                                                 "boolean parameter fallback");
+        success *= invalidParameterLeavesDefault(Params::T1,
+                                                 nan,
+                                                 "nonfinite parameter fallback");
+
+        // Expanded coefficients must remain finite.
+        {
+          auto data                    = makeOrderData(2);
+          data.parameters[Params::A1]  = std::numeric_limits<RealT>::max();
+          data.parameters[Params::A3]  = std::numeric_limits<RealT>::max();
+          success                     *= !verifies(data);
+        }
+
+        // Finite quadratic factors select order four even when their leading
+        // product overflows.
+        {
+          auto data                    = makeOrderData(0);
+          data.parameters[Params::A2]  = std::numeric_limits<RealT>::max();
+          data.parameters[Params::A4]  = static_cast<RealT>(2.0);
+          success                     *= !verifies(data);
+        }
+
+        // Exact nonzero factors still select order four; an underflowed a4 is invalid.
+        {
+          auto data                    = makeOrderData(0);
+          data.parameters[Params::A2]  = std::numeric_limits<RealT>::min();
+          data.parameters[Params::A4]  = std::numeric_limits<RealT>::min();
+          success                     *= !verifies(data);
+        }
+
+        // The numerator may not have a higher order than the denominator.
+        {
+          auto data                    = makeOrderData(0);
+          data.parameters[Params::A5]  = static_cast<RealT>(0.1);
+          success                     *= !verifies(data);
+          data                         = makeOrderData(0);
+          data.parameters[Params::A6]  = static_cast<RealT>(0.1);
+          success                     *= !verifies(data);
+          data                         = makeOrderData(1);
+          data.parameters[Params::A6]  = static_cast<RealT>(0.1);
+          success                     *= !verifies(data);
+          success                     *= verifies(makeOrderData(2));
+        }
+
+        for (const Params parameter : {Params::T2, Params::T4, Params::T6})
+        {
+          auto data                   = makeOrderData(4);
+          data.parameters[parameter]  = static_cast<RealT>(-0.1);
+          success                    *= !verifies(data);
+        }
+
+        {
+          auto data                       = makeOrderData(4);
+          data.parameters[Params::Lsmin]  = static_cast<RealT>(0.1);
+          data.parameters[Params::Lsmax]  = static_cast<RealT>(0.1);
+          success                        *= !verifies(data);
+          data.parameters[Params::Lsmin]  = static_cast<RealT>(0.2);
+          success                        *= !verifies(data);
+        }
+
+        // Finite unsupported inputs warn but do not invalidate the model.
+        {
+          auto data                        = makeOrderData(4);
+          data.parameters[Params::Vcl]     = static_cast<RealT>(-0.5);
+          data.parameters[Params::Vcu]     = static_cast<RealT>(0.5);
+          data.parameters[Params::Tdelay]  = static_cast<RealT>(0.1);
+          success                         *= verifies(data);
+        }
+
+        // Topology uses exact coefficient zeros, including signed zero.
+        success *= checkTopology(makeOrderData(0), 0, "exact zero");
+        {
+          auto data                    = makeOrderData(0);
+          data.parameters[Params::A1]  = -ZERO<RealT>;
+          data.parameters[Params::A2]  = -ZERO<RealT>;
+          data.parameters[Params::A3]  = -ZERO<RealT>;
+          data.parameters[Params::A4]  = -ZERO<RealT>;
+          success                     *= checkTopology(data, 0, "signed zero");
+        }
+        {
+          auto data                    = makeOrderData(1);
+          data.parameters[Params::A1]  = static_cast<RealT>(-0.4);
+          success                     *= checkTopology(data, 1, "negative coefficient");
+        }
+        {
+          auto data                    = makeOrderData(1);
+          data.parameters[Params::A1]  = static_cast<RealT>(1.0e-12);
+          success                     *= checkTopology(data, 1, "small nonzero coefficient");
+        }
 
         return success.report(__func__);
       }
 
-      /**
-       * @brief All states initialize to zero (stabilizer at rest).
-       * With u = 0, all residuals should be zero.
-       */
-      TestOutcome zeroInitialResidual()
+      TestOutcome initializationAndSignals()
       {
         TestStatus success = true;
 
-        // Create signal nodes for input (u) and output (Vss)
-        PhasorDynamics::SignalNode<ScalarT, IdxT> u_node;
-        PhasorDynamics::SignalNode<ScalarT, IdxT> vss_node;
-        ScalarT                                   u_value{0.0};
-        IdxT                                      u_index = 12; // beyond internal variables
-        ScalarT                                   vss_value{0.0};
-        IdxT                                      vss_index = INVALID_INDEX<IdxT>;
+        for (size_t order = 0; order <= 4; ++order)
+        {
+          Fixture<ScalarT> fixture(makeOrderData(order));
+          fixture.u()  = static_cast<ScalarT>(0.25);
+          success     *= fixture.initialize();
+          success     *= fixture.outputLinked();
+          success     *= (fixture.outputIndex() == static_cast<IdxT>(Internal::VSS));
+
+          const auto* y  = fixture.model.y().getData();
+          const auto* yp = fixture.model.yp().getData();
+          for (size_t row = 0; row < Utilities::enum_size<Internal>(); ++row)
+          {
+            RealT expected = 0.0;
+            if ((order > 0 && row == static_cast<size_t>(Internal::X1))
+                || (row >= static_cast<size_t>(Internal::X5)
+                    && row <= static_cast<size_t>(Internal::V6)))
+            {
+              expected = 0.25;
+            }
+            success *= rowMatches(static_cast<RealT>(y[row]), expected, "state", row, "initialization");
+            success *= rowMatches(static_cast<RealT>(yp[row]), 0.0, "derivative", row, "initialization");
+          }
 
-        // Link signal nodes to backing storage
-        u_node.link(&u_value, &u_index);
-        vss_node.link(&vss_value, &vss_index);
+          success *= (fixture.evaluate() == 0);
+          success *= allResidualsZero(fixture.model);
+          success *= scalarMatches(static_cast<RealT>(fixture.output()),
+                                   static_cast<RealT>(y[static_cast<size_t>(Internal::VSS)]),
+                                   "assigned VSS output");
+        }
 
-        using namespace GridKit::PhasorDynamics::Stabilizer;
+        // Zero denominator time constants use the documented floor and retain
+        // a consistent initial condition.
+        {
+          auto data                   = makeOrderData(4);
+          data.parameters[Params::T2] = ZERO<RealT>;
+          data.parameters[Params::T4] = ZERO<RealT>;
+          data.parameters[Params::T6] = ZERO<RealT>;
+          Fixture<ScalarT> floors(data);
+          floors.u()  = static_cast<ScalarT>(0.25);
+          success    *= floors.initialize();
+          success    *= (floors.evaluate() == 0);
+          success    *= allResidualsZero(floors.model);
+        }
+        {
+          auto data                   = makeOrderData(4);
+          data.parameters[Params::T2] = static_cast<RealT>(5.0e-4);
+          data.parameters[Params::T4] = static_cast<RealT>(5.0e-4);
+          data.parameters[Params::T6] = static_cast<RealT>(5.0e-4);
+          Fixture<ScalarT> floors(data);
+          floors.u()  = static_cast<ScalarT>(0.25);
+          success    *= floors.initialize();
+          success    *= (floors.evaluate() == 0);
+          success    *= allResidualsZero(floors.model);
+        }
 
-        auto data = makeTestData();
-        auto stab = Ieeest<ScalarT, IdxT>(data);
+        // Initialization applies the output limit to the zero unlimited signal.
+        {
+          auto data                      = makeOrderData(4);
+          data.parameters[Params::Lsmin] = static_cast<RealT>(0.2);
+          data.parameters[Params::Lsmax] = static_cast<RealT>(0.6);
+          Fixture<ScalarT> limited(data);
+          limited.u()  = static_cast<ScalarT>(0.25);
+          success     *= limited.initialize();
+          success     *= scalarMatches(static_cast<RealT>(limited.output()),
+                                   static_cast<RealT>(0.2),
+                                   "limited initial VSS",
+                                   kClampTol);
+        }
 
-        // Wire: stabilizer reads u_node as input, writes vss_node as output
-        stab.getPorts().in.template port<IeeestSignalInputs::input>().connect(&u_node);
-        stab.getPorts().out.template port<IeeestSignalOutputs::output>().connect(&vss_node);
+        noteExpectedLogs(
+            "Testing IEEEST rejected initialization paths. "
+            "The logged errors below are expected.");
 
-        stab.allocate();
-        success *= (stab.verify() == 0);
-        stab.initialize();
-        stab.evaluateResidual();
+        // A nonfinite input rejects initialization without changing state.
+        for (const RealT input : {std::numeric_limits<RealT>::quiet_NaN(),
+                                  std::numeric_limits<RealT>::infinity(),
+                                  -std::numeric_limits<RealT>::infinity()})
+        {
+          Fixture<ScalarT> fixture(makeOrderData(4));
+          success *= fixture.prepare();
+          poisonState(fixture.model);
+          const auto y_before   = copyVector(fixture.model.y());
+          const auto yp_before  = copyVector(fixture.model.yp());
+          fixture.u()           = static_cast<ScalarT>(input);
+          success              *= (fixture.model.initialize() != 0);
+          success              *= vectorMatches(fixture.model.y(), y_before, "rejected state");
+          success              *= vectorMatches(fixture.model.yp(), yp_before, "rejected derivative");
+        }
 
-        auto        tol    = 10 * std::numeric_limits<RealT>::epsilon();
-        const auto& f      = stab.getResidual();
-        const auto* f_data = f.getData();
-        for (size_t i = 0; i < f.getSize(); ++i)
+        // Verification failure is also detected before any state mutation.
         {
-          if (!isEqual(f_data[i], 0.0, tol))
-          {
-            std::cout << "Non-zero residual at index " << i << ": " << f_data[i] << "\n";
-            success = false;
-          }
+          IeeestT model(makeOrderData(4));
+          success *= (model.allocate() == 0);
+          poisonState(model);
+          const auto y_before   = copyVector(model.y());
+          const auto yp_before  = copyVector(model.yp());
+          success              *= (model.initialize() != 0);
+          success              *= vectorMatches(model.y(), y_before, "unverified state");
+          success              *= vectorMatches(model.yp(), yp_before, "unverified derivative");
         }
 
-        // Verify output signal is linked and reads the correct value
-        success *= vss_node.linked();
-        success *= (vss_node.getVariableIndex() == 11);
-        success *= isEqual(vss_node.read(), static_cast<ScalarT>(0.0), tol);
+        return success.report(__func__);
+      }
+
+      TestOutcome residualEquations()
+      {
+        TestStatus success = true;
+
+        const std::array<RealT, 12> order_zero{{-0.01, -0.02, -0.03, -0.04, 0.25, 0.24, -0.01, -0.30, -0.25, -0.31, 1.15, 0.0}};
+        const std::array<RealT, 12> order_one{{0.99, -0.02, -0.03, -0.04, 0.25, 0.24, -0.01, -0.20, -0.25, -0.31, 1.15, 0.0}};
+        const std::array<RealT, 12> order_two{{0.19, 1.88, -0.03, -0.04, 0.25, 0.24, -0.01, 0.54, -0.25, -0.31, 1.15, 0.0}};
+        const std::array<RealT, 12> order_three{{0.19, 0.28, 4.153333333333333, -0.04, 0.25, 0.24, -0.01, -0.42, -0.25, -0.31, 1.15, 0.0}};
+        const std::array<RealT, 12> order_three_reverse{{0.19, 0.28, 4.745, -0.04, 0.25, 0.24, -0.01, -0.42, -0.25, -0.31, 1.15, 0.0}};
+        const std::array<RealT, 12> order_four{{0.19, 0.28, 0.37, 1.0975, 0.25, 0.24, -0.01, -0.42, -0.25, -0.31, 1.15, 0.0}};
+        const std::array<RealT, 12> symmetric_four{{0.19, 0.28, 0.37, 2.71, 0.25, 0.24, -0.01, -0.42, -0.25, -0.31, 1.15, 0.0}};
+
+        success *= checkResidual(makeOrderData(0), order_zero, "order 0");
+        success *= checkResidual(makeOrderData(1), order_one, "order 1");
+        success *= checkResidual(makeOrderData(2), order_two, "order 2");
+        success *= checkResidual(makeOrderData(3), order_three, "order 3 first quadratic");
+        success *= checkResidual(makeOrderThreeReverseData(), order_three_reverse, "order 3 second quadratic");
+        success *= checkResidual(makeOrderData(4), order_four, "order 4");
+        success *= checkResidual(makeSymmetricOrderFourData(), symmetric_four, "symmetric order 4");
 
         return success.report(__func__);
       }
 
-      /**
-       * @brief Residual evaluation against hand-computed answer key.
-       *
-       * Sets specific y/yp values and verifies residuals match
-       * pre-computed expected values. See plan for derivation.
-       */
-      TestOutcome residual()
+      TestOutcome monitor()
       {
         TestStatus success = true;
 
-        PhasorDynamics::SignalNode<ScalarT, IdxT> u_node;
-        PhasorDynamics::SignalNode<ScalarT, IdxT> vss_node;
-        ScalarT                                   u_value{0.5};
-        IdxT                                      u_index = 12;
-        ScalarT                                   vss_value{0.0};
-        IdxT                                      vss_index = INVALID_INDEX<IdxT>;
-
-        u_node.link(&u_value, &u_index);
-        vss_node.link(&vss_value, &vss_index);
-
-        using namespace GridKit::PhasorDynamics::Stabilizer;
-
-        auto data = makeTestData();
-        auto stab = Ieeest<ScalarT, IdxT>(data);
-
-        stab.getPorts().in.template port<IeeestSignalInputs::input>().connect(&u_node);
-        stab.getPorts().out.template port<IeeestSignalOutputs::output>().connect(&vss_node);
-
-        stab.allocate();
-        stab.initialize();
-        setStatePoint(stab);
-        stab.evaluateResidual();
-
-        // Hand-computed answer key (see plan for full derivation)
-        const std::vector<ScalarT> 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<RealT>(1.0e-4);
-        auto&       residual      = stab.getResidual();
-        const auto* residual_data = residual.getData();
-
-        for (size_t i = 0; i < res_answer.size(); ++i)
-        {
-          auto test_tol = (i == 11) ? loose_tol : static_cast<RealT>(10 * std::numeric_limits<ScalarT>::epsilon());
-          if (!isEqual(residual_data[i], res_answer[i], test_tol))
-          {
-            std::cout << "Incorrect result for residual " << i << ": "
-                      << std::setprecision(15) << residual_data[i]
-                      << " != " << res_answer[i] << "\n";
-            success = false;
-          }
+        Fixture<ScalarT> fixture(makeOrderData(4));
+        success *= fixture.prepare();
+
+        auto* y                               = fixture.model.y().getData();
+        y[static_cast<size_t>(Internal::VSS)] = static_cast<ScalarT>(0.075);
+        fixture.model.y().setDataUpdated();
+
+        RealT                                     time = 0.0;
+        Model::VariableMonitorController<ScalarT> controller(time);
+        controller.addMonitor(fixture.model.getMonitor());
+
+        std::stringstream os;
+        controller.addSink({Model::VariableMonitorFormat::CSV}, os);
+        controller.start();
+        controller.print();
+        controller.stop();
+
+        std::string header;
+        std::string values;
+        std::getline(os, header);
+        std::getline(os, values);
+        success              *= (header == "t,Ieeest_ieeest_test_vss");
+        const auto monitored  = Tokenizer<RealT>(values, ',')();
+        success              *= (monitored.size() == 2);
+        if (monitored.size() == 2)
+        {
+          success *= scalarMatches(monitored[1], static_cast<RealT>(0.075), "monitored VSS");
         }
 
-        // Verify output signal reads the stabilizer output
-        success *= isEqual(vss_node.read(), static_cast<ScalarT>(0.05), loose_tol);
-
         return success.report(__func__);
       }
 
 #ifdef GRIDKIT_ENABLE_ENZYME
-      /**
-       * @brief Compare DependencyTracking Jacobian against Enzyme Jacobian.
-       */
       TestOutcome jacobian()
       {
         TestStatus success = true;
 
-        auto data = makeTestData();
+        std::vector<std::pair<Data, size_t>> cases;
+        for (size_t order = 0; order <= 4; ++order)
+        {
+          cases.emplace_back(makeOrderData(order), order);
+        }
+        cases.emplace_back(makeOrderThreeReverseData(), 3);
+        cases.emplace_back(makeSymmetricOrderFourData(), 4);
 
-        std::vector<DependencyTracking::Variable::DependencyMap>
-            dependency_tracking_jacobian = DependencyTrackingJacobian(data);
+        for (const auto& [data, order] : cases)
+        {
+          const auto dependency_jacobian = dependencyTrackingJacobian(data, success);
+          const auto enzyme_jacobian     = enzymeJacobian(data, success);
 
-        std::vector<DependencyTracking::Variable::DependencyMap>
-            enzyme_jacobian = EnzymeJacobian(data);
+          success         *= (dependency_jacobian.size() == enzyme_jacobian.size());
+          const auto rows  = std::min(dependency_jacobian.size(), enzyme_jacobian.size());
+          for (size_t row = 0; row < rows; ++row)
+          {
+            const auto& dependency_map = dependency_jacobian[row].getDependencies();
+            if (!isEqual(dependency_map, enzyme_jacobian[row], kTol))
+            {
+              std::cout << "IEEEST Jacobian row " << row << " for order " << order
+                        << " differs between dependency tracking and Enzyme\n"
+                        << "  dependency tracking: "
+                        << std::setprecision(std::numeric_limits<RealT>::max_digits10);
+              dependency_jacobian[row].print(std::cout);
+
+              std::cout << "\n  Enzyme:";
+              for (const auto& [column, value] : enzyme_jacobian[row])
+              {
+                std::cout << " (" << column << ", " << value << ')';
+              }
+              std::cout << '\n';
+              success = false;
+            }
+          }
 
-        // Compare DependencyTracking dependencies to Enzyme's
-        auto tol = 10 * std::numeric_limits<RealT>::epsilon();
-        for (size_t i = 0; i < dependency_tracking_jacobian.size(); ++i)
-        {
-          success *= (GridKit::Testing::isEqual(dependency_tracking_jacobian[i], enzyme_jacobian[i], tol));
+          for (size_t row = order; row < 4; ++row)
+          {
+            const DependencyMap expected{{row, static_cast<RealT>(-1.0)}};
+            const auto&         dependency_map = dependency_jacobian[row].getDependencies();
+            if (!isEqual(dependency_map, expected, kTol)
+                || !isEqual(enzyme_jacobian[row], expected, kTol))
+            {
+              std::cout << "IEEEST inactive notch row " << row << " for order "
+                        << order << " is not the frozen derivative diagonal\n";
+              success = false;
+            }
+          }
         }
 
         return success.report(__func__);
       }
+#endif
 
     private:
-      std::vector<DependencyTracking::Variable::DependencyMap> DependencyTrackingJacobian(
-          PhasorDynamics::Stabilizer::IeeestData<RealT, IdxT> ieeestdata)
+      using Params        = PhasorDynamics::Stabilizer::IeeestParameters;
+      using Internal      = PhasorDynamics::Stabilizer::IeeestInternalVariables;
+      using External      = PhasorDynamics::Stabilizer::IeeestExternalVariables;
+      using Mon           = PhasorDynamics::Stabilizer::IeeestMonitorableVariables;
+      using Data          = PhasorDynamics::Stabilizer::IeeestData<RealT, IdxT>;
+      using IeeestT       = PhasorDynamics::Stabilizer::Ieeest<ScalarT, IdxT>;
+      using DependencyMap = DependencyTracking::Variable::DependencyMap;
+
+      static constexpr RealT kTol =
+          static_cast<RealT>(100.0) * std::numeric_limits<RealT>::epsilon();
+      static constexpr RealT kClampTol = static_cast<RealT>(1.0e-4);
+
+      static constexpr std::array<const char*, 12> kRowNames{{"X1", "X2", "X3", "X4", "X5", "X6", "X7", "V4", "V5", "V6", "V7", "VSS"}};
+
+      template <typename T>
+      class Fixture
       {
-        using DepVar = DependencyTracking::Variable;
+      private:
+        T                                   u_value_{0};
+        IdxT                                u_index_{static_cast<IdxT>(Utilities::enum_size<Internal>())};
+        PhasorDynamics::SignalNode<T, IdxT> u_node_;
+        PhasorDynamics::SignalNode<T, IdxT> vss_node_;
+
+      public:
+        explicit Fixture(const Data& data, bool assign_output = true)
+          : model(data)
+        {
+          u_node_.link(&u_value_, &u_index_);
+          model.getPorts().in.template port<PhasorDynamics::Stabilizer::IeeestSignalInputs::input>().connect(&u_node_);
+          if (assign_output)
+          {
+            model.getPorts().out.template port<PhasorDynamics::Stabilizer::IeeestSignalOutputs::output>().connect(&vss_node_);
+          }
+        }
 
-        // Set up signal nodes with DependencyTracking scalar type
-        PhasorDynamics::SignalNode<DepVar, IdxT> u_node;
-        PhasorDynamics::SignalNode<DepVar, IdxT> vss_node;
-        DepVar                                   u_value{0.5};
-        IdxT                                     u_index = 12;
-        DepVar                                   vss_value{0.0};
-        IdxT                                     vss_index = INVALID_INDEX<IdxT>;
+        Fixture(const Fixture&)            = delete;
+        Fixture& operator=(const Fixture&) = delete;
 
-        u_node.link(&u_value, &u_index);
-        vss_node.link(&vss_value, &vss_index);
+        bool prepare()
+        {
+          return model.allocate() == 0 && model.verify() == 0;
+        }
 
-        using namespace GridKit::PhasorDynamics::Stabilizer;
+        bool initialize()
+        {
+          return prepare() && model.initialize() == 0;
+        }
 
-        auto stab = Ieeest<DepVar, IdxT>(ieeestdata);
-        stab.getPorts().in.template port<IeeestSignalInputs::input>().connect(&u_node);
-        stab.getPorts().out.template port<IeeestSignalOutputs::output>().connect(&vss_node);
+        int evaluate()
+        {
+          return model.evaluateResidual();
+        }
 
-        stab.allocate();
-        stab.initialize();
+        T& u()
+        {
+          return u_value_;
+        }
 
-        // Tag external signal u as an additional independent variable
-        u_value.setVariableNumber(2 * stab.size()); // @todo avoid requiring knowledge of the numbering
-        setStatePointDep(stab);
+        IdxT uIndex() const
+        {
+          return u_index_;
+        }
 
-        stab.updateTime(0.0, 1.0); // alpha = 1.0 to verify d/dy' term
+        T output() const
+        {
+          return vss_node_.read();
+        }
 
-        stab.evaluateResidual();
-        stab.evaluateJacobian();
-        auto  model_jacobian = stab.getCsrJacobian();
-        auto& output_stream  = Log::misc();
-        output_stream << "Sparse Csr Matrix: Ieeest DependencyTracking Jacobian\n";
-        model_jacobian->print(output_stream);
+        bool outputLinked() const
+        {
+          return vss_node_.linked();
+        }
 
-        return GridKit::Testing::MapFromCsr(model_jacobian);
+        IdxT outputIndex() const
+        {
+          return vss_node_.getVariableIndex();
+        }
+
+        PhasorDynamics::Stabilizer::Ieeest<T, IdxT> model;
+      };
+
+      Data makeBaseData() const
+      {
+        Data data;
+        data.device_class          = "Ieeest";
+        data.disambiguation_string = "ieeest_test";
+        data.monitored_variables.insert(Mon::vss);
+
+        data.parameters[Params::A1]     = ZERO<RealT>;
+        data.parameters[Params::A2]     = ZERO<RealT>;
+        data.parameters[Params::A3]     = ZERO<RealT>;
+        data.parameters[Params::A4]     = ZERO<RealT>;
+        data.parameters[Params::A5]     = ZERO<RealT>;
+        data.parameters[Params::A6]     = ZERO<RealT>;
+        data.parameters[Params::T1]     = static_cast<RealT>(0.5);
+        data.parameters[Params::T2]     = static_cast<RealT>(1.0);
+        data.parameters[Params::T3]     = static_cast<RealT>(0.3);
+        data.parameters[Params::T4]     = static_cast<RealT>(1.0);
+        data.parameters[Params::T5]     = static_cast<RealT>(2.0);
+        data.parameters[Params::T6]     = static_cast<RealT>(5.0);
+        data.parameters[Params::Ks]     = static_cast<RealT>(10.0);
+        data.parameters[Params::Lsmin]  = static_cast<RealT>(-0.1);
+        data.parameters[Params::Lsmax]  = static_cast<RealT>(0.1);
+        data.parameters[Params::Vcl]    = ZERO<RealT>;
+        data.parameters[Params::Vcu]    = ZERO<RealT>;
+        data.parameters[Params::Tdelay] = ZERO<RealT>;
+        return data;
       }
 
-      std::vector<DependencyTracking::Variable::DependencyMap> EnzymeJacobian(
-          PhasorDynamics::Stabilizer::IeeestData<RealT, IdxT> ieeestdata)
+      Data makeOrderData(size_t order) const
       {
-        PhasorDynamics::SignalNode<ScalarT, IdxT> u_node;
-        PhasorDynamics::SignalNode<ScalarT, IdxT> vss_node;
-        ScalarT                                   u_value{0.5};
-        IdxT                                      u_index = 12;
-        ScalarT                                   vss_value{0.0};
-        IdxT                                      vss_index = INVALID_INDEX<IdxT>;
+        auto data = makeBaseData();
+        switch (order)
+        {
+        case 0:
+          break;
+        case 1:
+          data.parameters[Params::A1] = static_cast<RealT>(0.4);
+          data.parameters[Params::A5] = static_cast<RealT>(0.5);
+          break;
+        case 2:
+          data.parameters[Params::A1] = static_cast<RealT>(0.1);
+          data.parameters[Params::A2] = static_cast<RealT>(0.2);
+          data.parameters[Params::A5] = static_cast<RealT>(0.5);
+          data.parameters[Params::A6] = static_cast<RealT>(0.6);
+          break;
+        case 3:
+          data.parameters[Params::A1] = static_cast<RealT>(0.1);
+          data.parameters[Params::A2] = static_cast<RealT>(0.2);
+          data.parameters[Params::A3] = static_cast<RealT>(0.3);
+          data.parameters[Params::A5] = static_cast<RealT>(0.5);
+          data.parameters[Params::A6] = static_cast<RealT>(0.6);
+          break;
+        case 4:
+          data.parameters[Params::A1] = static_cast<RealT>(0.1);
+          data.parameters[Params::A2] = static_cast<RealT>(0.2);
+          data.parameters[Params::A3] = static_cast<RealT>(0.3);
+          data.parameters[Params::A4] = static_cast<RealT>(0.4);
+          data.parameters[Params::A5] = static_cast<RealT>(0.5);
+          data.parameters[Params::A6] = static_cast<RealT>(0.6);
+          break;
+        default:
+          break;
+        }
+        return data;
+      }
 
-        u_node.link(&u_value, &u_index);
-        vss_node.link(&vss_value, &vss_index);
+      Data makeOrderThreeReverseData() const
+      {
+        auto data                   = makeBaseData();
+        data.parameters[Params::A1] = static_cast<RealT>(0.1);
+        data.parameters[Params::A3] = static_cast<RealT>(0.3);
+        data.parameters[Params::A4] = static_cast<RealT>(0.4);
+        data.parameters[Params::A5] = static_cast<RealT>(0.5);
+        data.parameters[Params::A6] = static_cast<RealT>(0.6);
+        return data;
+      }
 
-        using namespace GridKit::PhasorDynamics::Stabilizer;
+      Data makeSymmetricOrderFourData() const
+      {
+        auto data                   = makeOrderData(4);
+        data.parameters[Params::A1] = ZERO<RealT>;
+        data.parameters[Params::A3] = ZERO<RealT>;
+        return data;
+      }
 
-        auto stab = Ieeest<ScalarT, IdxT>(ieeestdata);
-        stab.getPorts().in.template port<IeeestSignalInputs::input>().connect(&u_node);
-        stab.getPorts().out.template port<IeeestSignalOutputs::output>().connect(&vss_node);
+      bool verifies(const Data& data) const
+      {
+        Fixture<ScalarT> fixture(data, false);
+        return fixture.prepare();
+      }
 
-        stab.allocate();
-        stab.initialize();
-        setStatePoint(stab);
+      template <typename ValueT>
+      bool invalidParameterLeavesDefault(Params      parameter,
+                                         ValueT      invalid_value,
+                                         const char* label) const
+      {
+        auto expected_data = makeOrderData(4);
+        expected_data.parameters.erase(parameter);
+        auto invalid_data                  = expected_data;
+        invalid_data.parameters[parameter] = invalid_value;
+
+        Fixture<ScalarT> expected(expected_data);
+        Fixture<ScalarT> invalid(invalid_data);
+        expected.u() = static_cast<ScalarT>(0.5);
+        invalid.u()  = static_cast<ScalarT>(0.5);
+
+        bool success = expected.prepare();
+        success      = success && invalid.model.allocate() == 0;
+        success      = success && invalid.model.verify() > 0;
+        setAnswerKey(expected.model);
+        setAnswerKey(invalid.model);
+        success = success && expected.evaluate() == 0;
+        success = success && invalid.evaluate() == 0;
+        success = success && vectorMatches(invalid.model.getResidual(), copyVector(expected.model.getResidual()), label);
+        return success;
+      }
 
-        stab.updateTime(0.0, 1.0); // alpha = 1.0 to verify d/dy' term
+      bool checkTopology(const Data& data, size_t order, const char* label) const
+      {
+        Fixture<ScalarT> fixture(data);
+        fixture.u() = static_cast<ScalarT>(0.5);
+        if (!fixture.initialize())
+        {
+          std::cout << "IEEEST " << label << " topology failed to initialize\n";
+          return false;
+        }
+
+        setAnswerKey(fixture.model);
+        if (fixture.evaluate() != 0)
+        {
+          return false;
+        }
 
-        stab.evaluateResidual();
-        stab.evaluateJacobian();
-        stab.constructCsr();
-        auto  model_jacobian = stab.getCsrJacobian();
-        auto& output_stream  = Log::misc();
-        output_stream << "Sparse Csr Matrix: Ieeest Enzyme Jacobian\n";
-        model_jacobian->print(output_stream);
+        bool        success = true;
+        const auto* f       = fixture.model.getResidual().getData();
+        const auto* yp      = fixture.model.yp().getData();
+        for (size_t row = 0; row < 4; ++row)
+        {
+          const bool frozen = isEqual(static_cast<RealT>(f[row]),
+                                      -static_cast<RealT>(yp[row]),
+                                      kTol);
+          if (frozen != (row >= order))
+          {
+            std::cout << "IEEEST " << label << " topology row " << row << " mismatch\n";
+            success = false;
+          }
+        }
+        return success;
+      }
 
-        return GridKit::Testing::MapFromCsr(model_jacobian);
+      static constexpr std::array<RealT, 12> answerState()
+      {
+        return {{0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 0.05, 0.05}};
       }
-#endif
 
-    private:
-      static constexpr ScalarT tol_ = 10 * std::numeric_limits<ScalarT>::epsilon();
+      static constexpr std::array<RealT, 12> answerDerivative()
+      {
+        return {{0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.0, 0.0, 0.0, 0.0, 0.0}};
+      }
 
-      /**
-       * @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<RealT, IdxT>
+      template <typename T>
+      void setAnswerKey(PhasorDynamics::Stabilizer::Ieeest<T, IdxT>& model) const
       {
-        using Params = PhasorDynamics::Stabilizer::IeeestParameters;
+        const auto state      = answerState();
+        const auto derivative = answerDerivative();
+        auto*      y          = model.y().getData();
+        auto*      yp         = model.yp().getData();
+        for (size_t row = 0; row < state.size(); ++row)
+        {
+          y[row]  = static_cast<T>(state[row]);
+          yp[row] = static_cast<T>(derivative[row]);
+        }
+        model.y().setDataUpdated();
+        model.yp().setDataUpdated();
+      }
 
-        PhasorDynamics::Stabilizer::IeeestData<RealT, IdxT> data;
-        data.device_class          = "stabilizer";
-        data.disambiguation_string = "ieeest_test";
-        data.monitored_variables.insert(PhasorDynamics::Stabilizer::IeeestMonitorableVariables::vss);
-
-        data.parameters[Params::A1]     = 0.1;
-        data.parameters[Params::A2]     = 0.2;
-        data.parameters[Params::A3]     = 0.3;
-        data.parameters[Params::A4]     = 0.4;
-        data.parameters[Params::A5]     = 0.5;
-        data.parameters[Params::A6]     = 0.6;
-        data.parameters[Params::T1]     = 0.5;
-        data.parameters[Params::T2]     = 1.0;
-        data.parameters[Params::T3]     = 0.3;
-        data.parameters[Params::T4]     = 1.0;
-        data.parameters[Params::T5]     = 2.0;
-        data.parameters[Params::T6]     = 5.0;
-        data.parameters[Params::Ks]     = 10.0;
-        data.parameters[Params::Lsmin]  = -0.1;
-        data.parameters[Params::Lsmax]  = 0.1;
-        data.parameters[Params::Vcl]    = 0.0;
-        data.parameters[Params::Vcu]    = 0.0;
-        data.parameters[Params::Tdelay] = 0.0;
+      bool checkResidual(const Data&                  data,
+                         const std::array<RealT, 12>& expected,
+                         const char*                  label) const
+      {
+        Fixture<ScalarT> fixture(data);
+        fixture.u() = static_cast<ScalarT>(0.5);
+        if (!fixture.initialize())
+        {
+          return false;
+        }
+        setAnswerKey(fixture.model);
+        if (fixture.evaluate() != 0)
+        {
+          return false;
+        }
 
-        return data;
+        bool        success = true;
+        const auto* f       = fixture.model.getResidual().getData();
+        for (size_t row = 0; row < expected.size(); ++row)
+        {
+          RealT tolerance = kTol;
+          if (row == static_cast<size_t>(Internal::VSS))
+          {
+            tolerance = kClampTol;
+          }
+          if (!rowMatches(static_cast<RealT>(f[row]), expected[row], "residual", row, label, tolerance))
+          {
+            success = false;
+          }
+        }
+        return success;
+      }
+
+      template <typename ModelT>
+      bool allResidualsZero(const ModelT& model) const
+      {
+        bool        success = true;
+        const auto* f       = model.getResidual().getData();
+        for (size_t row = 0; row < static_cast<size_t>(model.getResidual().getSize()); ++row)
+        {
+          if (!rowMatches(static_cast<RealT>(f[row]), 0.0, "residual", row, "at rest", kClampTol))
+          {
+            success = false;
+          }
+        }
+        return success;
       }
 
-      /**
-       * @brief Set a non-trivial operating point for residual/Jacobian tests.
-       * Avoids zeros and ones to catch coefficient errors.
-       */
-      void setStatePoint(PhasorDynamics::Stabilizer::Ieeest<ScalarT, IdxT>& stab)
+      bool rowMatches(RealT       actual,
+                      RealT       expected,
+                      const char* what,
+                      size_t      row,
+                      const char* context,
+                      RealT       tolerance = kTol) const
       {
-        auto* y  = stab.y().getData();
-        auto* yp = stab.yp().getData();
-
-        y[0]  = 0.1;  // x1
-        y[1]  = 0.2;  // x2
-        y[2]  = 0.3;  // x3
-        y[3]  = 0.4;  // x4
-        y[4]  = 0.5;  // x5
-        y[5]  = 0.6;  // x6
-        y[6]  = 0.7;  // x7
-        y[7]  = 0.8;  // v4
-        y[8]  = 0.9;  // v5
-        y[9]  = 1.0;  // v6
-        y[10] = 0.05; // v7  (within limiter range)
-        y[11] = 0.05; // Vss (model output)
-
-        yp[0] = 0.01; // x1_dot
-        yp[1] = 0.02; // x2_dot
-        yp[2] = 0.03; // x3_dot
-        yp[3] = 0.04; // x4_dot
-        yp[4] = 0.05; // x5_dot
-        yp[5] = 0.06; // x6_dot
-        yp[6] = 0.07; // x7_dot
-
-        stab.y().setDataUpdated();
-        stab.yp().setDataUpdated();
+        if (isEqual(actual, expected, tolerance))
+        {
+          return true;
+        }
+        std::cout << "IEEEST " << what << " row " << kRowNames[row] << ' '
+                  << context << " mismatch: "
+                  << std::setprecision(std::numeric_limits<RealT>::max_digits10)
+                  << actual << " != " << expected << '\n';
+        return false;
       }
 
-      /**
-       * @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<DependencyTracking::Variable, IdxT>& stab)
+      bool scalarMatches(RealT       actual,
+                         RealT       expected,
+                         const char* label,
+                         RealT       tolerance = kTol) const
       {
-        auto* y  = stab.y().getData();
-        auto* yp = stab.yp().getData();
-
-        y[0].setValue(0.1);
-        y[1].setValue(0.2);
-        y[2].setValue(0.3);
-        y[3].setValue(0.4);
-        y[4].setValue(0.5);
-        y[5].setValue(0.6);
-        y[6].setValue(0.7);
-        y[7].setValue(0.8);
-        y[8].setValue(0.9);
-        y[9].setValue(1.0);
-        y[10].setValue(0.05);
-        y[11].setValue(0.05);
-
-        yp[0].setValue(0.01);
-        yp[1].setValue(0.02);
-        yp[2].setValue(0.03);
-        yp[3].setValue(0.04);
-        yp[4].setValue(0.05);
-        yp[5].setValue(0.06);
-        yp[6].setValue(0.07);
-
-        stab.y().setDataUpdated();
-        stab.yp().setDataUpdated();
+        if (isEqual(actual, expected, tolerance))
+        {
+          return true;
+        }
+        std::cout << "IEEEST " << label << " mismatch: "
+                  << std::setprecision(std::numeric_limits<RealT>::max_digits10)
+                  << actual << " != " << expected << '\n';
+        return false;
       }
-    }; // class StabilizerIeeestTests
+
+      template <typename VectorT>
+      std::vector<RealT> copyVector(const VectorT& vector) const
+      {
+        const auto*        values = vector.getData();
+        std::vector<RealT> copy(static_cast<size_t>(vector.getSize()));
+        for (size_t row = 0; row < copy.size(); ++row)
+        {
+          copy[row] = static_cast<RealT>(values[row]);
+        }
+        return copy;
+      }
+
+      template <typename VectorT>
+      bool vectorMatches(const VectorT&            vector,
+                         const std::vector<RealT>& expected,
+                         const char*               label) const
+      {
+        bool        success = true;
+        const auto* values  = vector.getData();
+        for (size_t row = 0; row < expected.size(); ++row)
+        {
+          success *= rowMatches(static_cast<RealT>(values[row]),
+                                expected[row],
+                                label,
+                                row,
+                                "atomic rejection");
+        }
+        return success;
+      }
+
+      template <typename ModelT>
+      void poisonState(ModelT& model) const
+      {
+        auto* y  = model.y().getData();
+        auto* yp = model.yp().getData();
+        for (size_t row = 0; row < static_cast<size_t>(model.size()); ++row)
+        {
+          const RealT offset = static_cast<RealT>(row);
+          y[row]             = static_cast<typename ModelT::ScalarT>(0.125 + 0.01 * offset);
+          yp[row]            = static_cast<typename ModelT::ScalarT>(-0.25 - 0.01 * offset);
+        }
+        model.y().setDataUpdated();
+        model.yp().setDataUpdated();
+      }
+
+      void noteExpectedLogs(const char* message) const
+      {
+        const auto previous_verbosity = Log::verbosity();
+        Log::setVerbosity(Log::Verbosity::EVERYTHING);
+        Log::misc() << message << '\n';
+        Log::setVerbosity(previous_verbosity);
+      }
+
+#ifdef GRIDKIT_ENABLE_ENZYME
+      std::vector<DependencyTracking::Variable> dependencyTrackingJacobian(
+          const Data& data,
+          TestStatus& success) const
+      {
+        using DepVar = DependencyTracking::Variable;
+
+        Fixture<DepVar> fixture(data);
+        fixture.u()  = static_cast<DepVar>(0.5);
+        success     *= fixture.initialize();
+        setAnswerKey(fixture.model);
+
+        auto* y  = fixture.model.y().getData();
+        auto* yp = fixture.model.yp().getData();
+        for (size_t row = 0; row < static_cast<size_t>(fixture.model.size()); ++row)
+        {
+          y[row].setVariableNumber(row);
+        }
+        fixture.u().setVariableNumber(fixture.uIndex());
+        fixture.model.y().setDataUpdated();
+        success *= (fixture.evaluate() == 0);
+
+        const auto*         f = fixture.model.getResidual().getData();
+        std::vector<DepVar> rows(f, f + static_cast<size_t>(fixture.model.size()));
+
+        // Dependency tracking cannot safely alias y and yp to the same
+        // variable number in one pass. Evaluate df/dyp separately and add it
+        // to df/dy, matching Enzyme's df/dy + alpha*df/dyp at alpha = 1.
+        setAnswerKey(fixture.model);
+        fixture.u() = 0.5;
+        for (size_t row = 0; row < static_cast<size_t>(fixture.model.size()); ++row)
+        {
+          yp[row].setVariableNumber(row);
+        }
+        fixture.model.yp().setDataUpdated();
+        success *= (fixture.evaluate() == 0);
+
+        f = fixture.model.getResidual().getData();
+        for (size_t row = 0; row < rows.size(); ++row)
+        {
+          // Add df/dyp dependencies without adding a second residual value.
+          DepVar yp_row = f[row];
+          yp_row.setValue(0.0);
+          rows[row] += yp_row;
+        }
+        // Check the new single-pass CSR assembly against Enzyme as well.
+        Fixture<DepVar> automatic(data);
+        automatic.u()  = 0.5;
+        success       *= automatic.initialize();
+        setAnswerKey(automatic.model);
+        auto* automatic_y  = automatic.model.y().getData();
+        auto* automatic_yp = automatic.model.yp().getData();
+        for (IdxT row = 0; row < automatic.model.size(); ++row)
+        {
+          const auto index = static_cast<size_t>(automatic.model.getVariableIndex(row));
+          automatic_y[row].setVariableNumber(2 * index);
+          automatic_yp[row].setVariableNumber(2 * index + 1);
+        }
+        automatic.u().setVariableNumber(2 * static_cast<size_t>(automatic.uIndex()));
+        automatic.model.updateTime(0.0, 1.0);
+        success                 *= (automatic.evaluate() == 0);
+        success                 *= (automatic.model.evaluateJacobian() == 0);
+        const auto assembled     = MapFromCsr(automatic.model.getCsrJacobian());
+        const auto expected_csr  = enzymeJacobian(data, success);
+        success                 *= (assembled.size() == expected_csr.size());
+        for (size_t row = 0; row < std::min(assembled.size(), expected_csr.size()); ++row)
+        {
+          success *= isEqual(assembled[row], expected_csr[row], kTol);
+        }
+        return rows;
+      }
+
+      std::vector<DependencyMap> enzymeJacobian(const Data& data,
+                                                TestStatus& success) const
+      {
+        Fixture<ScalarT> fixture(data);
+        fixture.u()  = static_cast<ScalarT>(0.5);
+        success     *= fixture.initialize();
+        setAnswerKey(fixture.model);
+        fixture.model.updateTime(0.0, 1.0);
+        success *= (fixture.evaluate() == 0);
+        success *= (fixture.model.evaluateJacobian() == 0);
+        success *= (fixture.model.constructCsr() == 0);
+        return MapFromCsr(fixture.model.getCsrJacobian());
+      }
+#endif
+    };
 
   } // namespace Testing
 } // namespace GridKit
diff --git a/tests/UnitTests/PhasorDynamics/runComponentConnectionTests.cpp b/tests/UnitTests/PhasorDynamics/runComponentConnectionTests.cpp
index a84c7cd899..757cc5afdc 100644
--- a/tests/UnitTests/PhasorDynamics/runComponentConnectionTests.cpp
+++ b/tests/UnitTests/PhasorDynamics/runComponentConnectionTests.cpp
@@ -13,6 +13,7 @@ int main()
   result += test.genClassicalEsdc1a();
   result += test.genClassicalHygov();
   result += test.regcaRepca();
+  result += test.ieeestIeeet1();
   result += gastpti.genrouGastPti();
   result += gastpti.gensalGastPti();
   result += gastpti.genClassicalGastPti();
diff --git a/tests/UnitTests/PhasorDynamics/runStabilizerIeeestTests.cpp b/tests/UnitTests/PhasorDynamics/runStabilizerIeeestTests.cpp
index c1999c24aa..2b4022e593 100644
--- a/tests/UnitTests/PhasorDynamics/runStabilizerIeeestTests.cpp
+++ b/tests/UnitTests/PhasorDynamics/runStabilizerIeeestTests.cpp
@@ -6,9 +6,11 @@ int main()
 
   GridKit::Testing::StabilizerIeeestTests<double, size_t> test;
 
-  result += test.constructor();
-  result += test.zeroInitialResidual();
-  result += test.residual();
+  result += test.validation();
+  result += test.initializationAndSignals();
+  result += test.residualEquations();
+  result += test.monitor();
+
 #ifdef GRIDKIT_ENABLE_ENZYME
   result += test.jacobian();
 #endif