Minimal corrections to Phasor Dynamics models for validation

GridKit/Model/PhasorDynamics/Converter/REGCA/README.md

@@ -36,9 +36,9 @@ $Q_0$                            | [p.u.]   | `q0`     | Initial reactive power
 $S^\mathrm{base}$                | [MVA]    | `mva`    | REGCA component power base                            | 100.0         |
 $T_\mathrm{g}$                   | [sec]    | `Tg`     | Converter current-control lag time constant           | 0.02          | Block name: `Tg`
 $T_M$                            | [sec]    | `TM`     | Terminal voltage sensor time constant                 | 0.02          | Block name: `Tfltr`
-$R_q^{\max}$                     | [p.u./s] | `Rqmax`  | Reactive-current recovery positive rate limit         | 999.0         | Block name: `Iqrmax`
-$R_q^{\min}$                     | [p.u./s] | `Rqmin`  | Reactive-current recovery negative rate limit         | -999.0        | Block name: `Iqrmin`
-$R_p^{\max}$                     | [p.u./s] | `Rpmax`  | Active-current magnitude recovery rate limit          | 999.0         | Block name: `rrpwr`
+$R_q^{\max}$                     | [p.u./s] | `Rqmax`  | Reactive-current recovery positive rate limit         | 999.0         | Block name: `Iqrmax`; disabled when $R_q^{\max}\le 0$
+$R_q^{\min}$                     | [p.u./s] | `Rqmin`  | Reactive-current recovery negative rate limit         | -999.0        | Block name: `Iqrmin`; disabled when $R_q^{\min}\ge 0$
+$R_p^{\max}$                     | [p.u./s] | `Rpmax`  | Active-current magnitude recovery rate limit          | 999.0         | Block name: `rrpwr`; must be nonnegative
 $s_L$                            | [binary] | `sL`     | LVPL switch                                           | 1             | Block name: `LPVLSW`
 $I_{L1}$                         | [p.u.]   | `IL1`    | LVPL upper-current ceiling                            | 1.1           | Block name: `LVPL1`
 $V_{L0}$                         | [p.u.]   | `VL0`    | LVPL zero-crossing voltage                            | 0.4           | Block name: `zerox`
@@ -65,13 +65,11 @@ every other condition is a configuration error.
   S^\mathrm{base}
     &> 0 \\
   R_p^{\max}
-    &> 0 \\
-  R_q^{\min}
-    &< 0 < R_q^{\max} \\
-  s_L
-    &\in \{0,1\} \\
+    &\ge 0 \\
   I_{L1}
     &\ge 0 \\
+  s_L
+    &\in \{0,1\} \\
   0
     &\le V_{L0} < V_{L1} \\
   0
@@ -171,16 +169,19 @@ f_\mathrm{p}^{\lim}
 
 ### Differential Equations
 
