Support zero IEEET1 transducer time constant

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

@@ -94,10 +94,6 @@ namespace GridKit
         using MonitorT   = Model::VariableMonitor<Ieeet1, Ieeet1Data>;
 
         Ieeet1(BusT* bus);
-        Ieeet1(SignalT*          efd_signal,
-               SignalT*          speed_signal,
-               BusT*             bus,
-               const ModelDataT& data);
         Ieeet1(BusT*             bus,
                const ModelDataT& data);
         ~Ieeet1();
@@ -127,10 +123,10 @@ namespace GridKit
             const ScalarT*, const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*);
 
       private:
+        static constexpr RealT TR_MINIMUM = 1.0e-3;
+
         // Signal pointers
-        SignalT* efd_signal_;
-        SignalT* speed_signal_;
-        BusT*    bus_;
+        BusT* bus_;
 
         // Model Input parameters
         RealT Tr_;      ///< Time constant for voltage sensing
@@ -149,7 +145,7 @@ namespace GridKit
         RealT Ispdlim_; ///< Speed limit flag indicator
 
         // Model Derived parameters
-        // TODO -> Need to be solved for in instantiation!
+        // Saturation coefficients derived from E1, E2, Se1, and Se2.
         RealT SA_{0};
         RealT SB_{0};
 
@@ -158,8 +154,6 @@ namespace GridKit
         ScalarT vref_{0}; // (Setpoint voltage, can be different from terminal voltage)
         ScalarT vUEL_{0};
         ScalarT vOEL_{0};
-        ScalarT vS_{0};
-        ScalarT Ec_{0}; // "Compensated" terminal measurment, currently unused
 
         /// Component signal extension
         ComponentSignals<ScalarT, IdxT, Ieeet1InternalVariables, Ieeet1ExternalVariables> signals_;

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

@@ -13,12 +13,9 @@
 #include <GridKit/Model/PhasorDynamics/Bus/Bus.hpp>
 #include <GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp>
 #include <GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Data.hpp>
-#include <GridKit/Model/PhasorDynamics/SignalNode/SignalNode.hpp>
 #include <GridKit/Model/VariableMonitorImpl.hpp>
 #include <GridKit/Utilities/Logger/Logger.hpp>
 
-#define _USE_MATH_DEFINES
-
 namespace GridKit
 {
   namespace PhasorDynamics
@@ -37,34 +34,6 @@ namespace GridKit
         size_ = 9;
       }
 
