Edge case gates for IEEET1 (Validation)

GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp

@@ -123,26 +123,26 @@ namespace GridKit
             const ScalarT*, const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*);
 
       private:
-        static constexpr RealT TR_MINIMUM = 1.0e-3;
+        static constexpr RealT TIME_CONSTANT_MINIMUM = static_cast<RealT>(1.0e-3);
 
         // Signal pointers
         BusT* bus_;
 
         // Model Input parameters
-        RealT Tr_;      ///< Time constant for voltage sensing
-        RealT Ka_;      ///< Coefficient for voltage regulation
-        RealT Ta_;      ///< Time constant for voltage regulation
-        RealT Ke_;      ///< Coefficient for excitation system
-        RealT Te_;      ///< Time constant for excitation system
-        RealT Kf_;      ///< Coefficient for feedback
-        RealT Tf_;      ///< Time constant for feedback
-        RealT Vrmin_;   ///< LL to voltage regulation
-        RealT Vrmax_;   ///< HH to voltage regulation
-        RealT E1_;      ///< Saturation parameter
-        RealT E2_;      ///< Saturation parameter
-        RealT Se1_;     ///< Saturation parameter
-        RealT Se2_;     ///< Saturation parameter
-        RealT Ispdlim_; ///< Speed limit flag indicator
+        RealT Tr_{0.0};      ///< Time constant for voltage sensing
+        RealT Ka_{50.0};     ///< Coefficient for voltage regulation
+        RealT Ta_{0.04};     ///< Time constant for voltage regulation
+        RealT Ke_{-0.06};    ///< Coefficient for excitation system
+        RealT Te_{0.6};      ///< Time constant for excitation system
+        RealT Kf_{0.09};     ///< Coefficient for feedback
+        RealT Tf_{1.46};     ///< Time constant for feedback
+        RealT Vrmin_{-1.0};  ///< LL to voltage regulation
+        RealT Vrmax_{1.0};   ///< HH to voltage regulation
+        RealT E1_{2.8};      ///< Saturation parameter
+        RealT E2_{3.73};     ///< Saturation parameter
+        RealT Se1_{0.04};    ///< Saturation parameter
+        RealT Se2_{0.33};    ///< Saturation parameter
+        RealT Ispdlim_{0.0}; ///< Speed limit flag indicator
 
         // Model Derived parameters
         // Saturation coefficients derived from E1, E2, Se1, and Se2.

GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp

@@ -7,6 +7,7 @@
  *
  */
 
+#include <algorithm>
 #include <cmath>
 #include <iostream>
 
@@ -29,9 +30,8 @@ namespace GridKit
        */
       template <typename scalar_type, typename index_type>
       Ieeet1<scalar_type, index_type>::Ieeet1(BusT* bus)
-        : bus_(bus)
+        : Ieeet1(bus, ModelDataT{})
       {
-        size_ = 9;
       }
 
       /**
@@ -46,7 +46,6 @@ namespace GridKit
         : bus_(bus),
           monitor_(std::make_unique<MonitorT>(data))
       {
-
         // Parse data struct into model
         this->initModelParams(data);
 
@@ -116,7 +115,7 @@ namespace GridKit
       }
 
       /**
-       * @brief verify method checks that attached signals are also linked
+       * @brief Verify parameter values and attached signal links
        */
       template <typename scalar_type, typename index_type>
       int Ieeet1<scalar_type, index_type>::verify() const
@@ -126,6 +125,33 @@ namespace GridKit
 
         int ret = 0;
 