-The $I_q$ limiter branch is selected by the initial reactive power $Q_0$.
+The $I_q$ limiter branch is selected by the initial reactive power $Q_0$ and
+the sign that enables the corresponding limit.
 
 ```math
 \begin{aligned}
   0 &= -\dot V_M + \dfrac{1}{T_M} (V_T - V_M) \\
   0 &= -\dot I_q +
     \begin{cases}
-      \text{min}(f_\mathrm{q}, R_q^{\max}) & Q_0 > 0 \\
-      f_\mathrm{q}                         & Q_0 = 0 \\
-      \text{max}(f_\mathrm{q}, R_q^{\min}) & Q_0 < 0
+      \text{min}(f_\mathrm{q}, R_q^{\max})
+        & Q_0 > 0 \land R_q^{\max} > 0 \\
+      \text{max}(f_\mathrm{q}, R_q^{\min})
+        & Q_0 < 0 \land R_q^{\min} < 0 \\
+      f_\mathrm{q} & \text{otherwise}
     \end{cases} \\
   0 &= -\dot I_p +
     \begin{cases}

GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp

@@ -89,11 +89,11 @@ namespace GridKit
 
         use_rqmax_ = ZERO<RealT>;
         use_rqmin_ = ZERO<RealT>;
-        if (q0_ > ZERO<RealT>)
+        if (q0_ > ZERO<RealT> && Rqmax_ > ZERO<RealT>)
         {
           use_rqmax_ = ONE<RealT>;
         }
-        else if (q0_ < ZERO<RealT>)
+        else if (q0_ < ZERO<RealT> && Rqmin_ < ZERO<RealT>)
         {
           use_rqmin_ = ONE<RealT>;
         }
@@ -436,8 +436,7 @@ namespace GridKit
         }
 
         check(mva_base_ > ZERO<RealT>, "mva must be positive");
-        check(Rpmax_ > ZERO<RealT>, "Rpmax must be positive");
-        check(Rqmin_ < ZERO<RealT> && ZERO<RealT> < Rqmax_, "Rqmin < 0 < Rqmax is required");
+        check(Rpmax_ >= ZERO<RealT>, "Rpmax must be non-negative");
         check(IL1_ >= ZERO<RealT>, "IL1 must be non-negative");
         check(KL_ > ZERO<RealT>, "LVPL release slope must be positive");
         check(ZERO<RealT> <= VL0_ && VL0_ < VL1_, "VL0/VL1 must satisfy 0 <= VL0 < VL1");

GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp

@@ -38,7 +38,7 @@ namespace GridKit
         EV,      ///< \f$e_V\f$ Algebraic voltage-error summing output [p.u.]
         VLL,     ///< \f$V_{\mathrm{LL}}\f$ Algebraic input lead-lag output [p.u.]
         VHV,     ///< \f$V_{\mathrm{HV}}\f$ Algebraic high-value gate output [p.u.]
-        SE,      ///< \f$S_E\f$ Algebraic exciter saturation coefficient [p.u.]
+        SE,      ///< \f$s_e\f$ Scaled-quadratic saturation contribution [p.u.]
         VFE,     ///< \f$V_{\mathrm{FE}}\f$ Algebraic exciter feedback drive [p.u.]
         EFD,     ///< \f$E_{\mathrm{fd}}\f$ Algebraic field-voltage output [p.u.]
         MAXIMUM, ///< Number of ESDC1A internal variables
@@ -166,6 +166,7 @@ namespace GridKit
         RealT lim_on_{1};
         RealT SA_{0};
         RealT SB_{0};
+        RealT Ke_eff_{Ke_};
 
         IdxT parameter_error_count_{0};
 

GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aData.hpp

@@ -69,7 +69,7 @@ namespace GridKit
         vc,  ///< \f$V_C\f$ Filtered terminal-voltage magnitude [p.u.]
         vr,  ///< \f$V_R\f$ Voltage-regulator output [p.u.]
         vf,  ///< \f$V_F\f$ Stabilizing feedback state [p.u.]
-        se,  ///< \f$S_E\f$ Exciter saturation coefficient [p.u.]
+        se,  ///< \f$s_e\f$ Scaled-quadratic saturation contribution [p.u.]
         vfe  ///< \f$V_{\mathrm{FE}}\f$ Exciter feedback drive [p.u.]
       };
 

GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp

@@ -162,7 +162,6 @@ namespace GridKit
         }
 
         check(Ka_ > ZERO<RealT>, "Ka must be positive");
-        check(Tc_ >= ZERO<RealT>, "Tc must be non-negative");
         check(Vrmin_ <= Vrmax_, "Vrmin must be less than or equal to Vrmax");
         check(UEL_ >= static_cast<IdxT>(0) && UEL_ <= static_cast<IdxT>(3),
               "UEL must be 0, 1, 2, or 3");
@@ -176,14 +175,11 @@ namespace GridKit
         {
           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");
-
-          const bool saturation_points_are_ordered =
-              (E2_ > E1_ && Se2_ > Se1_)
-              || (E2_ < E1_ && Se2_ < Se1_);
-          check(saturation_points_are_ordered,
-                "E1/E2 and Se1/Se2 must be ordered consistently");
+          check(Se1_ >= ZERO<RealT>, "Se1 must be non-negative when saturation is enabled");
+          check(Se2_ >= ZERO<RealT>, "Se2 must be non-negative when saturation is enabled");
+
+          const bool sat_ordered = (E2_ > E1_ && Se2_ > Se1_) || (E2_ < E1_ && Se2_ < Se1_);
+          check(sat_ordered, "E1/E2 and Se1/Se2 must be ordered consistently");
         }
 
         if (!signals_.template isAssigned<Esdc1aInternalVariables::EFD>())
@@ -226,6 +222,16 @@ namespace GridKit
        * preserving the seeded `efd`, latches attached Known inputs, and
        * publishes the reference to an attached `vref` signal.
        *
+       * @warning IEEE Std 421.5-2016 states: “In some programs, if
+       *          \f$K_{E}\f$ is entered as zero, \f$K_{E}\f$ is automatically
+       *          calculated by the program to represent a self-excited shunt
+       *          field and a trimmed rheostat as its initial condition.” GridKit
+       *          preserves the configured \f$K_{E}\f$ and resolves
+       *          \f$K_{E}^{\mathrm{eff}}\f$ using the PSS/E-compatible
+       *          \f$V_R = V_R^{\max}/10 = 0.1 V_R^{\max}\f$ rule. The divisor
+       *          10 is unitless and represents 10% of the maximum regulator
+       *          output.
+       *
        * @return Zero on success; nonzero when the configuration or operating point is rejected.
        */
       template <typename scalar_type, typename index_type>
@@ -314,8 +320,19 @@ namespace GridKit
           return 1;
         }
 
-        const ScalarT se0  = SB_ * Math::qramp(efdp0 - SA_);
-        const ScalarT vfe0 = (Ke_ + se0) * efdp0;
+        const ScalarT se0 = SB_ * Math::qramp(efdp0 - SA_);
+
+        if (Ke_ == ZERO<RealT>)
+        {
+          Ke_eff_ = (Vrmax_ / 10.0 - static_cast<RealT>(se0))
+                    / static_cast<RealT>(efdp0);
+        }
+        else
+        {
+          Ke_eff_ = Ke_;
+        }
+
+        const ScalarT vfe0 = Ke_eff_ * efdp0 + se0;
         const ScalarT vr0  = vfe0;
         const ScalarT vhv0 = vr0 / Ka_;
 
@@ -344,15 +361,14 @@ namespace GridKit
           vll0 = vuel0 + inverseRamp(gate_margin0);
         }
 
-        const ScalarT vf0   = ScalarT{ZERO<RealT>};
         const ScalarT ev0   = vll0;
         const ScalarT xll0  = ev0;
-        const ScalarT vref0 = ev0 + vc0 + vf0 - vs0 - uel_on_ * vuel0;
+        const ScalarT vref0 = ev0 + vc0 - vs0 - uel_on_ * vuel0;
 
         y[EFDP] = efdp0;
         y[VC]   = vc0;
         y[VR]   = vr0;
-        y[VF]   = vf0;
+        y[VF]   = ZERO<RealT>;
         y[XLL]  = xll0;
         y[EV]   = ev0;
         y[VLL]  = vll0;
@@ -560,7 +576,7 @@ namespace GridKit
 
         const ScalarT ec                = std::sqrt(wb[0] * wb[0] + wb[1] * wb[1]);
         const ScalarT ev_target         = vref + vs + uel_on_ * vuel - vc - vf;
-        const ScalarT vfe_target        = (Ke_ + se) * efdp;
+        const ScalarT vfe_target        = Ke_eff_ * efdp + se;
         const ScalarT efdp_rate         = (vr - vfe) / Te_;
         const ScalarT limited_efdp_rate = awmin(efdp, efdp_rate, ZERO<RealT>);
 
@@ -749,7 +765,7 @@ namespace GridKit
        * @brief Resolve the parameter-derived constants and selector masks
        *
        * Raises the transducer, regulator, lead-lag, exciter, and feedback
-       * lags to the well-posedness floor, fits the quadratic saturation
+       * lags to the well-posedness floor, fits the scaled-quadratic saturation
        * curve, and turns the three selectors into multiplicative masks. The
        * masks let the residual select signal routing without
        * parameter-dependent control flow, which keeps its structure fixed for
@@ -814,23 +830,37 @@ namespace GridKit
         // The disabled test matches the verify() predicate exactly.
         const bool saturation_disabled =
             Se1_ == ZERO<RealT> && Se2_ == ZERO<RealT>;
-        const bool saturation_points_are_ordered =
-            (E2_ > E1_ && Se2_ > Se1_)
-            || (E2_ < E1_ && Se2_ < Se1_);
+        const bool sat_ordered = (E2_ > E1_ && Se2_ > Se1_) || (E2_ < E1_ && Se2_ < Se1_);
         const bool saturation_consistent =
             E1_ > ZERO<RealT> && E2_ > ZERO<RealT>
-            && Se1_ > ZERO<RealT> && Se2_ > ZERO<RealT>
-            && saturation_points_are_ordered;
+            && Se1_ >= ZERO<RealT> && Se2_ >= ZERO<RealT>
+            && sat_ordered;
         if (saturation_disabled || !saturation_consistent)
         {
           SA_ = ZERO<RealT>;
           SB_ = ZERO<RealT>;
           return;
         }
 
-        const RealT C = std::sqrt(Se2_ / Se1_);
+        if (Se1_ == ZERO<RealT>)
+        {
+          const RealT dE = E2_ - E1_;
+          SA_            = E1_;
+          SB_            = Se2_ * E2_ / (dE * dE);
+          return;
+        }
+
+        if (Se2_ == ZERO<RealT>)
+        {
+          const RealT dE = E1_ - E2_;
+          SA_            = E2_;
+          SB_            = Se1_ * E1_ / (dE * dE);
+          return;
+        }
+
+        const RealT C = std::sqrt(Se2_ * E2_ / (Se1_ * E1_));
         SA_           = (C * E1_ - E2_) / (C - ONE<RealT>);
-        SB_           = Se1_ / ((E1_ - SA_) * (E1_ - SA_));
+        SB_           = Se1_ * E1_ / ((E1_ - SA_) * (E1_ - SA_));
       }
 
       /**
@@ -840,30 +870,17 @@ namespace GridKit
        * *input*, so the residual reproduces the requested output through the
        * same smooth ramp it evaluates.
        *
-       * For large positive values, the ramp is effectively equal to the input, so the
-       * inverse is effectively the output. In that regime this function returns `ramp_output` directly.
-       * This branching is numerically more robust.
-       *
        * @param[in] ramp_output Strictly positive requested ramp output.
        * @return The input the smooth ramp maps to the requested output.
        *
        * @pre @p ramp_output is finite and strictly positive.
-       * @warning This function contains conditional branching and may be used
-       *          during initialization, but not during residual or Jacobian
-       *          evaluation.
        */
       template <typename scalar_type, typename index_type>
       typename Esdc1a<scalar_type, index_type>::RealT
       Esdc1a<scalar_type, index_type>::inverseRamp(RealT ramp_output) const
       {
-        static constexpr RealT SOFTPLUS_WIDTH = static_cast<RealT>(50.0);
-
-        const RealT scaled_output = Math::MU<RealT> * ramp_output;
-        if (scaled_output > SOFTPLUS_WIDTH)
-        {
-          return ramp_output;
-        }
-        return std::log(std::expm1(scaled_output)) / Math::MU<RealT>;
+        const RealT mu = Math::MU<RealT>;
+        return ramp_output + std::log(-std::expm1(-mu * ramp_output)) / mu;
       }
 
       /**