Support zero IEEET1 transducer time constant

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

@@ -261,15 +261,15 @@ namespace GridKit
       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
-        tag_[3] = true;  // y3 - vfx  - Exciter feedback
-        tag_[4] = false; // y4 - vtr  - Term Volt Err
-        tag_[5] = false; // y5 - vf   - Feedback volt
-        tag_[6] = false; // y6 - ve   - Excit. Cntrl Volt
-        tag_[7] = false; // y7 - efd  - Efd
-        tag_[8] = false; // y8 - ksat - Saturation
+        tag_[0] = (Tr_ != 0.0); // y0 - vts  - Sensed term volt
+        tag_[1] = true;         // y1 - vr   - Voltage reg
+        tag_[2] = true;         // y2 - efdp - Efd pre mult
+        tag_[3] = true;         // y3 - vfx  - Exciter feedback
+        tag_[4] = false;        // y4 - vtr  - Term Volt Err
+        tag_[5] = false;        // y5 - vf   - Feedback volt
+        tag_[6] = false;        // y6 - ve   - Excit. Cntrl Volt
+        tag_[7] = false;        // y7 - efd  - Efd
+        tag_[8] = false;        // y8 - ksat - Saturation
 
         return 0;
       }
@@ -312,20 +312,20 @@ namespace GridKit
         ScalarT omega     = ws[0];
         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
+        // The 'pre-limit' target of Vr.
+        ScalarT func            = -vr + Ka_ * vtr;
+        ScalarT func_normalized = (func / Ta_) / 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_);
 
         // Internal Differential Equations
-        f[0] = -vts_dot + (Ec - vts) / Tr_;
-        f[1] = -vr_dot + vr_ind * func;
-        f[2] = -efdp_dot + (vr - ve - Ke_ * efdp) / Te_;
-        f[3] = -vfx_dot + vf / Tf_;
+        f[0] = -Tr_ * vts_dot + Ec - vts;
+        f[1] = -Ta_ * vr_dot + vr_ind * func;
+        f[2] = -Te_ * efdp_dot + vr - ve - Ke_ * efdp;
+        f[3] = -Tf_ * vfx_dot + vf;
 
         // 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_;
 

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

@@ -76,7 +76,7 @@ Generally, this system has two solutions. The non-extraneous solution is as foll
 
 Symbol    | Units  | Description                       | Note
 ----------|--------|-----------------------------------|-------
-$V_{ts}$  | [p.u.] | Sensed terminal voltage           |
+$V_{ts}$  | [p.u.] | Sensed terminal voltage           | Algebraic when $T_R=0$
 $V_R$     | [p.u.] | Voltage regulator                 |
 $E_{fd}'$ | [p.u.] | Field-current pre-speed multiplier|
 $V_{fx}$  | [p.u.] | Exciter feedback internal state   |
@@ -116,55 +116,41 @@ $V_S$           | [p.u.] | Input from stabilizer controller               |
 
 ### Differential Equations
 
-The IEEET1 differential equations, as derived from the model diagram. Define the pre-limit derivative of $V_R$
+The IEEET1 differential equations, as derived from the model diagram.
 
-```math
-f = \dfrac{1}{T_A}\left[-V_R + K_A V_{tr}\right]
-```
-
-so that $\dot V_R$ is the anti-windup limited derivative.
 ```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 &= -T_R \dot V_{ts} + E_C - V_{ts} \\
+   0 &= -T_A \dot V_R
+      + \text{indicator}
+        \left(
+          V_R,
+          \dfrac{-V_R + K_A V_{tr}}{500T_A};
+          V_R^{\min},
+          V_R^{\max}
+        \right)
+        \left(-V_R + K_A V_{tr}\right) \\
+   0 &= -T_E \dot E_{fd}' + V_R - V_E - K_E E_{fd}' \\
+   0 &= -T_F \dot V_{fx} + V_f
 \end{aligned}
 ```
 
-CommonMath defines the [Anti-Windup](../../../../CommonMath.md#anti-windup-indicator) target and smooth approximation.
+CommonMath defines the smooth [Anti-Windup](../../../../CommonMath.md#anti-windup-indicator) indicator.
 
 ### 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_\text{spdlm})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
@@ -187,8 +173,10 @@ The machine initializes $E_{fd}$ first. IEEET1 reads that value as $E_{fd,0}$, a
 ```
 
 All internal derivatives initialize to zero.
-## Model Outputs
 
-The field voltage, $E_{fd}$, is an internal model variable.
+## Model Outputs
 
-The magnetic saturation coefficient $k_\text{sat}$ is calculated from $E_{fd}$ using the smooth piecewise version above.
+Output  | 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/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": {"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} }
    ]
 }