Modify Smooth Approximation of Limits

CHANGELOG.md

@@ -38,6 +38,7 @@
 - Updated Jacobian value storage from `ScalarT` to `RealT`.
 - Added a header file defining constants to be used throughout the code.
 - Added `GridKitDocs` target for Doxygen documentation.
+- Added a header file defining common math functions (e.g., sigmoid) to be used throughout the code.
 
 ## v0.1
 

GridKit/CMakeLists.txt

@@ -28,6 +28,7 @@ add_subdirectory(Solver)
 install(
   FILES
     Constants.hpp
+    CommonMath.hpp
     ScalarTraits.hpp
   DESTINATION
     include/GridKit)

GridKit/CommonMath.hpp

@@ -0,0 +1,96 @@
+#pragma once
+
+#include <cmath>
+
+#include <GridKit/Constants.hpp>
+#include <GridKit/ScalarTraits.hpp>
+
+namespace GridKit
+{
+  namespace Math
+  {
+    /**
+     * @brief Scaled sigmoid activation function
+     *
+     * @note The sigmoid constant (mu) value is chosen to balance accuracy
+     * and finite derivatives. Large values more closely approximate a step
+     * function, but lead to inf or NaN derivatives.
+     * The value of 240 corresponds to a time step of 1/4 of a 60Hz cycle.
+     * From experiments, numerical derivatives break down for mu above 700.
+     *
+     * @tparam ScalarT - scalar data type
+     *
+     * @param[in] x
+     * @return value of the sigmoid function
+     */
+    template <class ScalarT>
+    __attribute__((always_inline)) inline ScalarT sigmoid(const ScalarT x)
+    {
+      using RealT               = typename GridKit::ScalarTraits<ScalarT>::RealT;
+      static constexpr RealT MU = 240.0;
+      return ONE<RealT> / (ONE<RealT> + std::exp(-MU * x));
+    }
+
+    /**
+     * @brief Low indicator function for regulator limits
+     *
+     * @tparam ScalarT - Scalar data type
+     * @tparam RealT - Real data type (see GridKit::ScalarTraits<ScalarT>::RealT)
+     *
+     * @param[in] limit_min - Minimum limit
+     * @param[in] x - State variable
+     * @param[in] f - Conditional derivative of state variable
+     * @return Scalar value indicating limit activation
+     */
+    template <class ScalarT, typename RealT>
+    __attribute__((always_inline)) inline ScalarT indicator_low(
+        const RealT   limit_min,
+        const ScalarT x,
+        const ScalarT f)
+    {
+      return sigmoid(limit_min - x) * sigmoid(-f);
+    }
+
+    /**
+     * @brief High indicator function for regulator limits
+     *
+     * @tparam ScalarT - Scalar data type
+     * @tparam RealT - Real data type (see GridKit::ScalarTraits<ScalarT>::RealT)
+     *
+     * @param[in] limit_max - Maximum limit
+     * @param[in] x - State variable
+     * @param[in] f - Conditional derivative of state variable
+     * @return Scalar value indicating limit activation
+     */
+    template <class ScalarT, typename RealT>
+    __attribute__((always_inline)) inline ScalarT indicator_high(
+        const RealT   limit_max,
+        const ScalarT x,
+        const ScalarT f)
+    {
+      return sigmoid(x - limit_max) * sigmoid(f);
+    }
+
+    /**
+     * @brief Net Indicator function for regulator limits
+     *
+     * @tparam ScalarT - Scalar data type
+     * @tparam RealT - Real data type (see GridKit::ScalarTraits<ScalarT>::RealT)
+     *
+     * @param[in] limit_min - Minimum limit
+     * @param[in] limit_max - Maximum limit
+     * @param[in] x - State variable
+     * @param[in] f - Conditional derivative of state variable
+     * @return Scalar value indicating limit activation
+     */
+    template <class ScalarT, typename RealT>
+    __attribute__((always_inline)) inline ScalarT indicator(
+        const RealT   limit_min,
+        const RealT   limit_max,
+        const ScalarT x,
+        const ScalarT f)
+    {
+      return (ONE<RealT> - indicator_low(limit_min, x, f)) * (ONE<RealT> - indicator_high(limit_max, x, f));
+    }
+  } // namespace Math
+} // namespace GridKit

