Add REECB phasor dynamics model

Author: lukelowry

CHANGELOG.md

@@ -62,6 +62,7 @@
 - Added cmake-format hooks, including in pre-commit.
 - Added off-nominal tap ratio and phase shift support to the PhasorDynamics `Branch` model.
 - Added portable Vector class to GridKit
+- Added `REECB` converter model implementation for PhasorDynamics.
 
 ## v0.1
 

GridKit/CommonMath.hpp

@@ -15,7 +15,7 @@ namespace GridKit
      *
      * @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.
+     * function, but can make the transition numerically stiff.
      *
      * @tparam ScalarT - scalar data type
      *
@@ -27,7 +27,7 @@ namespace GridKit
     {
       using RealT               = typename GridKit::ScalarTraits<ScalarT>::RealT;
       static constexpr RealT MU = 240.0;
-      return ONE<RealT> / (ONE<RealT> + std::exp(-MU * x));
+      return HALF<RealT> * (ONE<RealT> + std::tanh(HALF<RealT> * MU * x));
     }
 
     /**
@@ -153,7 +153,32 @@ namespace GridKit
     }
 
     /**
-     * @brief Smooth two-sided deadband function
+     * @brief Smooth Type 1 no-offset two-sided deadband function
+     *
+     * Smooth approximation to a deadband that returns zero inside the band and
+     * passes the input through unchanged outside the band.
+     *
+     * @tparam ScalarT - scalar data type
+     * @tparam LowerT - data type of the lower limit
+     * @tparam UpperT - data type of the upper limit
+     *
+     * @param[in] x - Input signal
+     * @param[in] lower - Lower breakpoint
+     * @param[in] upper - Upper breakpoint
+     * @return Smooth no-offset deadbanded value
+     */
+    template <class ScalarT, typename RealT>
+    __attribute__((always_inline)) inline ScalarT deadband1(
+        const ScalarT x,
+        const RealT   lower,
+        const RealT   upper)
+    {
+      assert(lower <= upper);
+      return x * (sigmoid(lower - x) + sigmoid(x - upper));
+    }
+
+    /**
+     * @brief Smooth Type 2 offset two-sided deadband function
      *
      * Smooth approximation to x - min(max(x, lower), upper), composed from the
      * smooth ramp function.
@@ -164,10 +189,10 @@ namespace GridKit
      * @param[in] x - Input signal
      * @param[in] lower - Lower breakpoint
      * @param[in] upper - Upper breakpoint
-     * @return Smooth deadbanded value
+     * @return Smooth offset deadbanded value
      */
     template <class ScalarT, typename RealT>
-    __attribute__((always_inline)) inline ScalarT deadband(
+    __attribute__((always_inline)) inline ScalarT deadband2(
         const ScalarT x,
         const RealT   lower,
         const RealT   upper)
@@ -331,13 +356,15 @@ namespace GridKit
      * @return Scalar value in [0, 1]: 1 when dynamics should pass through,
      *         0 when integration should be blocked.
      */
-    template <class ScalarT, typename RealT>
+    template <class ScalarT, typename LowerT, typename UpperT>
     __attribute__((always_inline)) inline ScalarT indicator(
         const ScalarT x,
         const ScalarT f,
-        const RealT   limit_min,
-        const RealT   limit_max)
+        const LowerT  limit_min,
+        const UpperT  limit_max)
     {
+      using RealT = typename GridKit::ScalarTraits<ScalarT>::RealT;
+
       assert(limit_min <= limit_max);
 
       ScalarT above_min = above(x, limit_min);
@@ -357,20 +384,21 @@ namespace GridKit
      * and blocks motion that would push further into saturation.
      *
      * @tparam ScalarT - Scalar data type
-     * @tparam RealT - Real data type (see GridKit::ScalarTraits<ScalarT>::RealT)
+     * @tparam LowerT - data type of the lower limit
+     * @tparam UpperT - data type of the upper limit
      *
      * @param[in] x - Limited state or limited output signal
      * @param[in] f - Pre-limit derivative
      * @param[in] limit_min - Minimum limit
      * @param[in] limit_max - Maximum limit
      * @return Smooth anti-windup limited derivative
      */
-    template <class ScalarT, typename RealT>
+    template <class ScalarT, typename LowerT, typename UpperT>
     __attribute__((always_inline)) inline ScalarT antiwindup(
         const ScalarT x,
         const ScalarT f,
-        const RealT   limit_min,
-        const RealT   limit_max)
+        const LowerT  limit_min,
+        const UpperT  limit_max)
     {
       return indicator(x, f, limit_min, limit_max) * f;
     }

GridKit/CommonMath.md

@@ -49,7 +49,8 @@ The scale $\mu=4\cdot f_{\text{sync}}=240$ is chosen so $\sigma$ behaves like a
 | `max` | Smooth binary maximum | `REECA` |
 | `min` | Smooth binary minimum | `REECA` |
 | `clamp` | Bounded saturation | `IEEEST` |
-| `deadband` | Signed two-sided deadband | `REECA` |
+| `deadband1` | Type 1 no-offset signed two-sided deadband | - |
+| `deadband2` | Type 2 offset signed two-sided deadband | `REECA` |
 | `slew` | Symmetric slew-rate limiter | - |
 | `linseg` | Saturated linear segment contribution | `REGCA`, `REECA` |
 | `above` | Above-lower-limit indicator | - |