-      /**
-       * @brief  Constructor for IEEET1 Exciter
-       *
-       * @param data  Data object to store parameters
-       * @param bus   Signal used for terminal reference vmag
-       * @param speed Signal used for machine relative speed
-       * @param efd   Signal used for E field
-       */
-      template <typename scalar_type, typename index_type>
-      Ieeet1<scalar_type, index_type>::Ieeet1(SignalT*          efd_signal,
-                                              SignalT*          speed_signal,
-                                              BusT*             bus,
-                                              const ModelDataT& data)
-        : efd_signal_(efd_signal),
-          speed_signal_(speed_signal),
-          bus_(bus),
-          monitor_(std::make_unique<MonitorT>(data))
-      {
-
-        // Parse data struct into model
-        this->initModelParams(data);
-
-        initializeMonitor();
-
-        // 9 Internal Variables
-        size_ = 9;
-      }
-
       /**
        * @brief  Constructor for IEEET1 Exciter
        *
@@ -186,7 +155,7 @@ namespace GridKit
        *
        * Inputs:
        *   - EFD assigned by the generator (read from y_[7]).
-       *   - Bus voltage, used to form the sensed terminal voltage Ec.
+       *   - 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.
@@ -261,7 +230,6 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       int Ieeet1<scalar_type, index_type>::tagDifferentiable()
       {
-
         tag_[0] = true;  // y0 - vts  - Sensed term volt
         tag_[1] = true;  // y1 - vr   - Voltage reg
         tag_[2] = true;  // y2 - efdp - Efd pre mult
@@ -306,7 +274,7 @@ namespace GridKit
           const ScalarT* ws,
           ScalarT*       f)
       {
-        // Read E comp (terminal voltage, unless compensation impedance)
+        // Read bus voltage components
         ScalarT vreal = wb[0];
         ScalarT vimag = wb[1];
         ScalarT Ec    = std::sqrt(vreal * vreal + vimag * vimag);
@@ -333,22 +301,19 @@ namespace GridKit
         ScalarT vs_signal = ws[1];
 
         // The 'pre-limit' derivative of Vr.
-        ScalarT func            = (-vr + Ka_ * vtr) / Ta_;
-        ScalarT func_normalized = func / static_cast<RealT>(500.0); // TODO This is arbitrary, need more general conditioning method that is fast
-        ScalarT vr_ind          = Math::indicator(vr, func_normalized, Vrmin_, Vrmax_);
+        ScalarT func = (-vr + Ka_ * vtr) / Ta_;
 
         // Internal Differential Equations
         f[0] = -vts_dot + (Ec - vts) / Tr_;
-        f[1] = -vr_dot + vr_ind * func;
+        f[1] = -vr_dot + Math::antiwindup(vr, func, Vrmin_, Vrmax_);
         f[2] = -efdp_dot + (vr - ve - Ke_ * efdp) / Te_;
         f[3] = -vfx_dot + vf / Tf_;
 
         // Internal Algebraic Equations
         f[4] = -vts + vref_ + vUEL_ + vOEL_ + vs_signal - vtr - vf;
-        f[5] = -vf + (efdp * Kf_) / Tf_ - vfx;
+        f[5] = -Tf_ * (vf + vfx) + Kf_ * efdp;
         f[6] = -ve + ksat * efdp;
         f[7] = -efd + efdp + omega * efdp * Ispdlim_;
-
         f[8] = -ksat + SB_ * Math::qramp(efdp - SA_);
 
         return 0;
@@ -361,15 +326,7 @@ namespace GridKit
       template <typename scalar_type, typename index_type>
       int Ieeet1<scalar_type, index_type>::evaluateResidual()
       {
-        // Set Input Variables
-        // Meta PR Note: This seems to be very slow,
-        // but I see how read/write ownership may require this
-        //
-        // I believe implementing the equivalent to signal->read()
-        // at the system level would address this, by routing
-        // external signals into a generic inputs_ vector
-        // at the same time as the internal state values y_
-        // are recieved from IDA.
+        // Set input variables.
         if (signals_.template isAttached<Ieeet1ExternalVariables::OMEGA>())
         {
           ws_[0]         = signals_.template readExternalVariable<Ieeet1ExternalVariables::OMEGA>();
@@ -402,10 +359,15 @@ namespace GridKit
       void Ieeet1<scalar_type, index_type>::initModelParams(const ModelDataT& data)
       {
 
+        Tr_ = TR_MINIMUM;
         if (data.parameters.contains(ModelDataT::Parameters::Tr))
         {
           Tr_ = std::get<RealT>(data.parameters.at(ModelDataT::Parameters::Tr));
         }
+        if (Tr_ < TR_MINIMUM)
+        {
+          Tr_ = TR_MINIMUM;
+        }
         if (data.parameters.contains(ModelDataT::Parameters::Ka))
         {
           Ka_ = std::get<RealT>(data.parameters.at(ModelDataT::Parameters::Ka));
@@ -478,7 +440,7 @@ namespace GridKit
       {
         using Variable = ModelDataT::MonitorableVariables;
         monitor_->set(Variable::efd, [this]
-                      { return efd_signal_->read(); });
+                      { return y_[7]; });
         monitor_->set(Variable::ksat, [this]
                       { return SB_ * Math::qramp(y_[2] - SA_); });
       }

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

@@ -5,9 +5,7 @@
 Standard model of the IEEET1 Exciter.
 
 Notes:
-- $E_C$ should be an external signal from the generator or machine, not computed through an exciter-owned `Bus` reference.
-- The current implementation uses its `Bus` reference as a proxy for $E_C$.
-- This direct coupling affects numerical conditioning; production models typically use a decoupling reactance for the exciter-current path that forms $E_C$.
+- The voltage-sensing input currently uses the positive bus-voltage magnitude $\sqrt{V_r^2 + V_i^2}$; a separate compensation-impedance path is not modeled here.
 
 ![](../../../../../docs/Figures/PhasorDynamics_IEEET1_Diagram.png)
 
@@ -31,7 +29,7 @@ $E_1$       | [p.u.] | Saturation Parameter                 | 2.8     |
 $E_2$       | [p.u.] | Saturation Parameter                 | 3.73    |
 $S_1$       | [p.u.] | Saturation Parameter                 | 0.04    |
 $S_2$       | [p.u.] | Saturation Parameter                 | 0.33    |
-$I_\text{spdlm}$ | [binary] | Speed limit flag indicator       | 0       |
+$I_{\mathrm{spdlim}}$ | [binary] | Speed limit flag indicator       | 0       |
 
 ### Model Derived Parameters
 
@@ -71,12 +69,12 @@ Generally, this system has two solutions. The non-extraneous solution is as foll
 
 #### Differential
 
-Symbol    | Units  | Description                       | Note
-----------|--------|-----------------------------------|-------
-$V_{ts}$  | [p.u.] | Sensed terminal voltage           |
-$V_R$     | [p.u.] | Voltage regulator                 |
-$E_{fd}'$ | [p.u.] | Field-current pre-speed multiplier|
-$V_{fx}$  | [p.u.] | Exciter feedback internal state   |
+Symbol    | Units  | Description                        | Note
+----------|--------|------------------------------------|-------
+$V_{ts}$  | [p.u.] | Sensed terminal voltage            |
+$V_R$     | [p.u.] | Voltage regulator                  |
+$E_{fd}'$ | [p.u.] | Field-current pre-speed multiplier |
+$V_{fx}$  | [p.u.] | Exciter feedback internal state    |
 
 
 #### Algebraic
@@ -92,100 +90,85 @@ $k_\text{sat}$  | [p.u.] | Saturation variable               |
 
 ### External Variables
 
-#### Differential
-Symbol          | Units  | Description                       | Note
-----------------|--------|-----------------------------------|-------
-$\omega$  | [p.u.] | Machine Speed Deviation                   | Read from a Machine Model
-
-#### Algebraic
-
-
 Symbol          | Units  | Description                       | Note
 ----------------|--------|-----------------------------------|-------
-$E_C$           | [p.u.] | Compensated machine terminal voltage magnitude |
+$V_r$           | [p.u.] | Real bus voltage component                     |
+$V_i$           | [p.u.] | Imaginary bus voltage component                |
 $V_\text{ref}$  | [p.u.] | Reference terminal voltage                     |
-$V_{UEL}$       | [p.u.] | Input from under excitation limiter            |
-$V_{OEL}$       | [p.u.] | Input from over excitation limiter             |
-$V_S$           | [p.u.] | Input from stabilizer controller               |
+$V_{UEL}$       | [p.u.] | Input from under excitation limiter            | Constant zero until modeled
+$V_{OEL}$       | [p.u.] | Input from over excitation limiter             | Constant zero until modeled
+$V_S$           | [p.u.] | Input from stabilizer controller               | Optional, defaults to zero
+$\omega$        | [p.u.] | Machine speed deviation                        | Optional, defaults to zero
 
 
 ## Model Equations
 
 ### Differential Equations
 
-The IEEET1 differential equations, as derived from the model diagram. Define the pre-limit derivative of $V_R$
+For readability, define the pre-limit derivative of $V_R$ and voltage-sensing input:
 
 ```math
-f = \dfrac{1}{T_A}\left[-V_R + K_A V_{tr}\right]
+\begin{aligned}
+f_R &:= \dfrac{1}{T_A}\left(-V_R + K_A V_{tr}\right) \\
+E_C &:= \sqrt{V_r^2 + V_i^2}
+\end{aligned}
 ```
 
-so that $\dot V_R$ is the anti-windup limited derivative.
+The IEEET1 differential equations, as derived from the model diagram, are:
+
 ```math
 \begin{aligned}
-   \dot V_{ts}   &= \dfrac{1}{T_R}(E_C-V_{ts}) \\
-   \dot V_R      &= \text{antiwindup}
-      \left(V_R, f;\ V_R^{\min}, V_R^{\max}\right) \\
-   \dot E_{fd}'  &= \dfrac{1}{T_E}(V_R-V_E-K_E E_{fd}') \\
-   \dot V_{fx}   &= \dfrac{1}{T_F}V_f \\
+   0 &= -\dot V_{ts} + \dfrac{1}{T_R}\left(E_C - V_{ts}\right) \\
+   0 &= -\dot V_R
+      + \text{antiwindup}
+        \left(V_R, f_R; V_R^{\min}, V_R^{\max}\right) \\
+   0 &= -\dot E_{fd}' + \dfrac{1}{T_E}\left(V_R - V_E - K_E E_{fd}'\right) \\
+   0 &= -\dot V_{fx} + \dfrac{1}{T_F}\left(V_f\right)
 \end{aligned}
 ```
 
-CommonMath defines the [Anti-Windup](../../../../CommonMath.md#antiwindup) target and smooth approximation.
+CommonMath defines the smooth [Anti-Windup](../../../../CommonMath.md#anti-windup-indicator) target and approximation.
 
 ### Algebraic Equations
 
 The algebraic equations of the exciter.
 ```math
 \begin{aligned}
-   V_{tr} &= V_\text{ref} +V_{UEL} + V_{OEL} + V_S - V_f- V_{ts}\\
-    V_f &= \dfrac{E_{fd}' K_F}{T_F} - V_{fx}\\
-    V_E &= k_\text{sat}\cdot E_{fd}' \\
-    E_{fd}&= \begin{cases}
-        (1+\omega)E_{fd}'           &  \text{if } I_\text{spdlm}=1\\
-         E_{fd}' &  \text{else} \\
-   \end{cases}\\
-    k_\text{sat}&= \begin{cases}
-        S_B(E_{fd}' -S_A)^2        &  \text{if } E_{fd}' >S_A\\
-        0  &  \text{else } \\
-   \end{cases} \\
+   0 &= -V_{ts} + V_\text{ref} + V_{UEL} + V_{OEL} + V_S - V_{tr} - V_f \\
+   0 &= -T_F(V_f + V_{fx}) + K_F E_{fd}' \\
+   0 &= -V_E + k_\text{sat} E_{fd}' \\
+   0 &= -E_{fd} + (1 + \omega I_{\mathrm{spdlim}})E_{fd}' \\
+   0 &= -k_\text{sat} + S_B\, q(E_{fd}' - S_A)
 \end{aligned}
 ```
-#### Smooth Piecewise Approximation (Algebraic)
 
-For the algebraic piecewise functions (non-flags), this implementation is straightforward when the approximation above is used.
 Here $q$ is GridKit's [Quadratic Ramp](../../../../CommonMath.md#primitives).
-```math
-\begin{aligned}
-    E_{fd}
-    &=(1 + \omega I_\text{spdlm})E_{fd}' \\
-    k_\text{sat}
-    &=S_B\, q(E_{fd}' -S_A)
-\end{aligned}
-```
 
 
 ## Initialization
 
-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$. There is no compensation impedance, so $E_C$ is taken as the terminal 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 auxiliary inputs.
+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.
 
 ```math
 \begin{aligned}
-   E_C      &= \sqrt{V_r^2 + V_i^2} \\
-   E_{fd}'  &= \dfrac{E_{fd,0}}{1 + I_\text{spdlm}\,\omega} \\
+   E_{C,0}  &:= \sqrt{V_r^2 + V_i^2} \\
+   E_{fd}'  &= \dfrac{E_{fd,0}}{1 + I_{\mathrm{spdlim}}\,\omega} \\
    k_\text{sat}  &= S_B\, q(E_{fd}' - S_A) \\
    V_E      &= k_\text{sat}\, E_{fd}' \\
    V_R      &= K_E\, E_{fd}' + V_E \\
    V_{tr}   &= \dfrac{V_R}{K_A} \\
    V_{fx}   &= \dfrac{K_F}{T_F}\, E_{fd}' \\
-   V_{ts}   &= E_C \\
+   V_{ts}   &= E_{C,0} \\
    V_f      &= 0 \\
-   V_\text{ref}  &= E_C + V_{tr} - V_{UEL} - V_{OEL} - V_S
+   V_\text{ref}  &= E_{C,0} + V_{tr} - V_{UEL} - V_{OEL} - V_S
 \end{aligned}
 ```
 
 All internal derivatives initialize to zero.
-## Model Outputs
 
-The field voltage, $E_{fd}$, is an internal model variable.
+## Monitorable Variables
 
-The magnetic saturation coefficient $k_\text{sat}$ is calculated from $E_{fd}$ using the smooth piecewise version above.
+Variable | Units  | Description                       | Note
+---------|--------|-----------------------------------|------
+`efd`    | [p.u.] | Field winding voltage             |
+`ksat`   | [p.u.] | Magnetic saturation coefficient   | $S_B\,q(E_{fd}'-S_A)$

GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp

@@ -217,12 +217,10 @@ namespace GridKit
         ScalarT Ec = std::sqrt(wb[0] * wb[0] + wb[1] * wb[1]);
         ScalarT vs = ws[0];
 
-        ScalarT func            = (-efd + (K_ / Tb_) * (-vr + Ta_ * vtr)) / Te_;
-        ScalarT func_normalized = func / static_cast<RealT>(100.0); // TODO arbitrary conditioning; see IEEET1
-        ScalarT efd_ind         = Math::indicator(efd, func_normalized, Efdmin_, Efdmax_);
+        ScalarT func = (-efd + (K_ / Tb_) * (-vr + Ta_ * vtr)) / Te_;
 
         f[0] = -vr_dot + (-vr + Ta_ * vtr) / Tb_ - vtr;
-        f[1] = -efd_dot + efd_ind * func;
+        f[1] = -efd_dot + Math::antiwindup(efd, func, Efdmin_, Efdmax_);
         f[2] = -vtr - Ec + vref_ + vs + vOEL_ + vUEL_;
 
         return 0;

GridKit/Model/PhasorDynamics/INPUT_FORMAT.md

@@ -185,7 +185,7 @@ side for off-nominal transformer branches.
        { "class": "Branch", "ports": {"bus1":1, "bus2":2}, "id": "BR1", "params": {"R":0.0, "X":0.1, "G":0.0, "B":0.0} },
        { "class": "Genrou", "ports": {"bus":1, "speed": 1, "pmech":2, "efd":3}, "id": "DV1", "params": {"p0":1.0, "q0":0.05013, "H":3.0, "D":0.0, "Ra":0.0, "Tdop":7.0, "Tdopp":0.04, "Tqopp":0.05, "Tqop":0.75, "Xd":2.1, "Xdp":0.2, "Xdpp":0.18, "Xq":0.5, "Xqp": 0.0, "Xqpp":0.18, "Xl":0.15, "S10":0.0, "S12":0.0}, "mon": ["delta", "omega"] },
        { "class": "Tgov1", "ports": {"bus":1, "speed": 1, "pmech":2}, "id": "DV2", "params": {"Trate":100.0, "R":0.05, "T1":0.5,"T2":2.5, "T3":7.5, "Pvmax":0, "Pvmin":1, "Dt":0}},
-       { "class": "Ieeet1", "ports": {"bus":1, "speed": 1, "efd":3}, "id": "DV3", "params": {"Tr":0.001, "Ka":50.0, "Ta":0.04, "Ke":-0.06, "Te":0.6, "Kf":0.09, "Tf":1.46, "Vrmin":-1, "Vrmax":1, "E1":2.8, "E2":3.373, "Se1":0.04, "Se2":0.33, "Ispdlim":0}},
+       { "class": "Ieeet1", "ports": {"bus":1, "speed": 1, "efd":3}, "id": "DV3", "params": {"Tr":0.0, "Ka":50.0, "Ta":0.04, "Ke":-0.06, "Te":0.6, "Kf":0.09, "Tf":1.46, "Vrmin":-1, "Vrmax":1, "E1":2.8, "E2":3.373, "Se1":0.04, "Se2":0.33, "Ispdlim":0}},
        { "class": "BusFault", "ports": {"bus":1}, "id": "EVT1", "params": {"state0": false, "R":0.0, "X":1e-3} }
    ]
 }