GridKit/Model/Evaluator.hpp

@@ -2,6 +2,7 @@
 
 #include <vector>
 
+#include <GridKit/CommonMath.hpp>
 #include <GridKit/Constants.hpp>
 #include <GridKit/LinearAlgebra/SparseMatrix/COO_Matrix.hpp>
 #include <GridKit/ScalarTraits.hpp>

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

@@ -150,15 +150,6 @@ namespace GridKit
         /// Component signal extension
         ComponentSignals<ScalarT, IdxT, Ieeet1InternalVariables, Ieeet1ExternalVariables> signals_;
 
-        // Scale of Sigmoid function
-        // (temporary local implementation)
-        // This value gave higher precision.
-        const RealT mu_{400000.0};
-
-        // Activation function (sigmoid approximation)
-        ScalarT sigmoid(ScalarT x);
-        ScalarT indicator(ScalarT x, ScalarT f);
-
         // Parameter initialization function
         void initModelParams(const model_data_type& data);
       };

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

@@ -236,48 +236,6 @@ namespace GridKit
         return 0;
       }
 
-      /**
-       * @brief  Scaled sigmoid activation function
-       *
-       * @param[in] x State variable
-       * @tparam ScalarT Scalar data type
-       * @tparam IdxT Index data type
-       * @return Sigmoid approximation evaluated at x
-       *
-       * @warning This needs to be abstracted to be used
-       *          across phasor dynamics. Identical pattern is
-       *          being used in TGOV1 model.
-       *
-       *   Algebraic approximation of transcendental sigmoid.
-       */
-      template <class ScalarT, typename IdxT>
-      ScalarT Ieeet1<ScalarT, IdxT>::sigmoid(ScalarT x)
-      {
-        return ((HALF<RealT> * mu_ * x) / (ONE<RealT> + std::abs(mu_ * x))) + HALF<RealT>;
-      }
-
-      /**
-       * @brief Net Indicator function for regulator limits
-       *
-       * @param[in] x State variable
-       * @param[in] f Conditional derivative of state variable
-       * @tparam ScalarT Scalar data type
-       * @tparam IdxT Index data type
-       * @return Scalar value indicating limit activation.
-       *
-       * @warning This needs to be abstracted to be used
-       *          across phasor dynamics. Identical pattern is
-       *          being used in TGOV1 model.
-       */
-      template <class ScalarT, typename IdxT>
-      ScalarT Ieeet1<ScalarT, IdxT>::indicator(ScalarT x, ScalarT f)
-      {
-
-        ScalarT ind_low  = (this->sigmoid(Vrmin_ - x)) * (this->sigmoid(-f));
-        ScalarT ind_high = (this->sigmoid(x - Vrmax_)) * (this->sigmoid(f));
-        return (ONE<RealT> - ind_low) * (ONE<RealT> - ind_high);
-      }
-
       /**
        * @brief Residuals of system equations
        *
@@ -326,7 +284,7 @@ namespace GridKit
 
         // The 'pre-limit' derivative of Pv
         ScalarT f      = -vr + Ka_ * vtr;
-        ScalarT vr_ind = this->indicator(vr, f);
+        ScalarT vr_ind = Math::indicator(Vrmin_, Vrmax_, vr, f);
 
         // Internal Differential Equations
         f_[0] = -vts_dot + (Ec_ - vts) / Tr_;
@@ -340,11 +298,8 @@ namespace GridKit
         f_[6] = -ve + ksat * efdp;
         f_[7] = -efd + efdp + omega * efdp * Ispdlim_;
 
-        // TODO smooth approaximation for Enzyme
-        // NOTE seems about double PW saturation.
-        f_[8] = -ksat;
-        if (efdp > SA_)
-          f_[8] += SB_ * (efdp - SA_) * (efdp - SA_);
+        ScalarT efd_sat = (efdp - SA_) * (Math::sigmoid(efdp - SA_));
+        f_[8]           = -ksat + SB_ * efd_sat * efd_sat;
 
         return 0;
       }

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

@@ -155,58 +155,11 @@ The indicator function $\phi$ can be defined in terms of a scaled activation fun
 \end{aligned}
 ```
 
-The scale of the sigmoid function ($\alpha$ on the order of $10^3$) should be chosen so that for all practical parameters of the IEEET1 model, the sigmoid acts as a step function. This is further approximated by an algebraic form to obtain a practical function during implementation.
+The scale of the sigmoid function ($\alpha$ on the order of $10^3$) should be chosen so that for all practical parameters of the IEEET1 model, the sigmoid acts as a step function.
 ```math
 \begin{aligned}
    \sigma(x) = 
       \dfrac{1}{1+\exp(-\alpha x)}
-   \approx 
-      \dfrac{1}{2}
-      \left(
-         \dfrac{\alpha x}{1+|\alpha x|} + 1
-      \right)
-\end{aligned}
-```
-
-Applying these approximations produces a smooth approximation of the explicit differential equation of $V_R$. The approximation written out below approaches an exact solution as $\alpha\to\infty$.
-
-<!--
-Very explicit version, if prefered
-```math
-\begin{aligned}
-   \dot{V}_R 
-   &\approx \dfrac{K_{a}V_{tr}-V_{R}}{T_A}
-      \left[
-         1 +
-         \dfrac{\alpha^2(V_{rmin}\text{-}V_R)(K_{a}V_{tr}-V_{R})}
-         {4(1+\alpha |V_{rmin}\text{-}V_R|)(T_A+\alpha |K_{a}V_{tr}-V_{R}|)}
-      \right]
-      \left[
-         1-
-            \dfrac{\alpha^2(V_R\text{-}V_{rmax})(K_{a}V_{tr}-V_{R})}
-            {4(1+\alpha| V_R\text{-}V_{rmax}|)(T_A+\alpha |K_{a}V_{tr}-V_{R}|)}
-      \right]\\
-\end{aligned}
-```
--->
-```math
-\begin{aligned}
-   \dot{V}_R 
-   &\approx f
-      \left[
-         1 +
-         \dfrac{\alpha^2}{4}
-         \dfrac{V_{rmin}\text{-}V_R}
-         {1+\alpha |V_{rmin}\text{-}V_R|}
-         \dfrac{f}{1+\alpha |f|}
-      \right]
-      \left[
-         1-
-         \dfrac{\alpha^2}{4}
-         \dfrac{V_R\text{-}V_{rmax}}
-            {1+\alpha| V_R\text{-}V_{rmax}|}
-            \dfrac{f}{1+\alpha |f|}
-      \right]\\
 \end{aligned}
 ```
 
@@ -228,23 +181,20 @@ The algebraic equations of the exciter.
    \end{cases} \\
 \end{aligned}
 ```
-
 #### Smooth Piecewise Approximation (Algebraic) 
 
 For the algebraic piecewise functions (non-flags), this implementation is straightforward when the approximation above is used.
 ```math
 \begin{aligned}
+    E_{fd}
+    &=(1 + \omega I_{spdlm})E_{fd}' \\ 
     k_{sat}
-    &=\sigma (E_{fd}' -S_A) \cdot S_B(E_{fd}' -S_A)^2   
-\end{aligned}
-```
-The approximation written out below approaches an exact solution as $\alpha\to\infty$.
-```math
-\begin{aligned}
-    k_{sat}&\approx\dfrac{\alpha S_B}{2} \dfrac{(E_{fd}' -S_A)^3}{1+\alpha|E_{fd}' -S_A|}
+    &=S_B\left[(E_{fd}' -S_A) \cdot \sigma (E_{fd}' -S_A)\right]^2   
 \end{aligned}
 ```
 
+The approximation approaches an exact solution as $\alpha\to\infty$.
+
 ## Initialization
 
 At steady state we assume that $V_R$ is at or within its limits, 

GridKit/Model/PhasorDynamics/Governor/Tgov1/README.md

@@ -101,10 +101,10 @@ The indicator function $\phi$ can be defined in terms of a scaled activation fun
 \end{aligned}
 ```
 
-The scale of the sigmoid function ($\alpha$ on the order of $10^3$) should be chosen so that for all practical parameters of the TGOV1 model, the sigmoid acts as a step function. This is further approximated by an algebraic form to obtain a practical function during implementation.
+The scale of the sigmoid function ($\alpha$ on the order of $10^3$) should be chosen so that for all practical parameters of the TGOV1 model, the sigmoid acts as a step function.
 ```math
 \begin{aligned}
-   \sigma(x) =\dfrac{1}{1+\exp(-\alpha x)}\approx \dfrac{1}{2}\left(\dfrac{\alpha x}{1+|\alpha x|}\right)
+   \sigma(x) =\dfrac{1}{1+\exp(-\alpha x)}
 \end{aligned}
 ```
 

GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp

@@ -121,20 +121,10 @@ namespace GridKit
         // Input States (which can be parameters)
         ScalarT pref_{0};
 
-        // Scale of Sigmoid function (temporary local implementation)
-        const ScalarT mu_{4000.0};
-
         /// Component signal extension
         ComponentSignals<ScalarT, IdxT, Tgov1InternalVariables, Tgov1ExternalVariables> signals_;
 
-        // Activation function (sigmoid approximation)
-        ScalarT sigmoid(ScalarT x);
-
-        // Indicator of Valve limit states
-        ScalarT indicator_low(ScalarT x, ScalarT f);
-        ScalarT indicator_high(ScalarT x, ScalarT f);
-        ScalarT indicator(ScalarT x, ScalarT f);
-
+        // Parameter initialization function
         void initializeParameters(const model_data_type& data);
 
         /* Local copies of external variables */

GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp

@@ -229,48 +229,6 @@ namespace GridKit
         return 0;
       }
 
-      /**
-       * @brief Scaled sigmoid activation function
-       *
-       * Temporary local implementation of smooth approximation
-       * of a piecewise differential equation. Ideally this is
-       * a more abstracted capability with GK.
-       *
-       * Algebraic approximation of transcendental sigmoid.
-       */
-      template <class ScalarT, typename IdxT>
-      ScalarT Tgov1<ScalarT, IdxT>::sigmoid(ScalarT x)
-      {
-        return ((HALF<RealT> * mu_ * x) / (ONE<RealT> + std::abs(mu_ * x))) + HALF<RealT>;
-      }
-
-      /**
-       * @brief Indicator function for lower valve limit violation.
-       */
-      template <class ScalarT, typename IdxT>
-      ScalarT Tgov1<ScalarT, IdxT>::indicator_low(ScalarT x, ScalarT f)
-      {
-        return (this->sigmoid(Pvmin_ - x)) * (this->sigmoid(-f));
-      }
-
-      /**
-       * @brief Indicator function for high valve limit violation.
-       */
-      template <class ScalarT, typename IdxT>
-      ScalarT Tgov1<ScalarT, IdxT>::indicator_high(ScalarT x, ScalarT f)
-      {
-        return (this->sigmoid(x - Pvmax_)) * (this->sigmoid(f));
-      }
-
-      /**
-       * @brief Net Indicator function for valve limits.
-       */
-      template <class ScalarT, typename IdxT>
-      ScalarT Tgov1<ScalarT, IdxT>::indicator(ScalarT x, ScalarT f)
-      {
-        return (ONE<RealT> - this->indicator_low(x, f)) * (ONE<RealT> - this->indicator_high(x, f));
-      }
-
       /**
        * @brief Internal residuals
        *
@@ -297,7 +255,7 @@ namespace GridKit
 
         // The 'pre-limit' derivative of Pv
         ScalarT func     = (-pv + (pref_ - omega) / R_) / T1_;
-        ScalarT valv_ind = this->indicator(pv, func);
+        ScalarT valv_ind = Math::indicator(Pvmin_, Pvmax_, pv, func);
 
         // Internal Differential Equations
         f[0] = -ptx_dot + pv - (ptx + T2_ * pv) / T3_;