+        auto check = [&](bool condition, const char* message)
+        {
+          if (!condition)
+          {
+            Log::error() << "Ieeet1: " << message << '\n';
+            ret += 1;
+          }
+        };
+
+        check(Ka_ > ZERO<RealT>, "Ka must be positive");
+        check(Vrmin_ <= Vrmax_, "Vrmin must be less than or equal to Vrmax");
+        check(Ispdlim_ == ZERO<RealT> || Ispdlim_ == ONE<RealT>,
+              "Ispdlim must be 0 or 1");
+
+        const bool saturation_disabled =
+            Se1_ == ZERO<RealT> && Se2_ == ZERO<RealT>;
+
+        if (!saturation_disabled)
+        {
+          check(E1_ > ZERO<RealT>, "E1 must be positive when saturation is enabled");
+          check(E2_ > ZERO<RealT>, "E2 must be positive when saturation is enabled");
+          check(Se1_ > ZERO<RealT>, "Se1 must be positive when saturation is enabled");
+          check(Se2_ > ZERO<RealT>, "Se2 must be positive when saturation is enabled");
+          check(E1_ != E2_, "E1 and E2 must differ when saturation is enabled");
+          check(Se1_ != Se2_, "Se1 and Se2 must differ when saturation is enabled");
+        }
+
         if (signals_.template isAttached<OMEGA>())
         {
           if (!signals_.template isLinked<OMEGA>())
@@ -158,7 +184,7 @@ namespace GridKit
        *   - Bus voltage, used to form the sensed terminal voltage magnitude.
        *   - Attached external signals (omega, V_S)
        *
-       * Saturation is included via ksat computed from efdp and SA, SB.
+       * Enabled saturation is included via ksat computed from efdp and SA, SB.
        */
       template <typename scalar_type, typename index_type>
       int Ieeet1<scalar_type, index_type>::initialize()
@@ -360,15 +386,10 @@ namespace GridKit
       {
         using Parameter = typename ModelDataT::Parameters;
 
-        Tr_ = TR_MINIMUM;
         if (data.parameters.contains(Parameter::Tr))
         {
           Tr_ = std::get<RealT>(data.parameters.at(Parameter::Tr));
         }
-        if (Tr_ < TR_MINIMUM)
-        {
-          Tr_ = TR_MINIMUM;
-        }
         if (data.parameters.contains(Parameter::Ka))
         {
           Ka_ = std::get<RealT>(data.parameters.at(Parameter::Ka));
@@ -422,12 +443,33 @@ namespace GridKit
           Ispdlim_ = std::get<RealT>(data.parameters.at(Parameter::Ispdlim));
         }
 
-        // Derived Parameters
-        RealT SR = std::sqrt(Se2_ / Se1_);
+        Tr_ = std::max(Tr_, TIME_CONSTANT_MINIMUM);
+        Ta_ = std::max(Ta_, TIME_CONSTANT_MINIMUM);
+        Te_ = std::max(Te_, TIME_CONSTANT_MINIMUM);
+        Tf_ = std::max(Tf_, TIME_CONSTANT_MINIMUM);
+
+        SA_ = ZERO<RealT>;
+        SB_ = ZERO<RealT>;
+
+        const bool saturation_disabled =
+            Se1_ == ZERO<RealT> && Se2_ == ZERO<RealT>;
+
+        if (saturation_disabled)
+        {
+          return;
+        }
+
+        if (E1_ <= ZERO<RealT> || E2_ <= ZERO<RealT> || E1_ == E2_
+            || Se1_ <= ZERO<RealT> || Se2_ <= ZERO<RealT> || Se1_ == Se2_)
+        {
+          return;
+        }
+
+        const RealT C = std::sqrt(Se2_ / Se1_);
 
         // Solution 1 (Aligned with PW)
-        SA_ = (SR * E1_ - E2_) / (SR - 1);
-        SB_ = Se1_ / (E1_ - SA_) / (E1_ - SA_);
+        SA_ = (C * E1_ - E2_) / (C - ONE<RealT>);
+        SB_ = Se1_ / ((E1_ - SA_) * (E1_ - SA_));
       }
 
       template <typename scalar_type, typename index_type>

GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md

@@ -31,19 +31,47 @@ $S_1$       | [p.u.] | Saturation Parameter                 | 0.04    |
 $S_2$       | [p.u.] | Saturation Parameter                 | 0.33    |
 $I_{\mathrm{spdlim}}$ | [binary] | Speed limit flag indicator       | 0       |
 
+### Parameter Validation
+
+Invalid IEEET1 parameter sets are rejected by the following checks. Let $\epsilon_T=10^{-3}$.
+
+```math
+\begin{aligned}
+  T &\leftarrow \max\!\left(T, \epsilon_T\right)
+    \quad T \in \{T_R, T_A, T_E, T_F\} \\
+  K_A
+    &> 0 \\
+  V_R^{\min}
+    &\le V_R^{\max} \\
+  I_{\mathrm{spdlim}}
+    &\in \{0,1\} \\
+  \left(S_1, S_2\right)
+    &=(0,0)
+      \quad\text{or}\quad
+      \begin{gathered}
+        E_1, E_2, S_1, S_2 > 0 \\
+        E_1 \ne E_2 \\
+        S_1 \ne S_2
+      \end{gathered}
+\end{aligned}
+```
+
 ### Model Derived Parameters
 
-The relationship of the derived parameters is defined by the following quadratic model. The parameters are chosen so that the quadratic model represents
-the expected saturation near the operating region.
-``` math
+When saturation is disabled, $S_A=0$ and $S_B=0$. Otherwise,
+the parameters are chosen so that the following quadratic model represents the
+expected saturation near the operating region:
+
+```math
 \begin{aligned}
   S_1 &= S_B(E_1-S_A)^2 \\
   S_2 &= S_B(E_2-S_A)^2 \\
 \end{aligned}
 ```
 
 Generally, this system has two solutions. The non-extraneous solution is as follows.
-``` math
+
+```math
 \begin{aligned}
   C &=  \sqrt{
    \dfrac
@@ -147,7 +175,14 @@ Here $q$ is GridKit's [Quadratic Ramp](../../../../CommonMath.md#primitives).
 
 ## Initialization
 
-The implementation first applies $T_R \leftarrow \max(T_R, 10^{-3})$. The machine initializes $E_{fd}$ first. IEEET1 reads that value as $E_{fd,0}$, along with any attached $\omega$ and $V_S$, and solves the steady-state algebraic chain so all residuals vanish with $\dot y = 0$. The sensed terminal voltage initializes from the positive bus-voltage magnitude. Saturation and the speed-limit flag are included directly; $V_\text{ref}$ is set to close the $V_{tr}$ equation with the current input values.
+The implementation first applies $T \leftarrow \max(T, 10^{-3})$ for
+$T \in \{T_R, T_A, T_E, T_F\}$. This should be replaced with a structural tempalte change in the future. The machine initializes $E_{fd}$ first. IEEET1
+reads that value as $E_{fd,0}$, along with any attached $\omega$ and $V_S$, and
+solves the steady-state algebraic chain so all residuals vanish with
+$\dot y = 0$. The sensed terminal voltage initializes from the positive
+bus-voltage magnitude. Saturation is included when enabled, and the speed-limit
+flag is included directly; $V_\text{ref}$ is set to close the $V_{tr}$ equation
+with the current input values.
 
 ```math
 \begin{aligned}

tests/UnitTests/PhasorDynamics/ExciterIeeet1Tests.hpp

@@ -74,34 +74,113 @@ namespace GridKit
         return success.report(__func__);
       }
 