@@ -84,7 +85,15 @@ x & \ell \le x \le u \\
 u & x > u
 \end{cases}
 \\
-\text{deadband}(x;\ell,u)
+\text{deadband1}(x;\ell,u)
+&=
+\begin{cases}
+x & x < \ell \\
+0 & \ell \le x \le u \\
+x & x > u
+\end{cases}
+\\
+\text{deadband2}(x;\ell,u)
 &=
 \begin{cases}
 x-\ell & x < \ell \\
@@ -156,7 +165,8 @@ f & x \ge u \land f < 0 \\
 \text{max}(x,y) &= y+\rho(x-y) \\
 \text{min}(x,y) &= x-\rho(x-y) \\
 \text{clamp}(x;\ell,u) &= \ell+\rho(x-\ell)-\rho(x-u) \\
-\text{deadband}(x;\ell,u) &= \rho(x-u)-\rho(\ell-x) \\
+\text{deadband1}(x;\ell,u) &= x\left[\sigma(\ell-x)+\sigma(x-u)\right] \\
+\text{deadband2}(x;\ell,u) &= \rho(x-u)-\rho(\ell-x) \\
 \text{slew}(f;r) &= -r+\rho(f+r)-\rho(f-r) \\
 \text{linseg}(x;a,b,h) &= \dfrac{h}{b-a}\left[\rho(x-a)-\rho(x-b)\right] \\
 \text{above}(x;\ell) &= \sigma(x-\ell) \\

GridKit/Model/PhasorDynamics/CMakeLists.txt

@@ -40,6 +40,7 @@ add_subdirectory(Branch)
 add_subdirectory(Bus)
 add_subdirectory(BusFault)
 add_subdirectory(BusToSignalAdapter)
+add_subdirectory(Converter)
 add_subdirectory(Exciter)
 add_subdirectory(Governor)
 add_subdirectory(Load)

GridKit/Model/PhasorDynamics/ComponentLibrary.hpp

@@ -4,6 +4,7 @@
 #include <GridKit/Model/PhasorDynamics/Bus/Bus.hpp>
 #include <GridKit/Model/PhasorDynamics/BusFault/BusFault.hpp>
 #include <GridKit/Model/PhasorDynamics/BusToSignalAdapter/BusToSignalAdapter.hpp>
+#include <GridKit/Model/PhasorDynamics/Converter/REECB/Reecb.hpp>
 #include <GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp>
 #include <GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPti.hpp>
 #include <GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp>

GridKit/Model/PhasorDynamics/Converter/CMakeLists.txt

@@ -0,0 +1,6 @@
+# [[
+#  Author(s):
+#    - Luke Lowery <lukel@tamu.edu>
+# ]]
+
+add_subdirectory(REECB)

GridKit/Model/PhasorDynamics/Converter/README.md

@@ -9,6 +9,6 @@ models and the bus equations, typically through commanded active and reactive cu
 
 The GridKit converter documentation includes:
 
-- Renewable Energy Generator/Converter Model REGCA (See [REGCA](REGCA/README.md))
 - Renewable Energy Generator/Converter Model REGCB (See [REGCB](REGCB/README.md))
 - Renewable Energy Electrical Control Model REECA (See [REECA](REECA/README.md))
+- Renewable Energy Electrical Control Model REECB (See [REECB](REECB/README.md))

GridKit/Model/PhasorDynamics/Converter/REECB/CMakeLists.txt

@@ -0,0 +1,54 @@
+# [[
+#  Author(s):
+#    - Luke Lowery <lukel@tamu.edu>
+# ]]
+
+set(_install_headers Reecb.hpp ReecbData.hpp)
+
+if(GRIDKIT_ENABLE_ENZYME)
+  gridkit_add_library(
+    phasor_dynamics_converter_reecb
+    SOURCES ReecbEnzyme.cpp
+    HEADERS ${_install_headers}
+    INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include
+    LINK_LIBRARIES
+      PUBLIC
+      GridKit::phasor_dynamics_core
+      PUBLIC
+      GridKit::phasor_dynamics_signal
+      PRIVATE
+      ClangEnzymeFlags
+    COMPILE_OPTIONS
+      PRIVATE
+      -mllvm
+      -enzyme-auto-sparsity=1
+      -fno-math-errno)
+else()
+  gridkit_add_library(
+    phasor_dynamics_converter_reecb
+    SOURCES Reecb.cpp
+    HEADERS ${_install_headers}
+    INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include
+    LINK_LIBRARIES
+      PUBLIC
+      GridKit::phasor_dynamics_core
+      PUBLIC
+      GridKit::phasor_dynamics_signal)
+endif()
+
+gridkit_add_library(
+  phasor_dynamics_converter_reecb_dependency_tracking
+  SOURCES ReecbDependencyTracking.cpp
+  INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include
+  LINK_LIBRARIES
+    PUBLIC
+    GridKit::phasor_dynamics_core
+    PUBLIC
+    GridKit::phasor_dynamics_signal_dependency_tracking)
+
+target_link_libraries(
+  phasor_dynamics_components
+  INTERFACE GridKit::phasor_dynamics_converter_reecb)
+target_link_libraries(
+  phasor_dynamics_components_dependency_tracking
+  INTERFACE GridKit::phasor_dynamics_converter_reecb_dependency_tracking)