-      TestOutcome zeroTrUsesMinimum()
+      TestOutcome zeroTimeConstantsAndDisabledSaturation()
       {
         TestStatus success = true;
 
-        auto data                                                      = makeTestData();
-        data.parameters[PhasorDynamics::Exciter::Ieeet1Parameters::Tr] = 0.0;
+        using Params = PhasorDynamics::Exciter::Ieeet1Parameters;
+
+        auto data                    = makeTestData();
+        data.parameters[Params::Tr]  = 0.0;
+        data.parameters[Params::Ta]  = 0.0;
+        data.parameters[Params::Ke]  = 0.0;
+        data.parameters[Params::Te]  = 0.0;
+        data.parameters[Params::Kf]  = 0.0;
+        data.parameters[Params::Tf]  = 0.0;
+        data.parameters[Params::E1]  = 0.0;
+        data.parameters[Params::E2]  = 0.0;
+        data.parameters[Params::Se1] = 0.0;
+        data.parameters[Params::Se2] = 0.0;
 
         PhasorDynamics::Bus<ScalarT, IdxT>             bus(3.0, 4.0);
         PhasorDynamics::Exciter::Ieeet1<ScalarT, IdxT> exciter(&bus, data);
+        PhasorDynamics::SignalNode<ScalarT, IdxT>      efd_node;
+        ScalarT                                        efd_value{0.0};
+        IdxT                                           efd_index = INVALID_INDEX<IdxT>;
+
+        efd_node.set(&efd_value, &efd_index);
+        exciter.getSignals()
+            .template assignSignalNode<PhasorDynamics::Exciter::Ieeet1InternalVariables::EFD>(&efd_node);
 
         bus.allocate();
         exciter.allocate();
 
         bus.initialize();
-        exciter.initialize();
+        efd_node.init(1.2);
+        success *= (exciter.initialize() == 0);
         exciter.tagDifferentiable();
 
         success *= (exciter.tag()[0]);
+        success *= (exciter.tag()[1]);
+        success *= (exciter.tag()[2]);
+        success *= (exciter.tag()[3]);
+
+        success *= isEqual(exciter.y()[2], static_cast<ScalarT>(1.2));
+        success *= isEqual(exciter.y()[6], static_cast<ScalarT>(0.0));
+        success *= isEqual(exciter.y()[7], static_cast<ScalarT>(1.2));
+        success *= isEqual(exciter.y()[8], static_cast<ScalarT>(0.0));
+
+        for (const auto& value : exciter.y())
+        {
+          success *= std::isfinite(value);
+        }
+
+        exciter.evaluateResidual();
+        for (const auto& value : exciter.getResidual())
+        {
+          success *= std::isfinite(value);
+          success *= isEqual(value, static_cast<ScalarT>(0.0));
+        }
+
+        exciter.y()[2] = 4.0;
+        exciter.evaluateResidual();
+        success        *= isEqual(exciter.getResidual()[8], static_cast<ScalarT>(0.0));
+        exciter.y()[2]  = 1.2;
 
         exciter.yp()[0] = 123.0;
         exciter.evaluateResidual();
-        success *= isEqual(exciter.getResidual()[0], static_cast<ScalarT>(-123.0));
+        success         *= isEqual(exciter.getResidual()[0], static_cast<ScalarT>(-123.0));
+        exciter.yp()[0]  = 0.0;
 
-        exciter.yp()[0] = 0.0;
-        exciter.y()[0]  = 4.0;
+        exciter.y()[0] = 4.0;
         exciter.evaluateResidual();
         success *= isEqual(exciter.getResidual()[0], static_cast<ScalarT>(1.0e3));
 
+        exciter.y()[0] = 5.0;
+        exciter.y()[4] = 0.02;
+        exciter.evaluateResidual();
+        success *= isEqual(exciter.getResidual()[1], static_cast<ScalarT>(1.0e3));
+
+        exciter.y()[4] = 0.0;
+        exciter.y()[1] = 1.0;
+        exciter.evaluateResidual();
+        success *= isEqual(exciter.getResidual()[2], static_cast<ScalarT>(1.0e3));
+
+        exciter.y()[1] = 0.0;
+        exciter.y()[5] = 1.0;
+        exciter.evaluateResidual();
+        success *= isEqual(exciter.getResidual()[3], static_cast<ScalarT>(1.0e3));
+
+        return success.report(__func__);
+      }
+
+      TestOutcome invalidSaturationParameters()
+      {
+        TestStatus success = true;
+
+        using Params = PhasorDynamics::Exciter::Ieeet1Parameters;
+
+        auto data                    = makeTestData();
+        data.parameters[Params::Se1] = 0.0;
+
+        PhasorDynamics::Bus<ScalarT, IdxT>             bus(3.0, 4.0);
+        PhasorDynamics::Exciter::Ieeet1<ScalarT, IdxT> exciter(&bus, data);
+
+        bus.allocate();
+        exciter.allocate();
+
+        success *= (exciter.verify() != 0);
+
         return success.report(__func__);
       }
 

tests/UnitTests/PhasorDynamics/runExciterIeeet1Tests.cpp

@@ -8,7 +8,8 @@ int main()
 
   result += test.constructor();
   result += test.zeroInitialResidual();
-  result += test.zeroTrUsesMinimum();
+  result += test.zeroTimeConstantsAndDisabledSaturation();
+  result += test.invalidSaturationParameters();
 #ifdef GRIDKIT_ENABLE_ENZYME
   result += test.jacobian();
 #endif