From 72abbdd8cd9ff44ffcd2f0348e99e06c7660be3e Mon Sep 17 00:00:00 2001 From: lukelowry Date: Thu, 13 Aug 2026 15:13:13 -0500 Subject: [PATCH 01/10] minimal corrections --- .../PhasorDynamics/Converter/REGCA/README.md | 25 +-- .../Converter/REGCA/RegcaImpl.hpp | 7 +- .../PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp | 2 +- .../Exciter/ESDC1A/Esdc1aData.hpp | 2 +- .../Exciter/ESDC1A/Esdc1aImpl.hpp | 81 +++++---- .../PhasorDynamics/Exciter/ESDC1A/README.md | 55 ++++-- .../PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp | 2 +- .../Exciter/IEEET1/Ieeet1Impl.hpp | 54 ++++-- .../PhasorDynamics/Exciter/IEEET1/README.md | 46 +++-- .../PhasorDynamics/Governor/Tgov1/README.md | 123 +++++++------ .../PhasorDynamics/Governor/Tgov1/Tgov1.hpp | 32 ++-- .../Governor/Tgov1/Tgov1Data.hpp | 21 +-- .../Governor/Tgov1/Tgov1Impl.hpp | 162 ++++++++++++------ .../SynchronousMachine/GENROU/Genrou.hpp | 12 +- .../SynchronousMachine/GENROU/GenrouImpl.hpp | 138 ++++++--------- .../SynchronousMachine/GENROU/README.md | 67 +++----- .../SynchronousMachine/GENSAL/Gensal.hpp | 40 ++--- .../SynchronousMachine/GENSAL/GensalData.hpp | 70 ++++---- .../SynchronousMachine/GENSAL/GensalImpl.hpp | 44 +++-- .../SynchronousMachine/GENSAL/README.md | 13 +- .../Model/PhasorDynamics/SystemModelImpl.hpp | 7 + .../PhasorDynamics/ConverterRegcaTests.hpp | 2 - .../PhasorDynamics/ExciterEsdc1aTests.hpp | 125 +++++++------- .../PhasorDynamics/ExciterIeeet1Tests.hpp | 8 +- .../UnitTests/PhasorDynamics/GensalTests.hpp | 2 +- .../PhasorDynamics/GovernorTgov1Tests.hpp | 6 +- .../SystemSingleComponentTests.hpp | 3 + tests/UnitTests/Utilities/CaseFormatTests.hpp | 7 +- 28 files changed, 640 insertions(+), 516 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md b/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md index f3c59c7f7..269d122cb 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md +++ b/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} diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp index d5a8ae03c..456e5ef75 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp @@ -89,11 +89,11 @@ namespace GridKit use_rqmax_ = ZERO; use_rqmin_ = ZERO; - if (q0_ > ZERO) + if (q0_ > ZERO && Rqmax_ > ZERO) { use_rqmax_ = ONE; } - else if (q0_ < ZERO) + else if (q0_ < ZERO && Rqmin_ < ZERO) { use_rqmin_ = ONE; } @@ -436,8 +436,7 @@ namespace GridKit } check(mva_base_ > ZERO, "mva must be positive"); - check(Rpmax_ > ZERO, "Rpmax must be positive"); - check(Rqmin_ < ZERO && ZERO < Rqmax_, "Rqmin < 0 < Rqmax is required"); + check(Rpmax_ >= ZERO, "Rpmax must be non-negative"); check(IL1_ >= ZERO, "IL1 must be non-negative"); check(KL_ > ZERO, "LVPL release slope must be positive"); check(ZERO <= VL0_ && VL0_ < VL1_, "VL0/VL1 must satisfy 0 <= VL0 < VL1"); diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp index dd42f55ac..b7f791886 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp +++ b/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$E_{\mathrm{fd}}'S_E(E_{\mathrm{fd}}')\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 diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aData.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aData.hpp index 0dfc3bc42..3de51f210 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aData.hpp +++ b/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$E_{\mathrm{fd}}'S_E(E_{\mathrm{fd}}')\f$ Scaled-quadratic saturation contribution [p.u.] vfe ///< \f$V_{\mathrm{FE}}\f$ Exciter feedback drive [p.u.] }; diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp index 11906d3d9..2a0e43961 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp @@ -162,7 +162,6 @@ namespace GridKit } check(Ka_ > ZERO, "Ka must be positive"); - check(Tc_ >= ZERO, "Tc must be non-negative"); check(Vrmin_ <= Vrmax_, "Vrmin must be less than or equal to Vrmax"); check(UEL_ >= static_cast(0) && UEL_ <= static_cast(3), "UEL must be 0, 1, 2, or 3"); @@ -176,14 +175,11 @@ namespace GridKit { check(E1_ > ZERO, "E1 must be positive when saturation is enabled"); check(E2_ > ZERO, "E2 must be positive when saturation is enabled"); - check(Se1_ > ZERO, "Se1 must be positive when saturation is enabled"); - check(Se2_ > ZERO, "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, "Se1 must be non-negative when saturation is enabled"); + check(Se2_ >= ZERO, "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()) @@ -314,8 +310,15 @@ 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_); + + RealT ke0 = Ke_; + if (ke0 == ZERO) + { + ke0 = (Vrmax_ / 10.0 - static_cast(se0)) / static_cast(efdp0); + } + + const ScalarT vfe0 = ke0 * efdp0 + se0; const ScalarT vr0 = vfe0; const ScalarT vhv0 = vr0 / Ka_; @@ -344,15 +347,16 @@ namespace GridKit vll0 = vuel0 + inverseRamp(gate_margin0); } - const ScalarT vf0 = ScalarT{ZERO}; 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; + + Ke_ = ke0; y[EFDP] = efdp0; y[VC] = vc0; y[VR] = vr0; - y[VF] = vf0; + y[VF] = ZERO; y[XLL] = xll0; y[EV] = ev0; y[VLL] = vll0; @@ -560,7 +564,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_ * efdp + se; const ScalarT efdp_rate = (vr - vfe) / Te_; const ScalarT limited_efdp_rate = awmin(efdp, efdp_rate, ZERO); @@ -749,7 +753,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,13 +818,11 @@ namespace GridKit // The disabled test matches the verify() predicate exactly. const bool saturation_disabled = Se1_ == ZERO && Se2_ == ZERO; - 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 && E2_ > ZERO - && Se1_ > ZERO && Se2_ > ZERO - && saturation_points_are_ordered; + && Se1_ >= ZERO && Se2_ >= ZERO + && sat_ordered; if (saturation_disabled || !saturation_consistent) { SA_ = ZERO; @@ -828,9 +830,25 @@ namespace GridKit return; } - const RealT C = std::sqrt(Se2_ / Se1_); + if (Se1_ == ZERO) + { + const RealT dE = E2_ - E1_; + SA_ = E1_; + SB_ = Se2_ * E2_ / (dE * dE); + return; + } + + if (Se2_ == ZERO) + { + 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); - SB_ = Se1_ / ((E1_ - SA_) * (E1_ - SA_)); + SB_ = Se1_ * E1_ / ((E1_ - SA_) * (E1_ - SA_)); } /** @@ -840,30 +858,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 Esdc1a::RealT Esdc1a::inverseRamp(RealT ramp_output) const { - static constexpr RealT SOFTPLUS_WIDTH = static_cast(50.0); - - const RealT scaled_output = Math::MU * ramp_output; - if (scaled_output > SOFTPLUS_WIDTH) - { - return ramp_output; - } - return std::log(std::expm1(scaled_output)) / Math::MU; + const RealT mu = Math::MU; + return ramp_output + std::log(-std::expm1(-mu * ramp_output)) / mu; } /** diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md index aed2f0d58..9840a4646 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md @@ -54,7 +54,7 @@ Invalid ESDC1A parameter sets are rejected by the following checks: \begin{aligned} K_A &> 0 \\ - T_R, T_A, T_B, T_C, T_E, T_{F1} + T_R, T_A, T_B, T_E, T_{F1} &\ge 0 \\ V_R^{\min} &\le V_R^{\max} \\ @@ -71,11 +71,12 @@ The saturation points are either disabled together, S_E(E_1) = S_E(E_2) = 0, ``` -or define a valid two-point quadratic fit: +or define a valid two-point scaled-quadratic fit: ```math \begin{aligned} - E_1, E_2, S_E(E_1), S_E(E_2) &> 0 \\ + E_1, E_2 &> 0 \\ + S_E(E_1), S_E(E_2) &\ge 0 \\ \left(E_2-E_1\right) \left[S_E(E_2)-S_E(E_1)\right] &> 0 \end{aligned} @@ -99,14 +100,30 @@ raised to that floor in place, so every equation below uses the raised value: \end{cases} \end{aligned} ``` - When saturation is disabled, $S_A = 0$ and $S_B = 0$. Otherwise, +```math +E S_E(E) = S_B q(E-S_A). +``` + +When one saturation value is zero, + +```math +\begin{aligned} + S_E(E_1)=0 &: \quad S_A=E_1,\qquad + S_B=\dfrac{E_2S_E(E_2)}{(E_2-E_1)^2} \\ + S_E(E_2)=0 &: \quad S_A=E_2,\qquad + S_B=\dfrac{E_1S_E(E_1)}{(E_1-E_2)^2}. +\end{aligned} +``` + +and when both saturation values are positive, + ```math \begin{aligned} - C &= \sqrt{\dfrac{S_E(E_2)}{S_E(E_1)}} \\ + C &= \sqrt{\dfrac{E_2S_E(E_2)}{E_1S_E(E_1)}} \\ S_A &= \dfrac{C E_1 - E_2}{C - 1} \\ - S_B &= \dfrac{S_E(E_1)}{(E_1 - S_A)^2} + S_B &= \dfrac{E_1S_E(E_1)}{(E_1 - S_A)^2} \end{aligned} ``` @@ -149,7 +166,7 @@ Symbol | Units | Description $e_V$ | [p.u.] | Voltage-error summing output | $V_{\mathrm{LL}}$ | [p.u.] | Input lead-lag output | $V_{\mathrm{HV}}$ | [p.u.] | High-value gate output | -$S_E$ | [p.u.] | Exciter saturation coefficient | Evaluated at $E_{\mathrm{fd}}'$ +$E_{\mathrm{fd}}'S_E(E_{\mathrm{fd}}')$ | [p.u.] | Scaled-quadratic saturation contribution | $V_{\mathrm{FE}}$ | [p.u.] | Exciter feedback drive | $E_{\mathrm{fd}}$ | [p.u.] | Field-voltage output | Published through `efd` @@ -246,11 +263,12 @@ of [Appendix A](#appendix-a-awmin). & s_{\mathrm{UEL}} = 1 \end{cases} \\ 0 &= - -S_E + -E_{\mathrm{fd}}'S_E(E_{\mathrm{fd}}') + S_B q\left(E_{\mathrm{fd}}' - S_A\right) \\ 0 &= -V_{\mathrm{FE}} - + \left(K_E + S_E\right)E_{\mathrm{fd}}' \\ + + K_E E_{\mathrm{fd}}' + + E_{\mathrm{fd}}'S_E(E_{\mathrm{fd}}') \\ 0 &= -E_{\mathrm{fd}} + \left(1 + s_{\mathrm{spd}}\omega\right)E_{\mathrm{fd}}' @@ -297,10 +315,17 @@ routed through the gate: E_{\mathrm{fd}}' &\leftarrow \dfrac{E_{\mathrm{fd}}}{1 + s_{\mathrm{spd}}\omega} \\ - S_E + E_{\mathrm{fd}}'S_E(E_{\mathrm{fd}}') &\leftarrow S_B q\left(E_{\mathrm{fd}}' - S_A\right) \\ + K_E + &\leftarrow + \begin{cases} + \dfrac{V_R^{\max}/10-E_{\mathrm{fd}}'S_E(E_{\mathrm{fd}}')}{E_{\mathrm{fd}}'} & K_E=0 \\ + K_E & K_E\ne 0 + \end{cases} \\ V_{\mathrm{FE}} - &\leftarrow \left(K_E + S_E\right)E_{\mathrm{fd}}' \\ + &\leftarrow K_E E_{\mathrm{fd}}' + + E_{\mathrm{fd}}'S_E(E_{\mathrm{fd}}') \\ V_R &\leftarrow V_{\mathrm{FE}} \\ V_{\mathrm{HV}} @@ -327,9 +352,9 @@ routed through the gate: Initialization rejects a non-finite or zero bus-voltage magnitude, a non-finite field-voltage seed, non-finite Known signal inputs, a nonpositive speed multiplier $1 + s_{\mathrm{spd}}\omega$, $E_{\mathrm{fd}}'<0$ while -$s_{\mathrm{lim}}=1$, $V_R$ outside -$[V_R^{\min},V_R^{\max}]$, and high-value-gate active starts with -$s_{\mathrm{UEL}} = 0$ and +$s_{\mathrm{lim}}=1$, initial $V_R$ outside $[V_R^{\min},V_R^{\max}]$, +and high-value-gate active +starts with $s_{\mathrm{UEL}} = 0$ and $V_{\mathrm{HV}}\le V_{\mathrm{UEL}}$. Every check resolves before any storage is written, so a rejected @@ -361,7 +386,7 @@ Output | Units | Description | Note `vc` | [p.u.] | Filtered terminal-voltage magnitude | $V_C$ `vr` | [p.u.] | Voltage-regulator output | $V_R$ `vf` | [p.u.] | Stabilizing feedback state | $V_F$ -`se` | [p.u.] | Exciter saturation coefficient | $S_E$ +`se` | [p.u.] | Scaled-quadratic saturation contribution | $E_{\mathrm{fd}}'S_E(E_{\mathrm{fd}}')$ `vfe` | [p.u.] | Exciter feedback drive | $V_{\mathrm{FE}}$ ## Testing diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp index 36d03a16f..de663f014 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp @@ -50,7 +50,7 @@ namespace GridKit VF, ///< Feedback voltage VE, ///< Exciter control voltage EFD, ///< Efd - KSAT, ///< Saturation + KSAT, ///< \f$E_{\mathrm{fd}}'S(E_{\mathrm{fd}}')\f$ Scaled-quadratic saturation contribution MAXIMUM, }; diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp index 58d016ae6..a7eb0eb4b 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp @@ -8,6 +8,7 @@ */ #include +#include #include #include @@ -148,10 +149,11 @@ namespace GridKit { check(E1_ > ZERO, "E1 must be positive when saturation is enabled"); check(E2_ > ZERO, "E2 must be positive when saturation is enabled"); - check(Se1_ > ZERO, "Se1 must be positive when saturation is enabled"); - check(Se2_ > ZERO, "Se2 must be positive when saturation is enabled"); - check(E1_ != E2_, "E1 and E2 must differ when saturation is enabled"); - check(Se1_ != Se2_, "Se1 and Se2 must differ when saturation is enabled"); + check(Se1_ >= ZERO, "Se1 must be non-negative when saturation is enabled"); + check(Se2_ >= ZERO, "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 isAttached()) @@ -191,6 +193,11 @@ namespace GridKit template int Ieeet1::initialize() { + if (verify() != 0) + { + Log::error() << "Ieeet1: cannot initialize with invalid configuration\n"; + return 1; + } // External Variables ScalarT efd0{0}; @@ -226,9 +233,14 @@ namespace GridKit ScalarT efdp = efd0 / (ONE + omega * Ispdlim_); ScalarT ksat = SB_ * Math::qramp(efdp - SA_); - ScalarT ve = ksat * efdp; - ScalarT vr = Ke_ * efdp + ve; - ScalarT vtr = vr / Ka_; + if (Ke_ == ZERO) + { + Ke_ = (Vrmax_ / 10.0 - static_cast(ksat)) / static_cast(efdp); + } + + ScalarT ve = ksat; + ScalarT vr = Ke_ * efdp + ve; + ScalarT vtr = vr / Ka_; ScalarT vf{0}; ScalarT vfx = (Kf_ / Tf_) * efdp; @@ -345,7 +357,7 @@ namespace GridKit // Internal Algebraic Equations f[4] = -vts + vref_ + vUEL_ + vOEL_ + vs_signal - vtr - vf; f[5] = -Tf_ * (vf + vfx) + Kf_ * efdp; - f[6] = -ve + ksat * efdp; + f[6] = -ve + ksat; f[7] = -efd + efdp + omega * efdp * Ispdlim_; f[8] = -ksat + SB_ * Math::qramp(efdp - SA_); @@ -471,17 +483,35 @@ namespace GridKit return; } - if (E1_ <= ZERO || E2_ <= ZERO || E1_ == E2_ - || Se1_ <= ZERO || Se2_ <= ZERO || Se1_ == Se2_) + const bool sat_ordered = (E2_ > E1_ && Se2_ > Se1_) || (E2_ < E1_ && Se2_ < Se1_); + if (E1_ <= ZERO || E2_ <= ZERO + || Se1_ < ZERO || Se2_ < ZERO + || !sat_ordered) + { + return; + } + + if (Se1_ == ZERO) + { + const RealT dE = E2_ - E1_; + SA_ = E1_; + SB_ = Se2_ * E2_ / (dE * dE); + return; + } + + if (Se2_ == ZERO) { + const RealT dE = E1_ - E2_; + SA_ = E2_; + SB_ = Se1_ * E1_ / (dE * dE); return; } - const RealT C = std::sqrt(Se2_ / Se1_); + const RealT C = std::sqrt(Se2_ * E2_ / (Se1_ * E1_)); // Solution 1 (Aligned with PW) SA_ = (C * E1_ - E2_) / (C - ONE); - SB_ = Se1_ / ((E1_ - SA_) * (E1_ - SA_)); + SB_ = Se1_ * E1_ / ((E1_ - SA_) * (E1_ - SA_)); } template diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md index 4985920df..4f689e567 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md @@ -49,9 +49,8 @@ Invalid IEEET1 parameter sets are rejected by the following checks. Let $\epsilo &=(0,0) \quad\text{or}\quad \begin{gathered} - E_1, E_2, S_1, S_2 > 0 \\ - E_1 \ne E_2 \\ - S_1 \ne S_2 + E_1, E_2 > 0,\quad S_1, S_2 \ge 0 \\ + (E_2-E_1)(S_2-S_1) > 0 \end{gathered} \end{aligned} ``` @@ -59,24 +58,36 @@ Invalid IEEET1 parameter sets are rejected by the following checks. Let $\epsilo ### Model Derived Parameters When saturation is disabled, $S_A=0$ and $S_B=0$. Otherwise, -the parameters are chosen so that the following quadratic model represents the -expected saturation near the operating region: +the parameters are chosen for the scaled-quadratic saturation model: ```math \begin{aligned} - S_1 &= S_B(E_1-S_A)^2 \\ - S_2 &= S_B(E_2-S_A)^2 \\ + E S(E) &= S_B q(E-S_A) \\ + E_1S_1 &= S_B(E_1-S_A)^2 \\ + E_2S_2 &= S_B(E_2-S_A)^2 \\ \end{aligned} ``` -Generally, this system has two solutions. The non-extraneous solution is as follows. +When exactly one saturation value is zero, the normal curve fit uses the +corresponding voltage as the quadratic knee: + +```math +\begin{aligned} + S_1=0 &: \quad S_A=E_1,\qquad + S_B=\dfrac{E_2S_2}{(E_2-E_1)^2} \\ + S_2=0 &: \quad S_A=E_2,\qquad + S_B=\dfrac{E_1S_1}{(E_1-E_2)^2}. +\end{aligned} +``` + +When both saturation values are positive, the non-extraneous solution is: ```math \begin{aligned} C &= \sqrt{ \dfrac - {S_2} - {S_1} + {E_2S_2} + {E_1S_1} } \\ S_A &= @@ -86,7 +97,7 @@ Generally, this system has two solutions. The non-extraneous solution is as foll \\ S_B &= \dfrac - {S_1} + {E_1S_1} {(E_1-S_A)^2} \end{aligned} ``` @@ -114,7 +125,7 @@ $V_{tr}$ | [p.u.] | Terminal Voltage Error | $V_f$ | [p.u.] | Feedback Voltage | $V_E$ | [p.u.] | Excitation control voltage | $E_{fd}$ | [p.u.] | Field winding voltage | -$k_\text{sat}$ | [p.u.] | Saturation variable | +$k_\text{sat}$ | [p.u.] | Scaled-quadratic saturation contribution | $E_{fd}'S(E_{fd}')$ ### External Variables @@ -164,7 +175,7 @@ The algebraic equations of the exciter. \begin{aligned} 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 &= -V_E + k_\text{sat} \\ 0 &= -E_{fd} + (1 + \omega I_{\mathrm{spdlim}})E_{fd}' \\ 0 &= -k_\text{sat} + S_B\, q(E_{fd}' - S_A) \end{aligned} @@ -190,7 +201,12 @@ with the current input values. 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}' \\ + K_E &\leftarrow + \begin{cases} + \dfrac{V_R^{\max}/10-k_\text{sat}}{E_{fd}'} & K_E=0 \\ + K_E & K_E\ne 0 + \end{cases} \\ + V_E &= k_\text{sat} \\ V_R &= K_E\, E_{fd}' + V_E \\ V_{tr} &= \dfrac{V_R}{K_A} \\ V_{fx} &= \dfrac{K_F}{T_F}\, E_{fd}' \\ @@ -207,4 +223,4 @@ All internal derivatives initialize to zero. Variable | Units | Description | Note ---------|--------|-----------------------------------|------ `efd` | [p.u.] | Field winding voltage | -`ksat` | [p.u.] | Magnetic saturation coefficient | $S_B\,q(E_{fd}'-S_A)$ +`ksat` | [p.u.] | Scaled-quadratic saturation contribution | $S_B\,q(E_{fd}'-S_A)$ diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/README.md b/GridKit/Model/PhasorDynamics/Governor/Tgov1/README.md index b03b9bd43..0315ffe2b 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/README.md +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/README.md @@ -12,16 +12,45 @@ Figure 1: Governor TGOV1 model. Figure courtesy of [PowerWorld](https://www.powe Symbol | Units | Description | Typical Value | Note ------------|--------|-----------------------------------|---------------| ------ -$T_{\mathrm{rate}}$ | [MW] | Turbine rating | 100.0 | -$R$ | [p.u.] | Droop Constant | 0.05 | -$T_1$ | [sec] | Valve Time Delay | 0.5 | -$T_2$ | [sec] | Turbine Numerator Time Constant | 2.5 | -$T_3$ | [sec] | Turbine Delay | 7.5 | -$P_v^{\max}$ | [p.u.] | Max Valve Position | 1 | -$P_v^{\min}$ | [p.u.] | Min Valve Position | 0 | -$D_t$ | [p.u.] | Turbine Damping Coefficient | 0 | +$T_{\mathrm{rate}}$ | [MVA] | Governor component power base | 100.0 | +$R$ | [p.u.] | Permanent droop | 0.05 | +$T_1$ | [sec] | Steam-bowl time constant | 0.5 | +$T_2$ | [sec] | Turbine numerator time constant | 2.5 | +$T_3$ | [sec] | Reheater time constant | 7.5 | +$P_v^{\max}$ | [p.u.] | Maximum valve position | 1 | +$P_v^{\min}$ | [p.u.] | Minimum valve position | 0 | +$D_t$ | [p.u.] | Turbine damping coefficient | 0 | -## Model Variables +### Parameter Validation + +The component and system power bases must be positive, $R$ must be nonzero, +and $P_v^{\min}\le P_v^{\max}$. + +Set $T_{\mathrm{rate}}$ equal to the connected machine MVA base. A zero +component power base is not supported. + +### Model Derived Parameters + +Let $\epsilon_T=10^{-3}\ \mathrm{s}$. A time constant below $\epsilon_T$ is +raised to that floor in place, so every equation below uses the raised value: + +```math +\begin{aligned} + T_x + &\leftarrow \max\!\left(T_x,\epsilon_T\right), + \quad x\in\{1,3\} +\end{aligned} +``` + +## Model Ports + +Name | Port | Init | Description +--------|--------|---------|------------ +`speed` | Input | Known | Machine speed deviation; optional, defaults to zero +`pref` | Input | Unknown | Governor reference; optional +`pmech` | Output | Known | Mechanical-power signal seeded by the machine + +## Model Variables ### Internal Variables @@ -29,82 +58,80 @@ $D_t$ | [p.u.] | Turbine Damping Coefficient | 0 | Symbol | Units | Description | Note ----------|--------|-----------------------------------|------- -$P_{tx}$ | [p.u.] | Turbine Power (State 1 in Fig. 1) | -$P_v$ | [p.u.] | Valve Position (State 2 in Fig. 1)| +$P_t$ | [p.u.] | Turbine-block output | Component base +$P_v$ | [p.u.] | Valve position | Component base #### Algebraic + Symbol | Units | Description | Note ----------------|--------|-----------------------------------|------- -$P_m$ | [p.u.] | Mechnical Power to Generator | Read by a Machine Model +$P_m$ | [p.u.] | Mechanical-power output | System base; read by the machine model ### External Variables #### Differential + Symbol | Units | Description | Note ----------------|--------|-----------------------------------|------- -$\omega$ | [p.u.] | Machine Speed Deviation | Read from a Machine Model +$\omega$ | [p.u.] | Machine speed deviation | Optional `speed` input; defaults to zero #### Algebraic + Symbol | Units | Description | Note ----------------|--------|-----------------------------------|------- -$P_{ref}$ | [p.u.] | Reference Power | Either a constant parameter or external variable +$P_\mathrm{ref}$ | [p.u.] | Governor reference | Component base; optional `pref` input, otherwise held internally ## Model Equations For readability, define: + ```math -f = -P_v + \dfrac{1}{R}(P_{ref} - \omega) +g_v=-P_v+\dfrac{P_\mathrm{ref}-\omega}{R}. ``` ### Differential Equations -The TGOV1 differential equations, as derived from the model diagram. + +The TGOV1 differential equations, as derived from the model diagram, are ```math \begin{aligned} - 0 &= -T_3 \dot P_{tx} - P_{tx} + (T_3 - T_2)P_v \\ - 0 &= -T_1 \dot P_v - + \text{antiwindup}( - P_v, - f, - P_v^{\min}, - P_v^{\max} - ) + 0 &= -\dot P_v + + \dfrac{1}{T_1}\text{antiwindup} + \left(P_v,g_v;P_v^{\min},P_v^{\max}\right) \\ + 0 &= -\dot P_t-\dfrac{P_t-P_v-T_2\dot P_v}{T_3}. \end{aligned} ``` -CommonMath defines the [Anti-Windup](../../../../CommonMath.md#antiwindup) +CommonMath defines the [anti-windup](../../../../CommonMath.md#antiwindup) target and smooth approximation. ### Algebraic Equations -The algebraic equation dictating the mechnical power output. + +The mechanical-power output is given by + ```math -\begin{aligned} - 0 &= -\dfrac{S_{\mathrm{sys}}}{T_{\mathrm{rate}}} P_m - + \dfrac{1}{T_3}(P_{tx}+T_2P_v) - D_t \omega \\ -\end{aligned} +0=-\dfrac{S_\mathrm{sys}}{T_\mathrm{rate}}P_m + +P_t-D_t\omega. ``` ## Initialization -At steady state we assume that $P_v$ is at or within its limits. This implies the initial conditions are a function of the initial mechanical power converted to the TGOV1 component base. -```math -\begin{aligned} - P^{\mathrm{tgov1}}_{m,0} - &= \dfrac{S_{\mathrm{sys}}}{T_{\mathrm{rate}}}P_{m,0} \\ - P_{tx,0} - &= (T_3 - T_2)P^{\mathrm{tgov1}}_{m,0} \\ - P_{v,0} - &= P^{\mathrm{tgov1}}_{m,0} \\ - \dot P_{tx,0} - &= 0 \\ - \dot P_{v,0} - &= 0 -\end{aligned} -``` -And if the reference power is a constant parameter, we can determine the value by solving the steady state equations. +TGOV1 preserves the machine-provided $P_{m,0}$ and initializes the steady +state in dependency order: + ```math \begin{aligned} - P_{ref,0} - &= R P^{\mathrm{tgov1}}_{m,0} + P_{m,0}^{\mathrm{TGOV1}} + &\leftarrow \dfrac{S_\mathrm{sys}}{T_\mathrm{rate}}P_{m,0} \\ + P_{v,0} + &\leftarrow P_{m,0}^{\mathrm{TGOV1}}+D_t\omega_0 \\ + P_{t,0} + &\leftarrow P_{v,0} \\ + P_{\mathrm{ref},0} + &\leftarrow \omega_0+RP_{v,0} \\ + \dot P_{v,0},\dot P_{t,0} + &\leftarrow 0. \end{aligned} ``` + +Initialization rejects $P_{v,0}$ outside the configured valve limits. diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp index bb8abc988..89ce12994 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp @@ -3,7 +3,7 @@ * @author Luke Lowery (lukel@tamu.edu) * @author Adam Birchfield (abirchfield@tamu.edu) * @author Wiktoria Zielinska (zielinskawa@ORNL.gov) - * @brief Declaration of a Turbine Governor Model (IEEET1). + * @brief Declaration of the TGOV1 turbine-governor model. * */ @@ -41,17 +41,17 @@ namespace GridKit /// Internal variables of a `Tgov1` enum class Tgov1InternalVariables : size_t { - PTX, ///< \f$P_{tx}\f$ - PV, ///< \f$P_v\f$ - PM, ///< \f$P_m\f$ + PTX, ///< \f$P_t\f$ Turbine-block output + PV, ///< \f$P_v\f$ Valve position + PM, ///< \f$P_m\f$ Mechanical-power output MAXIMUM, }; /// External variables of a `Tgov1` enum class Tgov1ExternalVariables : size_t { - DELTAOMEGA, ///< \f$\Delta_\omega\f$ - PREF, ///< \f$P_{ref}\f$ + DELTAOMEGA, ///< \f$\omega\f$ Machine speed deviation + PREF, ///< \f$P_\mathrm{ref}\f$ Governor reference MAXIMUM, }; @@ -117,20 +117,20 @@ namespace GridKit private: // Input parameters - RealT Trate_{0}; - RealT R_{0}; - RealT Pvmin_{0}; - RealT Pvmax_{0}; - RealT T1_{0}; - RealT T2_{0}; - RealT T3_{0}; - RealT Dt_{0}; + RealT Trate_{static_cast(100.0)}; + RealT R_{static_cast(0.05)}; + RealT Pvmin_{ZERO}; + RealT Pvmax_{ONE}; + RealT T1_{static_cast(0.5)}; + RealT T2_{static_cast(2.5)}; + RealT T3_{static_cast(7.5)}; + RealT Dt_{ZERO}; // Derived parameters RealT va_component_base_{0}; // Input States (which can be parameters) - ScalarT pref_{0}; + ScalarT pref_set_{0}; /// Component signal extension ComponentSignals signals_; @@ -141,6 +141,8 @@ namespace GridKit ScalarT toComponentBase(ScalarT value) const; ScalarT toSystemBase(ScalarT value) const; + static constexpr RealT TIME_CONSTANT_MINIMUM = static_cast(1.0e-3); + /* Local copies of signal variables */ std::vector ws_; std::vector ws_indices_; diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Data.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Data.hpp index 2052572de..46ab373ed 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Data.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Data.hpp @@ -22,14 +22,14 @@ namespace GridKit */ enum class Tgov1Parameters { - Trate, ///< Turbine-rating power base - R, ///< Droop Constant - T1, ///< Valve Time Delay - T2, ///< Turbine Numerator Time Constant - T3, ///< Turbine Delay - Pvmax, ///< Max Valve Power - Pvmin, ///< Min Valve Power - Dt ///< Damping Coefficient + Trate, ///< \f$T_\mathrm{rate}\f$ Turbine-rating power base + R, ///< \f$R\f$ Permanent droop + T1, ///< \f$T_1\f$ Steam-bowl time constant + T2, ///< \f$T_2\f$ Turbine numerator time constant + T3, ///< \f$T_3\f$ Reheater time constant + Pvmax, ///< \f$P_v^\mathrm{max}\f$ Maximum valve position + Pvmin, ///< \f$P_v^\mathrm{min}\f$ Minimum valve position + Dt ///< \f$D_t\f$ Turbine damping coefficient }; /** @@ -45,7 +45,8 @@ namespace GridKit */ enum class Tgov1SignalInputs : size_t { - speed, + speed, ///< \f$\omega\f$ Optional machine speed-deviation signal ID + pref, ///< \f$P_\mathrm{ref}\f$ Optional governor-reference signal ID SIZE, }; @@ -54,7 +55,7 @@ namespace GridKit */ enum class Tgov1SignalOutputs : size_t { - pmech, + pmech, ///< \f$P_m\f$ Required mechanical-power output signal ID SIZE, }; diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp index ee5bed430..1e582930a 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp @@ -1,14 +1,15 @@ #pragma once /** - * @file Tgov1Impl.cpp + * @file Tgov1Impl.hpp * @author Luke Lowery (lukel@tamu.edu) * @author Adam Birchfield (abirchfield@tamu.edu) * @author Wiktoria Zielinska (zielinskawa@ORNL.gov) - * @brief Definition of a Turbine Governor Model (IEEET1). + * @brief Definition of the TGOV1 turbine-governor model. */ -#include +#include +#include #include #include @@ -39,7 +40,7 @@ namespace GridKit * @brief Constructs a Tgov1 governor model from its parameters * * @param pmech $P_m$ internal variable signal node - * @param omega $\Delta_\omega$ external variable signal node + * @param omega $\omega$ external variable signal node */ template Tgov1::Tgov1(SignalT* pmech, SignalT* omega) @@ -130,6 +131,9 @@ namespace GridKit template void Tgov1::setDerivedParams() { + T1_ = std::max(T1_, TIME_CONSTANT_MINIMUM); + T3_ = std::max(T3_, TIME_CONSTANT_MINIMUM); + va_component_base_ = Trate_ * static_cast(1.0e6); } @@ -183,16 +187,17 @@ namespace GridKit } // Resize signal variable data - ws_.resize(1); - ws_indices_.resize(1); - ws_[0] = 0.0; - ws_indices_[0] = INVALID_INDEX; + const auto signal_size = static_cast(Tgov1ExternalVariables::MAXIMUM); + ws_.assign(signal_size, ScalarT{0}); + ws_indices_.assign(signal_size, INVALID_INDEX); // Set output signals if (signals_.template isAssigned()) { - auto* y = y_.getData(); - signals_.template getSignalNode()->set(&y[2], &(this->getVariableIndex(2))); + auto* y = y_.getData(); + const auto PM = static_cast(Tgov1InternalVariables::PM); + signals_.template getSignalNode()->set( + &y[static_cast(PM)], &(this->getVariableIndex(PM))); } allocated_ = true; @@ -206,14 +211,40 @@ namespace GridKit int Tgov1::verify() const { static constexpr auto DELTAOMEGA = Tgov1ExternalVariables::DELTAOMEGA; + static constexpr auto PREF = Tgov1ExternalVariables::PREF; int ret = 0; + auto check = [&](bool condition, const char* message) + { + if (!condition) + { + Log::error() << "Tgov1: " << message << '\n'; + ret += 1; + } + }; + + check(Trate_ > ZERO, "Trate must be positive"); + check(va_system_base_ > ZERO, "system power base must be positive"); + check(R_ != ZERO, "R must be nonzero"); + check(Pvmin_ <= Pvmax_, "Pvmin must be less than or equal to Pvmax"); + check(signals_.template isAssigned(), + "pmech output signal must be assigned"); + if (signals_.template isAttached()) { if (!signals_.template isLinked()) { - Log::error() << "Tgov1: deltaomega signal attached with no linked generator\n"; + Log::error() << "Tgov1: speed signal attached with no linked generator\n"; + ret += 1; + } + } + + if (signals_.template isAttached()) + { + if (!signals_.template isLinked()) + { + Log::error() << "Tgov1: pref signal attached with no linked source\n"; ret += 1; } } @@ -228,32 +259,49 @@ namespace GridKit template int Tgov1::initialize() { - ScalarT p0{0}; + if (verify() != 0) + { + Log::error() << "Tgov1: cannot initialize with invalid configuration\n"; + return 1; + } - // Initial mechanical = initial electric torque - auto* y = y_.getData(); - auto* yp = yp_.getData(); - if (signals_.template isAssigned()) + const auto PTX = static_cast(Tgov1InternalVariables::PTX); + const auto PV = static_cast(Tgov1InternalVariables::PV); + const auto PM = static_cast(Tgov1InternalVariables::PM); + + auto* y = y_.getData(); + + ScalarT omega0{ZERO}; + if (signals_.template isAttached()) { - // System base -> governor base for governor initialization. - p0 = toComponentBase(y[2]); ///<- generator needs to be initialized first + omega0 = signals_.template readExternalVariable(); } - // Input Variables (Parameter for now) - pref_ = R_ * p0; + const ScalarT pmech0 = y[PM]; + const ScalarT pm0 = toComponentBase(pmech0); + const ScalarT pv0 = pm0 + Dt_ * omega0; + const ScalarT pturb0 = pv0; + const ScalarT pref0 = omega0 + R_ * pv0; + + const RealT pv0_value = static_cast(pv0); + const RealT limit_tolerance = static_cast(4.0) * std::numeric_limits::epsilon(); + if (pv0_value < Pvmin_ - limit_tolerance || pv0_value > Pvmax_ + limit_tolerance) + { + Log::error() << "Tgov1: initial valve position is outside limits\n"; + return 1; + } - // Internal States - y[0] = (T3_ - T2_) * p0; // y0 - Ptx (Turbine Power ) - y[1] = p0; // y1 - Pv (Valve Position) - y[2] = toSystemBase(p0); // y2 - Pm (Mech Power, System Base) + y[PTX] = pturb0; + y[PV] = pv0; - // D.V. Derivative - yp[0] = 0.0; // Ptx - yp[1] = 0.0; // Pv - yp[2] = 0.0; // Pm + pref_set_ = pref0; + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(pref_set_); + } y_.setDataUpdated(); - yp_.setDataUpdated(); + yp_.setToConst(static_cast(ZERO)); return 0; } @@ -264,8 +312,8 @@ namespace GridKit template int Tgov1::tagDifferentiable() { - tag_[0] = true; // Pv - tag_[1] = true; // Ptx + tag_[0] = true; // Ptx + tag_[1] = true; // Pv tag_[2] = false; // Pmech return 0; @@ -302,28 +350,26 @@ namespace GridKit const ScalarT* ws, ScalarT* f) { - // Read Internal Variables - ScalarT ptx = y[0]; // y0 - Ptx - ScalarT pv = y[1]; // y1 - Pv - ScalarT pmech = y[2]; // y2 - Pmech + const auto PTX = static_cast(Tgov1InternalVariables::PTX); + const auto PV = static_cast(Tgov1InternalVariables::PV); + const auto PM = static_cast(Tgov1InternalVariables::PM); - // Read Internal Derivatives - ScalarT ptx_dot = yp[0]; - ScalarT pv_dot = yp[1]; + const auto DELTAOMEGA = static_cast(Tgov1ExternalVariables::DELTAOMEGA); + const auto PREF = static_cast(Tgov1ExternalVariables::PREF); - // Set signal variable aliases - ScalarT omega = ws[0]; + const ScalarT pturb = y[PTX]; + const ScalarT pv = y[PV]; + const ScalarT pmech = y[PM]; - // The 'pre-limit' target of Pv - ScalarT func = -pv + (pref_ - omega) / R_; + const ScalarT pturb_dot = yp[PTX]; + const ScalarT pv_dot = yp[PV]; - // Internal Differential Equations - f[0] = -T3_ * ptx_dot - ptx + (T3_ - T2_) * pv; - f[1] = -T1_ * pv_dot + Math::antiwindup(pv, func, Pvmin_, Pvmax_); + const ScalarT omega = ws[DELTAOMEGA]; + const ScalarT pref = ws[PREF]; - // Internal Algebraic Equations - // Convert pmech to component base from its system base value - f[2] = -toComponentBase(pmech) + (ptx + T2_ * pv) / T3_ - (Dt_ * omega); + f[PTX] = -pturb_dot - (pturb - pv - T2_ * pv_dot) / T3_; + f[PV] = -pv_dot + Math::antiwindup(pv, -pv + (pref - omega) / R_, Pvmin_, Pvmax_) / T1_; + f[PM] = -toComponentBase(pmech) + pturb - Dt_ * omega; return 0; } @@ -335,11 +381,25 @@ namespace GridKit template int Tgov1::evaluateResidual() { - // Input Variables + const auto DELTAOMEGA = static_cast(Tgov1ExternalVariables::DELTAOMEGA); + const auto PREF = static_cast(Tgov1ExternalVariables::PREF); + + ws_[DELTAOMEGA] = ScalarT{ZERO}; + ws_[PREF] = pref_set_; + std::fill(ws_indices_.begin(), ws_indices_.end(), INVALID_INDEX); + if (signals_.template isAttached()) { - ws_[0] = signals_.template readExternalVariable(); - ws_indices_[0] = signals_.template readExternalVariableIndex(); + ws_[DELTAOMEGA] = signals_.template readExternalVariable(); + ws_indices_[DELTAOMEGA] = + signals_.template readExternalVariableIndex(); + } + + if (signals_.template isAttached()) + { + ws_[PREF] = signals_.template readExternalVariable(); + ws_indices_[PREF] = + signals_.template readExternalVariableIndex(); } const auto* y = y_.getData(); diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/Genrou.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/Genrou.hpp index 75e765b24..1dac5ac96 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/Genrou.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/Genrou.hpp @@ -38,21 +38,21 @@ namespace GridKit { DELTA, ///< \f$\delta\f$ OMEGA, ///< \f$\omega\f$ + EPQ, ///< \f$E'_q\f$ PSIPD, ///< \f$\psi'_d\f$ PSIPQ, ///< \f$\psi'_q\f$ EPD, ///< \f$E'_d\f$ - EPQ, ///< \f$E'_q\f$ + PSIPPQ, ///< \f$\psi''_q\f$ + PSIPPD, ///< \f$\psi''_d\f$ + PSIPP, ///< \f$\psi''\f$ + KSAT, ///< \f$k_{sat}\f$ VD, ///< \f$V_d\f$ VQ, ///< \f$V_q\f$ + TE, ///< \f$T_e\f$ ID, ///< \f$I_d\f$ IQ, ///< \f$I_q\f$ IR, ///< \f$I_r\f$ II, ///< \f$I_i\f$ - PSIPPQ, ///< \f$\psi''_q\f$ - PSIPPD, ///< \f$\psi''_d\f$ - PSIPP, ///< \f$\psi''\f$ - TE, ///< \f$T_e\f$ - KSAT, ///< \f$k_{sat}\f$ MAXIMUM, }; diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouImpl.hpp index 2d933043d..2169ea443 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/GenrouImpl.hpp @@ -397,94 +397,61 @@ namespace GridKit template int Genrou::initialize() { - // Saturated initialization with ksat iteration. - // See README "With Saturation" and plan for derivation. - // Network Frame Terminal Values ScalarT vr = Vr(); ScalarT vi = Vi(); ScalarT p = toMachineBase(static_cast(p0_)); ScalarT q = toMachineBase(static_cast(q0_)); ScalarT vm2 = vr * vr + vi * vi; - ScalarT vm = std::sqrt(vm2); ScalarT ir = (p * vr + q * vi) / vm2; ScalarT ii = (p * vi - q * vr) / vm2; - // Initial ksat guess from |V| - ScalarT vm_sat = vm - SA_; - ScalarT ksat = (vm_sat > ZERO) ? SB_ * vm_sat * vm_sat : ScalarT{ZERO}; - - ScalarT delta, id, iq, vd, vq; - ScalarT psiqpp, psidpp, psipp; - ScalarT Edp, psiqp, psidp, Eqp; - - static constexpr int max_iter = 10; - static constexpr RealT ksat_tol = 1e-12; - ScalarT ksat_prev = ksat + ONE; // force first iteration - - for (int iter = 0; iter < max_iter; ++iter) - { - // Convergence check (skip first iteration) - ScalarT ksat_err = ksat - ksat_prev; - if (iter > 0 && ksat_err * ksat_err < ksat_tol * ksat_tol) - break; - ksat_prev = ksat; - - // Compute all variables consistent with current ksat - ScalarT ksat_prime = ONE + Xqd_ * ksat; - ScalarT Xsat_delta = ksat_prime * Xdpp_ + Xq_ - Xqpp_; - - delta = std::atan2((vi + Ra_ * ii) * ksat_prime + Xsat_delta * ir, - (vr + Ra_ * ir) * ksat_prime - Xsat_delta * ii); - id = ir * std::sin(delta) - ii * std::cos(delta); - iq = ir * std::cos(delta) + ii * std::sin(delta); - vd = vr * std::sin(delta) - vi * std::cos(delta) + id * Ra_ - iq * Xqpp_; - vq = vr * std::cos(delta) + vi * std::sin(delta) + id * Xqpp_ + iq * Ra_; - psiqpp = -vd; - psidpp = vq; - Edp = (Xq1_ - Xqd_ * (Xqp_ - Xqpp_) * ksat) * iq / (ONE + Xqd_ * ksat); - psiqp = Edp + Xq2_ * iq; - psidp = psidpp - (Xdpp_ - Xl_) * id; - Eqp = psidp + (Xdp_ - Xl_) * id; - - // Update ksat from flux-linkage psipp for next iteration - ScalarT psiqpp_fl = -psiqp * Xq4_ - Edp * Xq5_; - ScalarT psidpp_fl = psidp * Xd4_ + Eqp * Xd5_; - psipp = std::sqrt(psiqpp_fl * psiqpp_fl + psidpp_fl * psidpp_fl); - ScalarT psipp_sat = psipp - SA_; - ksat = (psipp_sat > ZERO) ? SB_ * psipp_sat * psipp_sat : ScalarT{ZERO}; - - if (iter == max_iter - 1) - { - Log::warning() << "Genrou: saturated initialization did not converge" - << " (ksat_err=" << ksat_err << ")\n"; - } - } + // The subtransient-flux magnitude is invariant under the rotor-frame + // rotation, so saturation is available directly from network quantities. + const ScalarT Vint_r = vr + Ra_ * ir - Xqpp_ * ii; + const ScalarT Vint_i = vi + Ra_ * ii + Xqpp_ * ir; + ScalarT psipp = std::sqrt(Vint_r * Vint_r + Vint_i * Vint_i); + ScalarT ksat = SB_ * Math::qramp(psipp - SA_); + + const ScalarT ksat_prime = ONE + Xqd_ * ksat; + const ScalarT xsat_delta = ksat_prime * Xdpp_ + Xq_ - Xqpp_; + + ScalarT delta = std::atan2((vi + Ra_ * ii) * ksat_prime + xsat_delta * ir, + (vr + Ra_ * ir) * ksat_prime - xsat_delta * ii); + ScalarT id = ir * std::sin(delta) - ii * std::cos(delta); + ScalarT iq = ir * std::cos(delta) + ii * std::sin(delta); + ScalarT vq = vr * std::cos(delta) + vi * std::sin(delta) + id * Xqpp_ + iq * Ra_; + ScalarT Edp = (Xq1_ - Xqd_ * (Xqp_ - Xqpp_) * ksat) * iq / ksat_prime; + ScalarT psiqp = Edp + Xq2_ * iq; + ScalarT psidp = vq - (Xdpp_ - Xl_) * id; + ScalarT Eqp = psidp + (Xdp_ - Xl_) * id; - // Assign from converged values using flux-linkage forms ScalarT omega(0.0); auto* y = y_.getData(); auto* yp = yp_.getData(); - y[0] = delta; - y[1] = omega; - y[2] = Eqp; - y[3] = psidp; - y[4] = psiqp; - y[5] = Edp; - y[6] = psiqpp = -psiqp * Xq4_ - Edp * Xq5_; - y[7] = psidpp = psidp * Xd4_ + Eqp * Xd5_; - y[8] = psipp = std::sqrt(psiqpp * psiqpp + psidpp * psidpp); - ScalarT psipp_sat = psipp - SA_; - y[9] = ksat = (psipp_sat > ZERO) ? SB_ * psipp_sat * psipp_sat : ScalarT{ZERO}; - y[10] = vd = -psiqpp * (ONE + omega); - y[11] = vq = psidpp * (ONE + omega); - y[12] = (psidpp - id * Xdpp_) * iq - (psiqpp - iq * Xdpp_) * id; - y[13] = id; - y[14] = iq; - y[15] = ir; - y[16] = ii; - y[17] = G_ * (vd * std::sin(delta) + vq * std::cos(delta)) + y[0] = delta; + y[1] = omega; + y[2] = Eqp; + y[3] = psidp; + y[4] = psiqp; + y[5] = Edp; + ScalarT psiqpp = -psiqp * Xq4_ - Edp * Xq5_; + ScalarT psidpp = psidp * Xd4_ + Eqp * Xd5_; + y[6] = psiqpp; + y[7] = psidpp; + y[8] = psipp = std::sqrt(psiqpp * psiqpp + psidpp * psidpp); + y[9] = ksat = SB_ * Math::qramp(psipp - SA_); + ScalarT vd = -psiqpp * (ONE + omega); + vq = psidpp * (ONE + omega); + y[10] = vd; + y[11] = vq; + y[12] = (psidpp - id * Xdpp_) * iq - (psiqpp - iq * Xdpp_) * id; + y[13] = id; + y[14] = iq; + y[15] = ir; + y[16] = ii; + y[17] = G_ * (vd * std::sin(delta) + vq * std::cos(delta)) - B_ * (vd * -std::cos(delta) + vq * std::sin(delta)); y[18] = B_ * (vd * std::sin(delta) + vq * std::cos(delta)) + G_ * (vd * -std::cos(delta) + vq * std::sin(delta)); @@ -606,18 +573,17 @@ namespace GridKit f[5] = Edp_dot - (ONE / Tqop_) * (-Edp + Xqd_ * psiqpp * ksat + Xq1_ * (iq - Xq3_ * (Edp + iq * Xq2_ - psiqp))); /* 11 Genrou algebraic equations */ - f[6] = psiqpp - (-psiqp * Xq4_ - Edp * Xq5_); - f[7] = psidpp - (psidp * Xd4_ + Eqp * Xd5_); - f[8] = psipp - std::sqrt((psidpp * psidpp) + (psiqpp * psiqpp)); - ScalarT psipp_sat = psipp - SA_; - f[9] = ksat - SB_ * psipp_sat * psipp_sat * Math::sigmoid(psipp_sat); - f[10] = vd + psiqpp * (ONE + omega); - f[11] = vq - psidpp * (ONE + omega); - f[12] = telec - ((psidpp - id * Xdpp_) * iq - (psiqpp - iq * Xdpp_) * id); - f[13] = id - (ir * std::sin(delta) - ii * std::cos(delta)); - f[14] = iq - (ir * std::cos(delta) + ii * std::sin(delta)); - f[15] = ir + G_ * vr - B_ * vi - inr; - f[16] = ii + B_ * vr + G_ * vi - ini; + f[6] = psiqpp - (-psiqp * Xq4_ - Edp * Xq5_); + f[7] = psidpp - (psidp * Xd4_ + Eqp * Xd5_); + f[8] = psipp - std::sqrt((psidpp * psidpp) + (psiqpp * psiqpp)); + f[9] = ksat - SB_ * Math::qramp(psipp - SA_); + f[10] = vd + psiqpp * (ONE + omega); + f[11] = vq - psidpp * (ONE + omega); + f[12] = telec - ((psidpp - id * Xdpp_) * iq - (psiqpp - iq * Xdpp_) * id); + f[13] = id - (ir * std::sin(delta) - ii * std::cos(delta)); + f[14] = iq - (ir * std::cos(delta) + ii * std::sin(delta)); + f[15] = ir + G_ * vr - B_ * vi - inr; + f[16] = ii + B_ * vr + G_ * vi - ini; /* 2 Genrou current source definitions */ f[17] = inr - (G_ * (std::sin(delta) * vd + std::cos(delta) * vq) - B_ * (-std::cos(delta) * vd + std::sin(delta) * vq)); diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/README.md b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/README.md index 96679af3f..4150b49a9 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/README.md +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENROU/README.md @@ -136,7 +136,7 @@ Note that for implementation purposes, some of these equations may be simplified 0 &= -V_{d} -\psi''_{q}(1+\omega)\\ 0 &= -V_{q} +\psi''_{d}(1+\omega)\\ 0 &= -T_{elec} +(\psi''_{d} - I_dX_d'')I_q-(\psi''_{q} - I_qX_d'')I_d \\ - 0 &= -k_{sat} + S_B(\psi''-S_A)^2\sigma(\psi''-S_A) \\ + 0 &= -k_{sat} + S_B q(\psi''-S_A) \\ 0 &= -I_d + I_r \sin(\delta) - I_i \cos(\delta) \\ 0 &= -I_q + I_r \cos(\delta) + I_i \sin(\delta) \\ 0 &= -I_r + G (V_d \sin(\delta) + V_q \cos(\delta) - V_r) - B (V_d \cos(\delta) + V_q \sin(\delta) - V_i) \\ @@ -144,54 +144,39 @@ Note that for implementation purposes, some of these equations may be simplified \end{aligned} ``` +CommonMath defines the primitive +[quadratic ramp](../../../../CommonMath.md#primitives) $q$. + ## Initialization -### Without Saturation -Presume there is no saturation to simplify the solution procedure for the initial -conditions. +The power-flow solution gives $V_r$, $V_i$, $I_r$, and $I_i$. At synchronous +speed, the total subtransient-flux magnitude is available directly in the +network frame and is independent of rotor angle: -Using the power-flow solution, we have explicit solutions for the following -variables. The internal variables $I_d$, $I_q$, $V_d$, and $V_q$ are calculated -from the network interface equations. The remaining are algebraically solved -from the steady-state initial conditions. ``` math \begin{aligned} -\omega &= 0 \\ -\delta &= \text{arg} \left[V_r + jV_i + (R_a + jX_q) (I_r + jI_i)\right] \\ - \psi^{''}_{d} &= V_q \\ - \psi^{''}_{q} &= -V_d \\ - \psi^{''} &= \sqrt{(\psi''_{d})^2+(\psi''_{q})^2} \\ - k_{sat} &= - \begin{cases} - S_B(\psi^{''}-S_A)^2, & \psi^{''} > S_A\\ - 0, & \psi^{''} \le S_A - \end{cases}\\ - T_{elec} &= (\psi''_{d} - I_dX_d^{''})I_q-(\psi''_{q} - I_qX_d^{''})I_d \\ - P_{m} &= T_{elec} \\ - \psi_d' &= \psi_d'' - (X_d'' - X_\ell)I_d \\ - \psi_q' &= (X_q'' - X_\ell)I_q - \psi_q'' \\ - E^{'}_d &=\psi^{'}_q - X_{q2}I_q \\ - E^{'}_q &=\psi^{'}_d + X_{d2}I_d \\ - E_{fd} &= E'_{q}+X_{d1}I_{d}+\psi^{''}_{d}k_{sat} \\ + \psi'' + &\leftarrow \sqrt{ + \left(V_r + R_a I_r - X_q'' I_i\right)^2 + +\left(V_i + R_a I_i + X_q'' I_r\right)^2 + } \\ + k_{sat} + &\leftarrow S_B q(\psi''-S_A) \\ + k'_{sat} + &\leftarrow 1 + X_{qd}k_{sat} \\ + X_{sat,\delta} + &\leftarrow k'_{sat}X_d'' + X_q-X_q'' \\ + \delta + &\leftarrow \operatorname{atan2}\left( + (V_i+R_aI_i)k'_{sat}+X_{sat,\delta}I_r, + (V_r+R_aI_r)k'_{sat}-X_{sat,\delta}I_i + \right) \end{aligned} ``` -### With Saturation -It is important to point out that finding the initial value of $\delta$ for -the model without saturation, the direct method can be used. In case saturation -is considered, some "clever" math is needed. Key insight for determining initial -$\delta$ is that the magnitude of the saturation, which depends upon the magnitude -of $\psi''$, which is independent of $\delta$. - -``` math -\begin{aligned} - \delta=\tan^{-1} - \left[ - \dfrac{(V_{i}+R_{a}I_{i})k_{sat}+(k_{sat}X''_{d}+X_{q}-X''_{q})I_{r}} - {(V_{r}+R_{a}I_{r})k_{sat}-(k_{sat}X''_{d}+X_{q}-X''_{q})I_{i}} - \right] -\end{aligned} -``` +With $\delta$ known, the rotor-frame currents, voltages, flux states, field +voltage, and mechanical power follow directly from the steady-state model +equations above. ## Model Outputs diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/Gensal.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/Gensal.hpp index 2060834ef..74655e7f1 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/Gensal.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/Gensal.hpp @@ -32,32 +32,32 @@ namespace GridKit /// Internal variables of a `Gensal` enum class GensalInternalVariables : size_t { - DELTA, ///< rotor angle - OMEGA, ///< speed deviation - EPQ, ///< q-axis transient voltage - PSIPD, ///< d-axis transient flux - PSIPPQ, ///< q-axis subtransient flux - PSIPPD, ///< d-axis subtransient flux - KSAT, ///< saturation signal - VD, ///< d-axis terminal voltage - VQ, ///< q-axis terminal voltage - TE, ///< electrical torque - ID, ///< d-axis current - IQ, ///< q-axis current - IR, ///< network real current - II, ///< network imaginary current - INR, ///< Norton source real current - INI, ///< Norton source imaginary current + DELTA, ///< \f$\delta\f$ rotor angle + OMEGA, ///< \f$\omega\f$ speed deviation + EPQ, ///< \f$E'_q\f$ q-axis transient voltage + PSIPD, ///< \f$\psi'_d\f$ d-axis transient flux + PSIPPQ, ///< \f$\psi''_q\f$ q-axis subtransient flux + PSIPPD, ///< \f$\psi''_d\f$ d-axis subtransient flux + KSAT, ///< \f$k_{sat}\f$ saturation factor + VD, ///< \f$V_d\f$ d-axis terminal voltage + VQ, ///< \f$V_q\f$ q-axis terminal voltage + TE, ///< \f$T_e\f$ electrical torque + ID, ///< \f$I_d\f$ d-axis current + IQ, ///< \f$I_q\f$ q-axis current + IR, ///< \f$I_r\f$ network real current + II, ///< \f$I_i\f$ network imaginary current + INR, ///< \f$I_{\mathrm{N},r}\f$ Norton source real current + INI, ///< \f$I_{\mathrm{N},i}\f$ Norton source imaginary current MAXIMUM, }; /// External variables of a `Gensal` enum class GensalExternalVariables : size_t { - VR, ///< network real voltage - VI, ///< network imaginary voltage - PM, ///< mechanical power - EFD, ///< field voltage + VR, ///< \f$V_r\f$ network real voltage + VI, ///< \f$V_i\f$ network imaginary voltage + PM, ///< \f$P_m\f$ mechanical power + EFD, ///< \f$E_{fd}\f$ field voltage MAXIMUM, }; diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalData.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalData.hpp index 35555665f..85dc16acc 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalData.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalData.hpp @@ -15,22 +15,22 @@ namespace GridKit /// Initial parameters for a Gensal generator model enum class GensalParameters { - p0, ///< Initial active power - q0, ///< Initial reactive power - H, ///< Rotor inertia - D, ///< Damping coefficient - Ra, ///< Winding resistance - Tdop, ///< Open circuit direct axis transient time - Tdopp, ///< Open circuit direct axis sub-transient time - Tqopp, ///< Open circuit quadrature axis sub-transient time - Xd, ///< Direct axis synchronous reactance - Xdp, ///< Direct axis transient reactance - Xdpp, ///< Direct axis sub-transient reactance - Xq, ///< Quadrature axis synchronous reactance - Xl, ///< Stator leakage reactance - S10, ///< Saturation factor at 1.0 pu flux - S12, ///< Saturation factor at 1.2 pu flux - mva, ///< MVA base of the gensal model + p0, ///< \f$P_0\f$ Initial active power + q0, ///< \f$Q_0\f$ Initial reactive power + H, ///< \f$H\f$ Rotor inertia + D, ///< \f$D\f$ Damping coefficient + Ra, ///< \f$R_a\f$ Winding resistance + Tdop, ///< \f$T'_{d0}\f$ Open circuit direct axis transient time + Tdopp, ///< \f$T''_{d0}\f$ Open circuit direct axis sub-transient time + Tqopp, ///< \f$T''_{q0}\f$ Open circuit quadrature axis sub-transient time + Xd, ///< \f$X_d\f$ Direct axis synchronous reactance + Xdp, ///< \f$X'_d\f$ Direct axis transient reactance + Xdpp, ///< \f$X''_d\f$ Direct axis sub-transient reactance + Xq, ///< \f$X_q\f$ Quadrature axis synchronous reactance + Xl, ///< \f$X_{\ell}\f$ Stator leakage reactance + S10, ///< \f$S_{10}\f$ Saturation factor at 1.0 pu flux + S12, ///< \f$S_{12}\f$ Saturation factor at 1.2 pu flux + mva, ///< \f$S_\mathrm{mach}\f$ MVA base of the gensal model }; /// Buses for a Gensal generator model @@ -43,37 +43,37 @@ namespace GridKit /// Signal inputs for a Gensal generator model enum class GensalSignalInputs : size_t { - pmech, ///< Unique ID of the signal providing mechanical power - efd, ///< Unique ID of the signal providing exciter field voltage + pmech, ///< \f$P_m\f$ Unique ID of the signal providing mechanical power + efd, ///< \f$E_{fd}\f$ Unique ID of the signal providing exciter field voltage SIZE }; /// Signal outputs for a Gensal generator model enum class GensalSignalOutputs : size_t { - speed, ///< Unique ID of the signal receiving speed deviation + speed, ///< \f$\omega\f$ Unique ID of the signal receiving speed deviation SIZE }; /// Variables able to be monitored for a Gensal generator model enum class GensalMonitorableVariables { - ir, - ii, - p, - q, - delta, - omega, - speed, - Eqp, - psidp, - psiqpp, - psidpp, - vd, - vq, - te, - id, - iq + ir, ///< \f$I_r\f$ Network-frame real terminal current + ii, ///< \f$I_i\f$ Network-frame imaginary terminal current + p, ///< \f$P\f$ Active power + q, ///< \f$Q\f$ Reactive power + delta, ///< \f$\delta\f$ Rotor angle + omega, ///< \f$\omega\f$ Speed deviation + speed, ///< \f$1+\omega\f$ Per-unit machine speed + Eqp, ///< \f$E'_q\f$ Q-axis transient voltage + psidp, ///< \f$\psi'_d\f$ D-axis transient flux + psiqpp, ///< \f$\psi''_q\f$ Q-axis subtransient flux + psidpp, ///< \f$\psi''_d\f$ D-axis subtransient flux + vd, ///< \f$V_d\f$ D-axis terminal voltage + vq, ///< \f$V_q\f$ Q-axis terminal voltage + te, ///< \f$T_e\f$ Electrical torque + id, ///< \f$I_d\f$ D-axis current + iq ///< \f$I_q\f$ Q-axis current }; /** diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalImpl.hpp b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalImpl.hpp index ebb3f1471..385eda4d6 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/GensalImpl.hpp @@ -286,17 +286,16 @@ namespace GridKit ScalarT delta = std::atan2(Ei, Er); ScalarT omega(0.0); - ScalarT id = ir * std::sin(delta) - ii * std::cos(delta); - ScalarT iq = ir * std::cos(delta) + ii * std::sin(delta); - ScalarT psiqpp = -Xq2_ * iq; - ScalarT vd = -psiqpp * (ONE + omega); - ScalarT vq = vr * std::cos(delta) + vi * std::sin(delta) + id * Xdpp_ + iq * Ra_; - ScalarT psidpp = vq / (ONE + omega); - ScalarT psidp = psidpp - (Xdpp_ - Xl_) * id; - ScalarT Eqp = psidp + Xd2_ * id; - ScalarT Eqp_sat = Eqp - SA_; - ScalarT ksat = SB_ * Eqp_sat * Eqp_sat * Math::sigmoid(Eqp_sat); - ScalarT Te = (psidpp - id * Xdpp_) * iq - (psiqpp - iq * Xdpp_) * id; + ScalarT id = ir * std::sin(delta) - ii * std::cos(delta); + ScalarT iq = ir * std::cos(delta) + ii * std::sin(delta); + ScalarT psiqpp = -Xq2_ * iq; + ScalarT vd = -psiqpp * (ONE + omega); + ScalarT vq = vr * std::cos(delta) + vi * std::sin(delta) + id * Xdpp_ + iq * Ra_; + ScalarT psidpp = vq / (ONE + omega); + ScalarT psidp = psidpp - (Xdpp_ - Xl_) * id; + ScalarT Eqp = psidp + Xd2_ * id; + ScalarT ksat = SB_ * Math::qramp(Eqp - SA_); + ScalarT Te = (psidpp - id * Xdpp_) * iq - (psiqpp - iq * Xdpp_) * id; auto* y = y_.getData(); auto* yp = yp_.getData(); @@ -327,7 +326,7 @@ namespace GridKit signals_.template writeExternalVariable(pmech_set_); } - efd_set_ = Eqp + Xd1_ * (id + Xd3_ * (Eqp - psidp - Xd2_ * id)) + ksat; + efd_set_ = Eqp + Xd1_ * (id + Xd3_ * (Eqp - psidp - Xd2_ * id)) + Eqp * ksat; if (signals_.template isAttached()) { signals_.template writeExternalVariable(efd_set_); @@ -426,21 +425,20 @@ namespace GridKit /* 5 Gensal differential equations */ f[0] = delta_dot - omega * (TWO * pi * freq_system_base_); f[1] = omega_dot - (ONE / (TWO * H_)) * ((pmech - D_ * omega) / (ONE + omega) - telec); - f[2] = Eqp_dot - (ONE / Tdop_) * (efd - (Eqp + Xd1_ * (id + Xd3_ * (Eqp - psidp - Xd2_ * id)) + ksat)); + f[2] = Eqp_dot - (ONE / Tdop_) * (efd - (Eqp + Xd1_ * (id + Xd3_ * (Eqp - psidp - Xd2_ * id)) + Eqp * ksat)); f[3] = psidp_dot - (ONE / Tdopp_) * (Eqp - psidp - Xd2_ * id); f[4] = psiqpp_dot - (ONE / Tqopp_) * (-psiqpp - Xq2_ * iq); /* 9 Gensal algebraic equations */ - f[5] = psidpp - (psidp * Xd4_ + Eqp * Xd5_); - ScalarT Eqp_sat = Eqp - SA_; - f[6] = ksat - SB_ * Eqp_sat * Eqp_sat * Math::sigmoid(Eqp_sat); - f[7] = vd + psiqpp * (ONE + omega); - f[8] = vq - psidpp * (ONE + omega); - f[9] = telec - ((psidpp - id * Xdpp_) * iq - (psiqpp - iq * Xdpp_) * id); - f[10] = id - (ir * std::sin(delta) - ii * std::cos(delta)); - f[11] = iq - (ir * std::cos(delta) + ii * std::sin(delta)); - f[12] = ir + G_ * vr - B_ * vi - inr; - f[13] = ii + B_ * vr + G_ * vi - ini; + f[5] = psidpp - (psidp * Xd4_ + Eqp * Xd5_); + f[6] = ksat - SB_ * Math::qramp(Eqp - SA_); + f[7] = vd + psiqpp * (ONE + omega); + f[8] = vq - psidpp * (ONE + omega); + f[9] = telec - ((psidpp - id * Xdpp_) * iq - (psiqpp - iq * Xdpp_) * id); + f[10] = id - (ir * std::sin(delta) - ii * std::cos(delta)); + f[11] = iq - (ir * std::cos(delta) + ii * std::sin(delta)); + f[12] = ir + G_ * vr - B_ * vi - inr; + f[13] = ii + B_ * vr + G_ * vi - ini; /* 2 Gensal current source definitions */ f[14] = inr - (G_ * (std::sin(delta) * vd + std::cos(delta) * vq) - B_ * (-std::cos(delta) * vd + std::sin(delta) * vq)); diff --git a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/README.md b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/README.md index 8bbf0a206..ab2c0b4e1 100644 --- a/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/README.md +++ b/GridKit/Model/PhasorDynamics/SynchronousMachine/GENSAL/README.md @@ -74,7 +74,7 @@ $\psi''_q$ | [p.u.] | Total q-axis subtransient flux | Symbol | Units | Description | Note ------------|--------|-----------------------------------| ------ $\psi''_d$ | [p.u.] | Total d-axis subtransient flux | -$k_{sat}$ | [p.u.] | Additive saturation signal | +$k_{sat}$ | [p.u.] | Saturation factor | $V_d$ | [p.u.] | Machine internal voltage, d-axis | $V_q$ | [p.u.] | Machine internal voltage, q-axis | $T_e$ | [p.u.] | Electrical torque | @@ -108,7 +108,7 @@ $E_{fd}$ | [p.u.] | Field winding voltage from the excitation system | Ow \left( E_{fd}-E'_q-X_{d1} (I_d+X_{d3}(E'_q-\psi'_d-X_{d2}I_d)) - -k_{sat} + -E'_q k_{sat} \right)\\ \dot{\psi}'_d &= \dfrac{1}{T''_{d0}}(E'_q-\psi'_d-X_{d2}I_d)\\ \dot{\psi}''_q &= \dfrac{1}{T''_{q0}}(-\psi''_q-X_{q2}I_q) @@ -119,7 +119,7 @@ $E_{fd}$ | [p.u.] | Field winding voltage from the excitation system | Ow ``` math \begin{aligned} 0 &= -\psi''_d + E'_qX_{d5}+\psi'_dX_{d4}\\ - 0 &= -k_{sat} + S_B(E'_q-S_A)^2\sigma(E'_q-S_A)\\ + 0 &= -k_{sat} + S_B q(E'_q-S_A)\\ 0 &= -V_d -\psi''_q(1+\omega)\\ 0 &= -V_q +\psi''_d(1+\omega)\\ 0 &= -T_e +(\psi''_d-I_dX_d'')I_q-(\psi''_q-I_qX_d'')I_d\\ @@ -130,6 +130,9 @@ $E_{fd}$ | [p.u.] | Field winding voltage from the excitation system | Ow \end{aligned} ``` +CommonMath defines the primitive +[quadratic ramp](../../../../CommonMath.md#primitives) $q$. + ## Initialization Using the power-flow solution, initial currents are calculated from active and @@ -148,10 +151,10 @@ steady-state GENSAL equations. \psi''_d &= V_q\\ \psi'_d &= \psi''_d-(X_d''-X_\ell)I_d\\ E'_q &= \psi'_d+X_{d2}I_d\\ - k_{sat} &= S_B(E'_q-S_A)^2\sigma(E'_q-S_A)\\ + k_{sat} &= S_B q(E'_q-S_A)\\ T_e &= (\psi''_d-I_dX_d'')I_q-(\psi''_q-I_qX_d'')I_d\\ P_m &= T_e\\ - E_{fd} &= E'_q+X_{d1}(I_d+X_{d3}(E'_q-\psi'_d-X_{d2}I_d))+k_{sat} + E_{fd} &= E'_q+X_{d1}(I_d+X_{d3}(E'_q-\psi'_d-X_{d2}I_d))+E'_q k_{sat} \end{aligned} ``` diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index 1c60b16b0..3dd5a2e68 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -364,6 +364,13 @@ namespace GridKit gov->getSignals().template attachSignalNode(getSignal(speed)); } + if (govdata.signal_inputs.contains(Tgov1SignalInputs::pref)) + { + IdxT pref = govdata.signal_inputs.at(Tgov1SignalInputs::pref); + constexpr auto PREF = Tgov1ExternalVariables::PREF; + gov->getSignals().template attachSignalNode(getSignal(pref)); + } + if (govdata.signal_outputs.contains(Tgov1SignalOutputs::pmech)) { IdxT pmech = govdata.signal_outputs.at(Tgov1SignalOutputs::pmech); diff --git a/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp index 427648e1d..68a5def2f 100644 --- a/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp @@ -73,8 +73,6 @@ namespace GridKit success *= (bad_switch_model.verify() > 0); success *= invalidParameterCase(bus, Params::mva, 0.0); - success *= invalidParameterCase(bus, Params::Rpmax, 0.0); - success *= invalidParameterCase(bus, Params::Rqmin, 0.0); success *= invalidParameterCase(bus, Params::IL1, -0.1); success *= invalidParameterCase(bus, Params::VL1, 0.3); success *= invalidParameterCase(bus, Params::VA1, 0.3); diff --git a/tests/UnitTests/PhasorDynamics/ExciterEsdc1aTests.hpp b/tests/UnitTests/PhasorDynamics/ExciterEsdc1aTests.hpp index 881cd48ce..e701a7afb 100644 --- a/tests/UnitTests/PhasorDynamics/ExciterEsdc1aTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ExciterEsdc1aTests.hpp @@ -77,7 +77,6 @@ namespace GridKit success *= invalidParameterCase(Params::Ka, 0.0); success *= invalidParameterCase(Params::Ta, -0.1); success *= invalidParameterCase(Params::Te, -0.1); - success *= invalidParameterCase(Params::Tc, -0.1); success *= invalidParameterCase(Params::Tr, -0.1); success *= invalidParameterCase(Params::Tb, -0.1); success *= invalidParameterCase(Params::Tf1, -0.1); @@ -86,7 +85,7 @@ namespace GridKit success *= invalidParameterCase(Params::UEL, static_cast(2.0)); success *= invalidParameterCase(Params::UEL, static_cast(2.5)); success *= invalidParameterCase(Params::UEL, true); - success *= invalidParameterCase(Params::Se1, 0.0); + success *= invalidParameterCase(Params::Se1, -0.1); success *= invalidParameterCase(Params::E2, 2.8); success *= invalidParameterCase(Params::Se2, 0.08); success *= invalidParameterCase(Params::E1, -1.0); @@ -495,8 +494,7 @@ namespace GridKit return success.report(__func__); } - /// A fixed numerical answer key for all 11 ESDC1A equations. The - /// expected values are literals, not a second implementation of ESDC1A. + /// A fixed, analytically simple answer key for all 11 ESDC1A equations. TestOutcome residualEquations() { TestStatus success = true; @@ -508,20 +506,18 @@ namespace GridKit setAnswerKeyState(fixture.esdc1a); success *= (fixture.evaluate() == 0); - // Values are pinned after an independent one-time evaluation of the - // documented equations at setAnswerKeyState()/setAnswerKeyInputs(). const std::array(Internal::MAXIMUM)> expected{{ - {Internal::EFDP, 0.04000000000000004}, - {Internal::VC, 0.19442890089805262}, - {Internal::VR, 0.13666666666666663}, - {Internal::VF, -0.022222222222222213}, - {Internal::XLL, 0.024999999999999994}, - {Internal::EV, -0.27999999999999986}, - {Internal::VLL, -0.007500000000000031}, - {Internal::VHV, 0.31535073999664665}, - {Internal::SE, -0.07541019662496842}, - {Internal::VFE, 0.1600000000000001}, - {Internal::EFD, 0.9100000000000001}, + {Internal::EFDP, 0.04}, + {Internal::VC, 0.27}, + {Internal::VR, 0.97}, + {Internal::VF, -0.02}, + {Internal::XLL, 0.025}, + {Internal::EV, -0.26}, + {Internal::VLL, -0.0075}, + {Internal::VHV, 0.3}, + {Internal::SE, 0.16}, + {Internal::VFE, 0.14}, + {Internal::EFD, 1.32}, }}; success *= (static_cast(fixture.esdc1a.getResidual().getSize()) == expected.size()); @@ -553,7 +549,7 @@ namespace GridKit setDerivative(fixture.esdc1a, {{Internal::EFDP, 0.1}, {Internal::VF, 0.05}}); success *= (fixture.evaluate() == 0); success *= residualsMatch(fixture.esdc1a, - {{Internal::EFDP, 0.7}, {Internal::VF, -0.13571428571428573}}, + {{Internal::EFDP, 0.7}, {Internal::VF, -19.0 / 140.0}}, "field-voltage and feedback drive"); // Summing junction: UEL < 2 excludes the UEL input from the error. @@ -631,15 +627,14 @@ namespace GridKit return success.report(__func__); } - /// High-value gate selection, quadratic saturation, field-voltage-state + /// High-value gate selection, scaled-quadratic saturation, field-voltage-state /// limiting, and the speed multiplier at driven states with literal /// expectations. TestOutcome excitationLimits() { TestStatus success = true; - // The gate passes the larger of VLL and VUEL when UEL < 2; the - // smooth maximum keeps two-sided sensitivity at the tie. + // The gate passes the larger of VLL and VUEL when UEL < 2. struct GateCase { const char* label; @@ -647,10 +642,9 @@ namespace GridKit RealT expected; }; - const std::array gate_cases{{ + const std::array gate_cases{{ {"gate selects the lead-lag branch", -0.5, 0.3}, - {"gate selects the UEL branch", 0.8, 0.6000000000000001}, - {"gate tie point", 0.5, 0.30288811325233306}, + {"gate selects the UEL branch", 0.8, 0.6}, }}; Fixture gate(makeData()); @@ -664,8 +658,8 @@ namespace GridKit success *= residualsMatch(gate.esdc1a, {{Internal::VHV, test_case.expected}}, test_case.label); } - // Both valid point orderings produce the same quadratic curve at the - // supplied points and on either side of the fitted knee. + // Both valid point orderings recover E S_E(E) = 0.25(E - 1)^2, + // including the supplied points, an interior point, and the knee. struct SaturationOrderCase { const char* label; @@ -676,10 +670,17 @@ namespace GridKit }; const std::array saturation_order_cases{{ - {"ascending saturation points", 2.4, 0.1, 3.2, 0.5}, - {"descending saturation points", 3.2, 0.5, 2.4, 0.1}, + {"ascending saturation points", 2.0, 0.125, 4.0, 0.5625}, + {"descending saturation points", 4.0, 0.5625, 2.0, 0.125}, }}; + struct SaturationEvaluation + { + RealT efdp; + RealT se; + RealT expected; + }; + for (const auto& test_case : saturation_order_cases) { auto data = makeResidualData(); @@ -691,29 +692,23 @@ namespace GridKit saturation.attachAllInputs(); success *= saturation.initialize(1.2); - setState(saturation.esdc1a, {{Internal::EFDP, 2.4}, {Internal::SE, 0.0}}); - success *= (saturation.evaluate() == 0); - success *= residualsMatch(saturation.esdc1a, - {{Internal::SE, 0.1}}, - test_case.label); - - setState(saturation.esdc1a, {{Internal::EFDP, 3.2}}); - success *= (saturation.evaluate() == 0); - success *= residualsMatch(saturation.esdc1a, - {{Internal::SE, 0.5}}, - test_case.label); - - setState(saturation.esdc1a, {{Internal::EFDP, 2.0}, {Internal::SE, 0.05}}); - success *= (saturation.evaluate() == 0); - success *= residualsMatch(saturation.esdc1a, - {{Internal::SE, -0.035410196624968436}}, - test_case.label); + const std::array evaluations{{ + {test_case.e1, 0.0, test_case.e1 * test_case.se1}, + {test_case.e2, 0.0, test_case.e2 * test_case.se2}, + {3.0, 0.25, 0.75}, + {1.0, 0.05, -0.05}, + }}; - setState(saturation.esdc1a, {{Internal::EFDP, 1.0}}); - success *= (saturation.evaluate() == 0); - success *= residualsMatch(saturation.esdc1a, - {{Internal::SE, -0.05}}, - test_case.label); + for (const auto& evaluation : evaluations) + { + setState(saturation.esdc1a, + {{Internal::EFDP, evaluation.efdp}, + {Internal::SE, evaluation.se}}); + success *= (saturation.evaluate() == 0); + success *= residualsMatch(saturation.esdc1a, + {{Internal::SE, evaluation.expected}}, + test_case.label); + } } auto disabled_data = makeResidualData(); @@ -814,7 +809,7 @@ namespace GridKit // enabled. for (const auto& [enabled, expected] : std::array, 2>{{ {false, 0.0}, - {true, 0.06000000000000005}, + {true, 0.06}, }}) { auto data = makeData(); @@ -988,8 +983,10 @@ namespace GridKit PhasorDynamics::Exciter::Esdc1a esdc1a; }; - static constexpr RealT kStateVr = 0.9; - static constexpr RealT kStateVi = 0.4; + static constexpr RealT kStateVr = 0.8; + static constexpr RealT kStateVi = 0.6; + // Keeps both smooth-maximum inputs active in the Jacobian comparison. + static constexpr RealT kJacobianVuel = 0.334; Data makeMinimalData() const { @@ -1088,36 +1085,34 @@ namespace GridKit template void setAnswerKeyInputs(Fixture& fixture) const { - fixture.input(External::OMEGA) = 0.03; + fixture.input(External::OMEGA) = 0.05; fixture.input(External::VREF) = 1.05; fixture.input(External::VS) = 0.04; - fixture.input(External::VUEL) = 0.334; + fixture.input(External::VUEL) = 0.2; } /// The rich state shared by the residual answer key and the Jacobian - /// comparison. Every row is distinct so a swapped index cannot pass, - /// and VLL sits close enough to VUEL that the smooth gate keeps - /// two-sided sensitivity. + /// comparison. Every residual row is distinct so a swapped index cannot pass. template void setAnswerKeyState(PhasorDynamics::Exciter::Esdc1a& esdc1a) const { setState(esdc1a, - {{Internal::EFDP, 2.00}, + {{Internal::EFDP, 2.4}, {Internal::VC, 0.95}, {Internal::VR, 0.45}, - {Internal::VF, 0.06}, + {Internal::VF, 0.04}, {Internal::XLL, 0.30}, {Internal::EV, 0.36}, {Internal::VLL, 0.33}, - {Internal::VHV, 0.02}, - {Internal::SE, 0.09}, + {Internal::VHV, 0.03}, + {Internal::SE, 0.08}, {Internal::VFE, 0.42}, - {Internal::EFD, 1.15}}); + {Internal::EFD, 1.2}}); setDerivative(esdc1a, {{Internal::EFDP, 0.01}, {Internal::VC, -0.02}, {Internal::VR, 0.03}, - {Internal::VF, -0.04}, + {Internal::VF, -0.02}, {Internal::XLL, 0.05}}); } @@ -1482,6 +1477,7 @@ namespace GridKit fixture.attachAllInputs(); success *= fixture.initialize(1.2); setAnswerKeyInputs(fixture); + fixture.input(External::VUEL) = kJacobianVuel; setAnswerKeyState(fixture.esdc1a); numberVariables(fixture); success *= (fixture.evaluate() == 0); @@ -1510,6 +1506,7 @@ namespace GridKit } setAnswerKeyInputs(fixture); + fixture.input(External::VUEL) = kJacobianVuel; setAnswerKeyState(fixture.esdc1a); fixture.esdc1a.updateTime(0.0, 1.0); success *= (fixture.evaluate() == 0); diff --git a/tests/UnitTests/PhasorDynamics/ExciterIeeet1Tests.hpp b/tests/UnitTests/PhasorDynamics/ExciterIeeet1Tests.hpp index 57d49660b..8f40ab78b 100644 --- a/tests/UnitTests/PhasorDynamics/ExciterIeeet1Tests.hpp +++ b/tests/UnitTests/PhasorDynamics/ExciterIeeet1Tests.hpp @@ -135,7 +135,7 @@ namespace GridKit for (IdxT i = 0; i < exciter.getResidual().getSize(); ++i) { success *= std::isfinite(f[i]); - success *= isEqual(f[i], static_cast(0.0)); + success *= isEqual(f[i], static_cast(0.0), 100.0 * std::numeric_limits::epsilon()); } y[2] = 4.0; @@ -161,13 +161,13 @@ namespace GridKit y[4] = 0.02; exciter.y().setDataUpdated(); exciter.evaluateResidual(); - success *= isEqual(f[1], static_cast(1.0e3)); + success *= isEqual(f[1], static_cast(900.0)); y[4] = 0.0; y[1] = 1.0; exciter.y().setDataUpdated(); exciter.evaluateResidual(); - success *= isEqual(f[2], static_cast(1.0e3)); + success *= isEqual(f[2], static_cast(900.0)); y[1] = 0.0; y[5] = 1.0; @@ -190,7 +190,7 @@ namespace GridKit using Params = PhasorDynamics::Exciter::Ieeet1Parameters; auto data = makeTestData(); - data.parameters[Params::Se1] = 0.0; + data.parameters[Params::Se1] = -0.1; PhasorDynamics::Bus bus(3.0, 4.0); PhasorDynamics::Exciter::Ieeet1 exciter(&bus, data); diff --git a/tests/UnitTests/PhasorDynamics/GensalTests.hpp b/tests/UnitTests/PhasorDynamics/GensalTests.hpp index 10e84e0c7..ae15c4e9f 100644 --- a/tests/UnitTests/PhasorDynamics/GensalTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GensalTests.hpp @@ -267,7 +267,7 @@ namespace GridKit const std::vector res_answer = { 0.0, 0.0, - 2.1083333333333334, + 2.2083333333333335, -1.028125, 0.65, 0.0, diff --git a/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp b/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp index 473d333cd..532a21d81 100644 --- a/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp +++ b/tests/UnitTests/PhasorDynamics/GovernorTgov1Tests.hpp @@ -101,9 +101,9 @@ namespace GridKit PhasorDynamics::Governor::Tgov1 gov(&pmech, &omega); // Test answer keys - const std::vector res_answer = {0.0, - -1.0, - -0.2}; + const std::vector res_answer = {static_cast(2.0) / static_cast(15.0), + -2.0, + static_cast(1.0) / static_cast(3.0)}; bus.allocate(); gen.allocate(); diff --git a/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp b/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp index 5c8386b6f..234c13455 100644 --- a/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp +++ b/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp @@ -418,7 +418,10 @@ namespace GridKit PhasorDynamics::SystemModel* system = new PhasorDynamics::SystemModel(); + PhasorDynamics::SignalNode pmech; PhasorDynamics::Governor::Tgov1 tgov1; + tgov1.getSignals() + .template assignSignalNode(&pmech); system->addComponent(&tgov1); success *= system->allocate() == 0; diff --git a/tests/UnitTests/Utilities/CaseFormatTests.hpp b/tests/UnitTests/Utilities/CaseFormatTests.hpp index 74d60227f..4cbf66b08 100644 --- a/tests/UnitTests/Utilities/CaseFormatTests.hpp +++ b/tests/UnitTests/Utilities/CaseFormatTests.hpp @@ -264,7 +264,7 @@ namespace GridKit "devices": [ { "class": "Branch", "ports": {"bus1":1, "bus2":2}, "id": "BR1", "params": {"R":0.0, "X":0.1, "G":0.0, "B":0.0, "tap":1.05, "phase":0.1} }, { "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": {"speed": 1, "pmech":2}, "id": "DV2", "params": {"R":0.05, "T1":0.5,"T2":2.5, "T3":7.5, "Pvmax":0.0, "Pvmin":1.0, "Dt":0.0}}, + { "class": "Tgov1", "ports": {"speed": 1, "pref":10, "pmech":2}, "id": "DV2", "params": {"R":0.05, "T1":0.5,"T2":2.5, "T3":7.5, "Pvmax":1.0, "Pvmin":0.0, "Dt":0.0}}, { "class": "Esdc1a", "ports": {"bus":1, "speed":1, "vref":4, "vs":5, "vuel":6, "efd":3}, "id": "DV5", "params": {"Tr":0.0, "Ka":40.0, "Ta":0.1, "Tb":0.0, "Tc":0.0, "Vrmax":1.0, "Vrmin":-1.0, "Ke":0.1, "Te":0.5, "Kf":0.05, "Tf1":0.7, "Spdmlt":false, "E1":2.8, "Se1":0.08, "E2":3.7, "Se2":0.33, "UEL":0, "exclim":true}, "mon": ["efd", "vc", "vr", "vf", "se", "vfe"] }, { "class": "Hygov", "ports": {"speed": 1, "pmech": 7, "pref": 8, "paux": 9}, "id": "DV6", "params": {"Trate": 80.0, "Rperm": 0.05, "Rtemp": 0.35, "Tr": 5.0, "Tf": 0.05, "Tg": 0.5, "Velm": 0.2, "Gmax": 0.98, "Gmin": 0.02, "Tw": 1.2, "At": 1.1, "Dturb": 0.4, "Qnl": 0.08, "Tn": 0.7, "Tnp": 1.4, "db1": 0.01, "db2": 0.02, "Hdam": 1.05, @@ -386,10 +386,11 @@ namespace GridKit success *= std::get(result.gov[0].parameters[Governor::Tgov1Parameters::T1]) == 0.5; success *= std::get(result.gov[0].parameters[Governor::Tgov1Parameters::T2]) == 2.5; success *= std::get(result.gov[0].parameters[Governor::Tgov1Parameters::T3]) == 7.5; - success *= std::get(result.gov[0].parameters[Governor::Tgov1Parameters::Pvmax]) == 0; - success *= std::get(result.gov[0].parameters[Governor::Tgov1Parameters::Pvmin]) == 1; + success *= std::get(result.gov[0].parameters[Governor::Tgov1Parameters::Pvmax]) == 1; + success *= std::get(result.gov[0].parameters[Governor::Tgov1Parameters::Pvmin]) == 0; success *= std::get(result.gov[0].parameters[Governor::Tgov1Parameters::Dt]) == 0; success *= result.gov[0].signal_inputs[Governor::Tgov1SignalInputs::speed] == 1; + success *= result.gov[0].signal_inputs[Governor::Tgov1SignalInputs::pref] == 10; success *= result.gov[0].signal_outputs[Governor::Tgov1SignalOutputs::pmech] == 2; success *= result.gov[0].disambiguation_string == "DV2"; From f18cc720a8258e117f5ff9a929fafbb4905bcc6f Mon Sep 17 00:00:00 2001 From: lukelowry Date: Thu, 13 Aug 2026 19:42:20 -0500 Subject: [PATCH 02/10] Effective saturation parameter --- .../PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp | 3 ++- .../Exciter/ESDC1A/Esdc1aData.hpp | 2 +- .../Exciter/ESDC1A/Esdc1aImpl.hpp | 6 ++++-- .../PhasorDynamics/Exciter/ESDC1A/README.md | 19 +++++++++---------- .../PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp | 2 ++ .../Exciter/IEEET1/Ieeet1Impl.hpp | 9 ++++++--- .../PhasorDynamics/Exciter/IEEET1/README.md | 6 +++--- .../PhasorDynamics/Governor/Tgov1/README.md | 2 +- 8 files changed, 28 insertions(+), 21 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp index b7f791886..5d9e8f36d 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp +++ b/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$E_{\mathrm{fd}}'S_E(E_{\mathrm{fd}}')\f$ Scaled-quadratic saturation contribution [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}; diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aData.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aData.hpp index 3de51f210..d0bfa5d7c 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aData.hpp +++ b/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$E_{\mathrm{fd}}'S_E(E_{\mathrm{fd}}')\f$ Scaled-quadratic saturation contribution [p.u.] + se, ///< \f$s_e\f$ Scaled-quadratic saturation contribution [p.u.] vfe ///< \f$V_{\mathrm{FE}}\f$ Exciter feedback drive [p.u.] }; diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp index 2a0e43961..099d673d8 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp @@ -351,7 +351,7 @@ namespace GridKit const ScalarT xll0 = ev0; const ScalarT vref0 = ev0 + vc0 - vs0 - uel_on_ * vuel0; - Ke_ = ke0; + Ke_eff_ = ke0; y[EFDP] = efdp0; y[VC] = vc0; @@ -564,7 +564,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_ * efdp + se; + 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); @@ -762,6 +762,8 @@ namespace GridKit template void Esdc1a::setDerivedParameters() { + Ke_eff_ = Ke_; + // The lags are raised to the floor in place, so a negative value is // rejected here while the value as read is still available. verify() // reports the count. diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md index 9840a4646..164ada637 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md @@ -166,7 +166,7 @@ Symbol | Units | Description $e_V$ | [p.u.] | Voltage-error summing output | $V_{\mathrm{LL}}$ | [p.u.] | Input lead-lag output | $V_{\mathrm{HV}}$ | [p.u.] | High-value gate output | -$E_{\mathrm{fd}}'S_E(E_{\mathrm{fd}}')$ | [p.u.] | Scaled-quadratic saturation contribution | +$s_e$ | [p.u.] | Scaled-quadratic saturation contribution | $E_{\mathrm{fd}}'S_E(E_{\mathrm{fd}}')$ $V_{\mathrm{FE}}$ | [p.u.] | Exciter feedback drive | $E_{\mathrm{fd}}$ | [p.u.] | Field-voltage output | Published through `efd` @@ -263,12 +263,12 @@ of [Appendix A](#appendix-a-awmin). & s_{\mathrm{UEL}} = 1 \end{cases} \\ 0 &= - -E_{\mathrm{fd}}'S_E(E_{\mathrm{fd}}') + -s_e + S_B q\left(E_{\mathrm{fd}}' - S_A\right) \\ 0 &= -V_{\mathrm{FE}} - + K_E E_{\mathrm{fd}}' - + E_{\mathrm{fd}}'S_E(E_{\mathrm{fd}}') \\ + + K_E^{\mathrm{eff}} E_{\mathrm{fd}}' + + s_e \\ 0 &= -E_{\mathrm{fd}} + \left(1 + s_{\mathrm{spd}}\omega\right)E_{\mathrm{fd}}' @@ -315,17 +315,16 @@ routed through the gate: E_{\mathrm{fd}}' &\leftarrow \dfrac{E_{\mathrm{fd}}}{1 + s_{\mathrm{spd}}\omega} \\ - E_{\mathrm{fd}}'S_E(E_{\mathrm{fd}}') + s_e &\leftarrow S_B q\left(E_{\mathrm{fd}}' - S_A\right) \\ - K_E + K_E^{\mathrm{eff}} &\leftarrow \begin{cases} - \dfrac{V_R^{\max}/10-E_{\mathrm{fd}}'S_E(E_{\mathrm{fd}}')}{E_{\mathrm{fd}}'} & K_E=0 \\ + \dfrac{V_R^{\max}/10-s_e}{E_{\mathrm{fd}}'} & K_E=0 \\ K_E & K_E\ne 0 \end{cases} \\ V_{\mathrm{FE}} - &\leftarrow K_E E_{\mathrm{fd}}' - + E_{\mathrm{fd}}'S_E(E_{\mathrm{fd}}') \\ + &\leftarrow K_E^{\mathrm{eff}} E_{\mathrm{fd}}' + s_e \\ V_R &\leftarrow V_{\mathrm{FE}} \\ V_{\mathrm{HV}} @@ -386,7 +385,7 @@ Output | Units | Description | Note `vc` | [p.u.] | Filtered terminal-voltage magnitude | $V_C$ `vr` | [p.u.] | Voltage-regulator output | $V_R$ `vf` | [p.u.] | Stabilizing feedback state | $V_F$ -`se` | [p.u.] | Scaled-quadratic saturation contribution | $E_{\mathrm{fd}}'S_E(E_{\mathrm{fd}}')$ +`se` | [p.u.] | Scaled-quadratic saturation contribution | $s_e$ `vfe` | [p.u.] | Exciter feedback drive | $V_{\mathrm{FE}}$ ## Testing diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp index de663f014..783483e82 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp @@ -146,6 +146,8 @@ namespace GridKit RealT Ispdlim_{0.0}; ///< Speed limit flag indicator // Model Derived parameters + RealT Ke_eff_{Ke_}; + // Saturation coefficients derived from E1, E2, Se1, and Se2. RealT SA_{0}; RealT SB_{0}; diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp index a7eb0eb4b..e994ea579 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp @@ -233,13 +233,14 @@ namespace GridKit ScalarT efdp = efd0 / (ONE + omega * Ispdlim_); ScalarT ksat = SB_ * Math::qramp(efdp - SA_); + Ke_eff_ = Ke_; if (Ke_ == ZERO) { - Ke_ = (Vrmax_ / 10.0 - static_cast(ksat)) / static_cast(efdp); + Ke_eff_ = (Vrmax_ / 10.0 - static_cast(ksat)) / static_cast(efdp); } ScalarT ve = ksat; - ScalarT vr = Ke_ * efdp + ve; + ScalarT vr = Ke_eff_ * efdp + ve; ScalarT vtr = vr / Ka_; ScalarT vf{0}; ScalarT vfx = (Kf_ / Tf_) * efdp; @@ -351,7 +352,7 @@ namespace GridKit // Internal Differential Equations f[0] = -vts_dot + (Ec - vts) / Tr_; f[1] = -vr_dot + Math::antiwindup(vr, func, Vrmin_, Vrmax_); - f[2] = -efdp_dot + (vr - ve - Ke_ * efdp) / Te_; + f[2] = -efdp_dot + (vr - ve - Ke_eff_ * efdp) / Te_; f[3] = -vfx_dot + vf / Tf_; // Internal Algebraic Equations @@ -467,6 +468,8 @@ namespace GridKit Ispdlim_ = std::get(data.parameters.at(Parameter::Ispdlim)); } + Ke_eff_ = Ke_; + Tr_ = std::max(Tr_, TIME_CONSTANT_MINIMUM); Ta_ = std::max(Ta_, TIME_CONSTANT_MINIMUM); Te_ = std::max(Te_, TIME_CONSTANT_MINIMUM); diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md index 4f689e567..ff89ac45f 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md @@ -161,7 +161,7 @@ The IEEET1 differential equations, as derived from the model diagram, are: 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 E_{fd}' + \dfrac{1}{T_E}\left(V_R - V_E - K_E^{\mathrm{eff}} E_{fd}'\right) \\ 0 &= -\dot V_{fx} + \dfrac{1}{T_F}\left(V_f\right) \end{aligned} ``` @@ -201,13 +201,13 @@ with the current input values. 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) \\ - K_E &\leftarrow + K_E^{\mathrm{eff}} &\leftarrow \begin{cases} \dfrac{V_R^{\max}/10-k_\text{sat}}{E_{fd}'} & K_E=0 \\ K_E & K_E\ne 0 \end{cases} \\ V_E &= k_\text{sat} \\ - V_R &= K_E\, E_{fd}' + V_E \\ + V_R &= K_E^{\mathrm{eff}} E_{fd}' + V_E \\ V_{tr} &= \dfrac{V_R}{K_A} \\ V_{fx} &= \dfrac{K_F}{T_F}\, E_{fd}' \\ V_{ts} &= E_{C,0} \\ diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/README.md b/GridKit/Model/PhasorDynamics/Governor/Tgov1/README.md index 0315ffe2b..f18a5d2e5 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/README.md +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/README.md @@ -102,7 +102,7 @@ The TGOV1 differential equations, as derived from the model diagram, are \end{aligned} ``` -CommonMath defines the [anti-windup](../../../../CommonMath.md#antiwindup) +CommonMath defines the [Antiwindup](../../../../CommonMath.md#antiwindup) target and smooth approximation. ### Algebraic Equations From 27d2ce8888929b028bec60d0219abdcd24a7268f Mon Sep 17 00:00:00 2001 From: lukelowry Date: Thu, 13 Aug 2026 21:05:52 -0500 Subject: [PATCH 03/10] derived parameter ke --- .../PhasorDynamics/Exciter/ESDC1A/README.md | 17 +++++++++++------ .../PhasorDynamics/Exciter/IEEET1/README.md | 16 +++++++++++----- 2 files changed, 22 insertions(+), 11 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md index 164ada637..b2102fa31 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md @@ -127,6 +127,17 @@ and when both saturation values are positive, \end{aligned} ``` +During initialization, the effective +exciter constant is derived from the initial operating point: + +```math +K_E^{\mathrm{eff}} = + \begin{cases} + \dfrac{V_R^{\max}/10-s_e}{E_{\mathrm{fd}}'} & K_E=0 \\ + K_E & K_E\ne 0. + \end{cases} +``` + ## Model Ports Name | Port | Init | Description @@ -317,12 +328,6 @@ routed through the gate: \dfrac{E_{\mathrm{fd}}}{1 + s_{\mathrm{spd}}\omega} \\ s_e &\leftarrow S_B q\left(E_{\mathrm{fd}}' - S_A\right) \\ - K_E^{\mathrm{eff}} - &\leftarrow - \begin{cases} - \dfrac{V_R^{\max}/10-s_e}{E_{\mathrm{fd}}'} & K_E=0 \\ - K_E & K_E\ne 0 - \end{cases} \\ V_{\mathrm{FE}} &\leftarrow K_E^{\mathrm{eff}} E_{\mathrm{fd}}' + s_e \\ V_R diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md index ff89ac45f..ed1a44ac5 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md @@ -102,6 +102,17 @@ When both saturation values are positive, the non-extraneous solution is: \end{aligned} ``` +During initialization, the effective +exciter constant is derived from the initial operating point: + +```math +K_E^{\mathrm{eff}} = + \begin{cases} + \dfrac{V_R^{\max}/10-k_\text{sat}}{E_{fd}'} & K_E=0 \\ + K_E & K_E\ne 0. + \end{cases} +``` + ## Model Variables ### Internal Variables @@ -201,11 +212,6 @@ with the current input values. 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) \\ - K_E^{\mathrm{eff}} &\leftarrow - \begin{cases} - \dfrac{V_R^{\max}/10-k_\text{sat}}{E_{fd}'} & K_E=0 \\ - K_E & K_E\ne 0 - \end{cases} \\ V_E &= k_\text{sat} \\ V_R &= K_E^{\mathrm{eff}} E_{fd}' + V_E \\ V_{tr} &= \dfrac{V_R}{K_A} \\ From d7fffaba09225cdaedb9bcf5b3efa68e8ede05ea Mon Sep 17 00:00:00 2001 From: lukelowry Date: Mon, 17 Aug 2026 10:01:50 -0500 Subject: [PATCH 04/10] claarify exciter slope parameter --- .../PhasorDynamics/Exciter/ESDC1A/README.md | 31 ++++++++++---- .../PhasorDynamics/Exciter/IEEET1/README.md | 41 +++++++++++++------ 2 files changed, 51 insertions(+), 21 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md index b2102fa31..502d5724d 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md @@ -30,7 +30,7 @@ $T_B$ | [sec] | `Tb` | Input lead-lag den $T_C$ | [sec] | `Tc` | Input lead-lag numerator time constant | 0.0 $V_R^{\max}$ | [p.u.] | `Vrmax` | Maximum voltage-regulator output | 1.0 $V_R^{\min}$ | [p.u.] | `Vrmin` | Minimum voltage-regulator output | -1.0 -$K_E$ | [p.u.] | `Ke` | Exciter constant | 0.1 +$K_E$ | [p.u.] | `Ke` | Exciter field resistance line slope margin | 0.1 $T_E$ | [sec] | `Te` | Exciter time constant | 0.5 $K_F$ | [p.u.] | `Kf` | Stabilizing feedback gain | 0.05 $T_{F1}$ | [sec] | `Tf1` | Stabilizing feedback time constant | 0.7 @@ -127,17 +127,32 @@ and when both saturation values are positive, \end{aligned} ``` -During initialization, the effective -exciter constant is derived from the initial operating point: +$K_E$ is the configured exciter field resistance line slope margin. A nonzero configured value is used directly. A configured $K_E=0$ selects the +[PSS/E-compatible model](https://www.powerworld.com/WebHelp/Content/MainDocumentation_HTML/Transient_Stability_Dialog_Options_Power_System_Model.htm). +At initialization, requiring $\dot E_{\mathrm{fd}}'=0$ and +$V_R=V_R^{\max}/10$ gives ```math -K_E^{\mathrm{eff}} = - \begin{cases} - \dfrac{V_R^{\max}/10-s_e}{E_{\mathrm{fd}}'} & K_E=0 \\ - K_E & K_E\ne 0. - \end{cases} +\begin{aligned} +0 + &= V_R + -K_E^{\mathrm{eff}}E_{\mathrm{fd}}' + -s_e, \\ +K_E^{\mathrm{eff}} + &= + \begin{cases} + \dfrac{1}{E_{\mathrm{fd}}'} + \left(\dfrac{V_R^{\max}}{10}-s_e\right) + & K_E=0, \\ + K_E + & K_E\ne 0. + \end{cases} +\end{aligned} ``` +Thus $K_E^{\mathrm{eff}}$ is the resolved value of the same exciter +coefficient, not an additional model input. + ## Model Ports Name | Port | Init | Description diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md index ed1a44ac5..7384d6b04 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md @@ -19,7 +19,7 @@ Symbol | Units | Description | Typical Value | No $T_R$ | [sec] | Time constant for voltage sensing | 0 | $K_A$ | [p.u.] | Coefficient for voltage regulation | 50 | $T_A$ | [sec] | Time constant for voltage regulation | 0.04 | -$K_E$ | [p.u.] | Coefficient for excitation system | -0.06 | +$K_E$ | [p.u.] | Exciter field resistance line slope margin | -0.06 | $T_E$ | [sec] | Time constant for excitation system | 0.6 | $K_F$ | [p.u.] | Coefficient for feedback | 0.09 | $T_F$ | [sec] | Time constant for feedback | 1.46 | @@ -102,17 +102,32 @@ When both saturation values are positive, the non-extraneous solution is: \end{aligned} ``` -During initialization, the effective -exciter constant is derived from the initial operating point: +$K_E$ is the configured exciter field resistance line slope margin. A nonzero configured value is used directly. A configured $K_E=0$ selects the +[PSS/E-compatible model](https://www.powerworld.com/WebHelp/Content/MainDocumentation_HTML/Transient_Stability_Dialog_Options_Power_System_Model.htm). +At initialization, requiring $\dot E_{fd}'=0$ and +$V_R=V_R^{\max}/10$ gives ```math -K_E^{\mathrm{eff}} = - \begin{cases} - \dfrac{V_R^{\max}/10-k_\text{sat}}{E_{fd}'} & K_E=0 \\ - K_E & K_E\ne 0. - \end{cases} +\begin{aligned} +0 + &= V_R + -k_\text{sat} + -K_E^{\mathrm{eff}}E_{fd}', \\ +K_E^{\mathrm{eff}} + &= + \begin{cases} + \dfrac{1}{E_{fd}'} + \left(\dfrac{V_R^{\max}}{10}-k_\text{sat}\right) + & K_E=0, \\ + K_E + & K_E\ne 0. + \end{cases} +\end{aligned} ``` +Thus $K_E^{\mathrm{eff}}$ is the resolved value of the same exciter +coefficient, not an additional model input. + ## Model Variables ### Internal Variables @@ -200,7 +215,7 @@ Here $q$ is GridKit's [Quadratic Ramp](../../../../CommonMath.md#primitives). The implementation first applies $T \leftarrow \max(T, 10^{-3})$ for $T \in \{T_R, T_A, T_E, T_F\}$. This should be replaced with a structural template change in the future. The machine initializes $E_{fd}$ first. IEEET1 -reads that value as $E_{fd,0}$, along with any attached $\omega$ and $V_S$, and +reads that value, 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 is included when enabled, and the speed-limit @@ -209,16 +224,16 @@ with the current input values. ```math \begin{aligned} - E_{C,0} &:= \sqrt{V_r^2 + V_i^2} \\ - E_{fd}' &= \dfrac{E_{fd,0}}{1 + I_{\mathrm{spdlim}}\,\omega} \\ + E_C &:= \sqrt{V_r^2 + V_i^2} \\ + E_{fd}' &= \dfrac{E_{fd}}{1 + I_{\mathrm{spdlim}}\,\omega} \\ k_\text{sat} &= S_B\, q(E_{fd}' - S_A) \\ V_E &= k_\text{sat} \\ V_R &= K_E^{\mathrm{eff}} E_{fd}' + V_E \\ V_{tr} &= \dfrac{V_R}{K_A} \\ V_{fx} &= \dfrac{K_F}{T_F}\, E_{fd}' \\ - V_{ts} &= E_{C,0} \\ + V_{ts} &= E_C \\ V_f &= 0 \\ - V_\text{ref} &= E_{C,0} + V_{tr} - V_{UEL} - V_{OEL} - V_S + V_\text{ref} &= E_C + V_{tr} - V_{UEL} - V_{OEL} - V_S \end{aligned} ``` From 776322ca67b3739af50936cbafd9392f6fd57538 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Mon, 17 Aug 2026 10:52:18 -0500 Subject: [PATCH 05/10] justification and further references --- .../PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp | 10 ++++++++++ .../PhasorDynamics/Exciter/ESDC1A/README.md | 16 +++++++++++----- .../PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp | 10 ++++++++++ .../PhasorDynamics/Exciter/IEEET1/README.md | 16 +++++++++++----- 4 files changed, 42 insertions(+), 10 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp index 099d673d8..5b8711916 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp @@ -222,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.4-2014 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 diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md index 502d5724d..75bf163e8 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md @@ -30,7 +30,7 @@ $T_B$ | [sec] | `Tb` | Input lead-lag den $T_C$ | [sec] | `Tc` | Input lead-lag numerator time constant | 0.0 $V_R^{\max}$ | [p.u.] | `Vrmax` | Maximum voltage-regulator output | 1.0 $V_R^{\min}$ | [p.u.] | `Vrmin` | Minimum voltage-regulator output | -1.0 -$K_E$ | [p.u.] | `Ke` | Exciter field resistance line slope margin | 0.1 +$K_E$ | [p.u.] | `Ke` | Exciter field resistance line slope margin; 0 requests automatic calculation, not a zero coefficient | 0.1 $T_E$ | [sec] | `Te` | Exciter time constant | 0.5 $K_F$ | [p.u.] | `Kf` | Stabilizing feedback gain | 0.05 $T_{F1}$ | [sec] | `Tf1` | Stabilizing feedback time constant | 0.7 @@ -127,10 +127,16 @@ and when both saturation values are positive, \end{aligned} ``` -$K_E$ is the configured exciter field resistance line slope margin. A nonzero configured value is used directly. A configured $K_E=0$ selects the -[PSS/E-compatible model](https://www.powerworld.com/WebHelp/Content/MainDocumentation_HTML/Transient_Stability_Dialog_Options_Power_System_Model.htm). -At initialization, requiring $\dot E_{\mathrm{fd}}'=0$ and -$V_R=V_R^{\max}/10$ gives +$K_E$ is the configured exciter field resistance line slope margin. A nonzero configured value is used directly. The +[IEEE Std 421.4-2014](https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=6803835) states: + +> In some programs, if $K_E$ is entered as zero, $K_E$ is automatically calculated by the program to represent a self-excited shunt field and a trimmed rheostat as its initial condition. + +GridKit implements the +[PSS/E-compatible automatic-parameter rule](https://www.powerworld.com/WebHelp/Content/MainDocumentation_HTML/Transient_Stability_Dialog_Options_Power_System_Model.htm). +The divisor $10$ in $V_R=V_R^{\max}/10=0.1V_R^{\max}$ is unitless and +sets $V_R$ to 10% of the maximum regulator output, retaining the per-unit +units of $V_R^{\max}$: ```math \begin{aligned} diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp index e994ea579..570795a7d 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp @@ -189,6 +189,16 @@ namespace GridKit * - Attached external signals (omega, V_S) * * Enabled saturation is included via ksat computed from efdp and SA, SB. + * + * @warning IEEE Std 421.4-2014 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. */ template int Ieeet1::initialize() diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md index 7384d6b04..0ddbede74 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md @@ -19,7 +19,7 @@ Symbol | Units | Description | Typical Value | No $T_R$ | [sec] | Time constant for voltage sensing | 0 | $K_A$ | [p.u.] | Coefficient for voltage regulation | 50 | $T_A$ | [sec] | Time constant for voltage regulation | 0.04 | -$K_E$ | [p.u.] | Exciter field resistance line slope margin | -0.06 | +$K_E$ | [p.u.] | Exciter field resistance line slope margin; 0 requests automatic calculation, not a zero coefficient | -0.06 | $T_E$ | [sec] | Time constant for excitation system | 0.6 | $K_F$ | [p.u.] | Coefficient for feedback | 0.09 | $T_F$ | [sec] | Time constant for feedback | 1.46 | @@ -102,10 +102,16 @@ When both saturation values are positive, the non-extraneous solution is: \end{aligned} ``` -$K_E$ is the configured exciter field resistance line slope margin. A nonzero configured value is used directly. A configured $K_E=0$ selects the -[PSS/E-compatible model](https://www.powerworld.com/WebHelp/Content/MainDocumentation_HTML/Transient_Stability_Dialog_Options_Power_System_Model.htm). -At initialization, requiring $\dot E_{fd}'=0$ and -$V_R=V_R^{\max}/10$ gives +$K_E$ is the configured exciter field resistance line slope margin. A nonzero configured value is used directly. The +[IEEE Std 421.4-2014](https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=6803835) states: + +> In some programs, if $K_E$ is entered as zero, $K_E$ is automatically calculated by the program to represent a self-excited shunt field and a trimmed rheostat as its initial condition. + +GridKit implements the +[PSS/E-compatible automatic-parameter rule](https://www.powerworld.com/WebHelp/Content/MainDocumentation_HTML/Transient_Stability_Dialog_Options_Power_System_Model.htm). +The divisor $10$ in $V_R=V_R^{\max}/10=0.1V_R^{\max}$ is unitless and +sets $V_R$ to 10% of the maximum regulator output, retaining the per-unit +units of $V_R^{\max}$: ```math \begin{aligned} From d75babab5a0d6351b79f5a1bb6d9daa0a1cdec33 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Mon, 17 Aug 2026 11:20:35 -0500 Subject: [PATCH 06/10] final adjust --- GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp | 2 +- GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md | 5 +++-- GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp | 2 +- GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md | 5 +++-- 4 files changed, 8 insertions(+), 6 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp index 5b8711916..e4567a165 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp @@ -222,7 +222,7 @@ namespace GridKit * preserving the seeded `efd`, latches attached Known inputs, and * publishes the reference to an attached `vref` signal. * - * @warning IEEE Std 421.4-2014 states: “In some programs, if + * @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 diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md index 75bf163e8..ff7cae735 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md @@ -128,12 +128,13 @@ and when both saturation values are positive, ``` $K_E$ is the configured exciter field resistance line slope margin. A nonzero configured value is used directly. The -[IEEE Std 421.4-2014](https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=6803835) states: +[IEEE Std 421.5-2016](https://standards.ieee.org/ieee/421.5/5356/) states: > In some programs, if $K_E$ is entered as zero, $K_E$ is automatically calculated by the program to represent a self-excited shunt field and a trimmed rheostat as its initial condition. GridKit implements the -[PSS/E-compatible automatic-parameter rule](https://www.powerworld.com/WebHelp/Content/MainDocumentation_HTML/Transient_Stability_Dialog_Options_Power_System_Model.htm). +[PSS/E-compatible automatic-parameter rule](https://www.powerworld.com/WebHelp/Content/MainDocumentation_HTML/Transient_Stability_Dialog_Options_Power_System_Model.htm), as specified in the +[PSS/E 33.4 Program Application Guide, Vol. II](https://pdfcoffee.com/pagv2-pdf-free.html). The divisor $10$ in $V_R=V_R^{\max}/10=0.1V_R^{\max}$ is unitless and sets $V_R$ to 10% of the maximum regulator output, retaining the per-unit units of $V_R^{\max}$: diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp index 570795a7d..6615417d7 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp @@ -190,7 +190,7 @@ namespace GridKit * * Enabled saturation is included via ksat computed from efdp and SA, SB. * - * @warning IEEE Std 421.4-2014 states: “In some programs, if + * @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 diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md index 0ddbede74..b5d7d3a6d 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md @@ -103,12 +103,13 @@ When both saturation values are positive, the non-extraneous solution is: ``` $K_E$ is the configured exciter field resistance line slope margin. A nonzero configured value is used directly. The -[IEEE Std 421.4-2014](https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=6803835) states: +[IEEE Std 421.5-2016](https://standards.ieee.org/ieee/421.5/5356/) states: > In some programs, if $K_E$ is entered as zero, $K_E$ is automatically calculated by the program to represent a self-excited shunt field and a trimmed rheostat as its initial condition. GridKit implements the -[PSS/E-compatible automatic-parameter rule](https://www.powerworld.com/WebHelp/Content/MainDocumentation_HTML/Transient_Stability_Dialog_Options_Power_System_Model.htm). +[PSS/E-compatible automatic-parameter rule](https://www.powerworld.com/WebHelp/Content/MainDocumentation_HTML/Transient_Stability_Dialog_Options_Power_System_Model.htm), as specified in the +[PSS/E 33.4 Program Application Guide, Vol. II](https://pdfcoffee.com/pagv2-pdf-free.html). The divisor $10$ in $V_R=V_R^{\max}/10=0.1V_R^{\max}$ is unitless and sets $V_R$ to 10% of the maximum regulator output, retaining the per-unit units of $V_R^{\max}$: From 4507e138ca00b1632a668344141633298b5c423e Mon Sep 17 00:00:00 2001 From: lukelowry Date: Mon, 17 Aug 2026 12:58:52 -0500 Subject: [PATCH 07/10] move derived parameter and smallest possible initialization order --- .../Exciter/ESDC1A/Esdc1aImpl.hpp | 31 ++++++------- .../PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp | 1 + .../Exciter/IEEET1/Ieeet1Impl.hpp | 29 +++++++----- .../Stabilizer/IEEEST/IeeestImpl.hpp | 33 ++++++++++++-- .../Model/PhasorDynamics/SystemModelImpl.hpp | 45 ++++++++++--------- .../PhasorDynamics/StabilizerIeeestTests.hpp | 2 +- 6 files changed, 90 insertions(+), 51 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp index e4567a165..f7ba5a917 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp @@ -222,16 +222,6 @@ 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 @@ -322,13 +312,14 @@ namespace GridKit const ScalarT se0 = SB_ * Math::qramp(efdp0 - SA_); - RealT ke0 = Ke_; - if (ke0 == ZERO) + RealT ke_eff0 = Ke_eff_; + if (Ke_ == ZERO) { - ke0 = (Vrmax_ / 10.0 - static_cast(se0)) / static_cast(efdp0); + ke_eff0 = (Vrmax_ / 10.0 - static_cast(se0)) + / static_cast(efdp0); } - const ScalarT vfe0 = ke0 * efdp0 + se0; + const ScalarT vfe0 = ke_eff0 * efdp0 + se0; const ScalarT vr0 = vfe0; const ScalarT vhv0 = vr0 / Ka_; @@ -361,7 +352,7 @@ namespace GridKit const ScalarT xll0 = ev0; const ScalarT vref0 = ev0 + vc0 - vs0 - uel_on_ * vuel0; - Ke_eff_ = ke0; + Ke_eff_ = ke_eff0; y[EFDP] = efdp0; y[VC] = vc0; @@ -768,6 +759,16 @@ namespace GridKit * masks let the residual select signal routing without * parameter-dependent control flow, which keeps its structure fixed for * sparse automatic differentiation. + * + * @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. */ template void Esdc1a::setDerivedParameters() diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp index 783483e82..a86518cca 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp @@ -166,6 +166,7 @@ namespace GridKit // Parameter initialization function void initModelParams(const ModelDataT& data); + void setDerivedParameters(); /// Associate variable getter functions with enum values void initializeMonitor(); diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp index 6615417d7..a06e39833 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp @@ -190,15 +190,6 @@ namespace GridKit * * Enabled saturation is included via ksat computed from efdp and SA, SB. * - * @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. */ template int Ieeet1::initialize() @@ -243,7 +234,6 @@ namespace GridKit ScalarT efdp = efd0 / (ONE + omega * Ispdlim_); ScalarT ksat = SB_ * Math::qramp(efdp - SA_); - Ke_eff_ = Ke_; if (Ke_ == ZERO) { Ke_eff_ = (Vrmax_ / 10.0 - static_cast(ksat)) / static_cast(efdp); @@ -478,6 +468,25 @@ namespace GridKit Ispdlim_ = std::get(data.parameters.at(Parameter::Ispdlim)); } + setDerivedParameters(); + } + + /** + * @brief Resolve the parameter-derived constants + * + * @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. + */ + template + void Ieeet1::setDerivedParameters() + { Ke_eff_ = Ke_; Tr_ = std::max(Tr_, TIME_CONSTANT_MINIMUM); diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp index 970907ab0..0b4454775 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp @@ -219,6 +219,15 @@ namespace GridKit template int Ieeest::initialize() { + if (verify() != 0) + { + Log::error() << "Ieeest: cannot initialize with invalid configuration\n"; + return 1; + } + + const ScalarT u = + signals_.template readExternalVariable(); + auto* y = y_.getData(); auto* yp = yp_.getData(); @@ -228,6 +237,22 @@ namespace GridKit yp[static_cast(i)] = 0.0; } + ws_[0] = u; + ws_indices_[0] = + signals_.template readExternalVariableIndex(); + + y[0] = use_notch_ * u; + y[4] = u; + y[5] = u; + y[6] = u; + y[7] = u; + y[8] = u; + y[9] = u; + + // Preserve the current T6 = 0 bypass behavior. + y[10] = bypass_T6_block_ * Ks_ * u; + y[11] = Math::clamp(y[10], Lsmin_, Lsmax_); + y_.setDataUpdated(); yp_.setDataUpdated(); @@ -303,16 +328,18 @@ namespace GridKit ScalarT u = ws[0]; + const ScalarT x2_rhs = (use_4th_order_ + use_3rd_order_) * x3 + + use_2nd_order_ * (-a0_ * x1 - a1_ * x2 + u) * safe_inv_a2_; + f[0] = -x1_dot + use_notch_ * x2; - f[1] = -x2_dot + (use_4th_order_ + use_3rd_order_) * x3 - + use_2nd_order_ * (-a0_ * x1 - a1_ * x2 + u) * safe_inv_a2_; + f[1] = -x2_dot + x2_rhs; f[2] = -x3_dot + use_4th_order_ * x4 + use_3rd_order_ * (-a0_ * x1 - a1_ * x2 - a2_ * x3 + u) * safe_inv_a3_; f[3] = -x4_dot + use_4th_order_ * (-a0_ * x1 - a1_ * x2 - a2_ * x3 - a3_ * x4 + u) * safe_inv_a4_; f[4] = -T2_ * x5_dot - x5 + v4; f[5] = -T4_ * x6_dot - x6 + v5; f[6] = -T6_ * x7_dot - x7 + v6; - f[7] = -v4 + bypass_notch_ * u + use_notch_ * (x1 + A5_ * x2 + (use_4th_order_ + use_3rd_order_) * A6_ * x3); + f[7] = -v4 + bypass_notch_ * u + use_notch_ * (x1 + A5_ * x2 + A6_ * x2_rhs); f[8] = use_T2_block_ * (-T2_ * (v5 - x5) + T1_ * (v4 - x5)) + bypass_T2_block_ * (v4 - v5); f[9] = use_T4_block_ * (-T4_ * (v6 - x6) + T3_ * (v5 - x6)) + bypass_T4_block_ * (v5 - v6); f[10] = use_T6_block_ * (-T6_ * v7 + Ks_ * T5_ * (v6 - x7)) + bypass_T6_block_ * (Ks_ * v6 - v7); diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index 3dd5a2e68..8e640eac1 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -446,6 +446,29 @@ namespace GridKit addComponent(hygov); } + // Add IEEEST stabilizers before exciters that consume their output during + // initialization. + for (const auto& stabdata : data.stabilizer) + { + auto* stabilizer = new Ieeest(stabdata); + + if (stabdata.signal_inputs.contains(IeeestSignalInputs::input)) + { + IdxT input = stabdata.signal_inputs.at(IeeestSignalInputs::input); + constexpr auto U = IeeestExternalVariables::U; + stabilizer->getSignals().template attachSignalNode(getSignal(input)); + } + + if (stabdata.signal_outputs.contains(IeeestSignalOutputs::output)) + { + IdxT output = stabdata.signal_outputs.at(IeeestSignalOutputs::output); + constexpr auto VSS = IeeestInternalVariables::VSS; + stabilizer->getSignals().template assignSignalNode(getSignal(output)); + } + + addComponent(stabilizer); + } + for (const auto& excitedata : data.exciter) { IdxT bus_index = 0; @@ -555,28 +578,6 @@ namespace GridKit addComponent(exciter); } - // Add IEEEST stabilizers - for (const auto& stabdata : data.stabilizer) - { - auto* stabilizer = new Ieeest(stabdata); - - if (stabdata.signal_inputs.contains(IeeestSignalInputs::input)) - { - IdxT input = stabdata.signal_inputs.at(IeeestSignalInputs::input); - constexpr auto U = IeeestExternalVariables::U; - stabilizer->getSignals().template attachSignalNode(getSignal(input)); - } - - if (stabdata.signal_outputs.contains(IeeestSignalOutputs::output)) - { - IdxT output = stabdata.signal_outputs.at(IeeestSignalOutputs::output); - constexpr auto VSS = IeeestInternalVariables::VSS; - stabilizer->getSignals().template assignSignalNode(getSignal(output)); - } - - addComponent(stabilizer); - } - // Add REPCA plant controllers after the signal producers they read at // initialization for (const auto& repcadata : data.repca) diff --git a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp index 14cd9cd88..04662ba54 100644 --- a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp +++ b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp @@ -230,6 +230,7 @@ namespace GridKit std::vector residual_y(residual_y_view.getData(), residual_y_view.getData() + residual_y_view.getSize()); // --- d/dy': tag derivatives as independent --- + u_value = 0.5; stab.initialize(); auto* yp = stab.yp().getData(); for (size_t i = 0; i < stab.size(); ++i) @@ -237,7 +238,6 @@ namespace GridKit yp[i].setVariableNumber(i); } - u_value = 0.5; setStatePointDep(stab); stab.evaluateResidual(); From c5b70ae0bcc905cc9704889d1364c69497ef9ed9 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Mon, 17 Aug 2026 14:38:59 -0500 Subject: [PATCH 08/10] consistant and mroe readable --- .../Exciter/ESDC1A/Esdc1aImpl.hpp | 33 +++++++++---------- .../Exciter/IEEET1/Ieeet1Impl.hpp | 28 ++++++++-------- 2 files changed, 31 insertions(+), 30 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp index f7ba5a917..4789ad481 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp @@ -222,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 @@ -312,14 +322,17 @@ namespace GridKit const ScalarT se0 = SB_ * Math::qramp(efdp0 - SA_); - RealT ke_eff0 = Ke_eff_; if (Ke_ == ZERO) { - ke_eff0 = (Vrmax_ / 10.0 - static_cast(se0)) + Ke_eff_ = (Vrmax_ / 10.0 - static_cast(se0)) / static_cast(efdp0); } + else + { + Ke_eff_ = Ke_; + } - const ScalarT vfe0 = ke_eff0 * efdp0 + se0; + const ScalarT vfe0 = Ke_eff_ * efdp0 + se0; const ScalarT vr0 = vfe0; const ScalarT vhv0 = vr0 / Ka_; @@ -352,8 +365,6 @@ namespace GridKit const ScalarT xll0 = ev0; const ScalarT vref0 = ev0 + vc0 - vs0 - uel_on_ * vuel0; - Ke_eff_ = ke_eff0; - y[EFDP] = efdp0; y[VC] = vc0; y[VR] = vr0; @@ -759,22 +770,10 @@ namespace GridKit * masks let the residual select signal routing without * parameter-dependent control flow, which keeps its structure fixed for * sparse automatic differentiation. - * - * @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. */ template void Esdc1a::setDerivedParameters() { - Ke_eff_ = Ke_; - // The lags are raised to the floor in place, so a negative value is // rejected here while the value as read is still available. verify() // reports the count. diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp index a06e39833..b241f2abf 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp @@ -190,6 +190,15 @@ namespace GridKit * * Enabled saturation is included via ksat computed from efdp and SA, SB. * + * @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. */ template int Ieeet1::initialize() @@ -236,7 +245,12 @@ namespace GridKit ScalarT ksat = SB_ * Math::qramp(efdp - SA_); if (Ke_ == ZERO) { - Ke_eff_ = (Vrmax_ / 10.0 - static_cast(ksat)) / static_cast(efdp); + Ke_eff_ = (Vrmax_ / 10.0 - static_cast(ksat)) + / static_cast(efdp); + } + else + { + Ke_eff_ = Ke_; } ScalarT ve = ksat; @@ -473,22 +487,10 @@ namespace GridKit /** * @brief Resolve the parameter-derived constants - * - * @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. */ template void Ieeet1::setDerivedParameters() { - Ke_eff_ = Ke_; - Tr_ = std::max(Tr_, TIME_CONSTANT_MINIMUM); Ta_ = std::max(Ta_, TIME_CONSTANT_MINIMUM); Te_ = std::max(Te_, TIME_CONSTANT_MINIMUM); From ddc4c314794e856e0d69f79c266092a2168a7ed2 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Mon, 17 Aug 2026 16:41:25 -0500 Subject: [PATCH 09/10] no op, documentation adjustment --- .../Stabilizer/IEEEST/Ieeest.hpp | 30 +++++++-------- .../Stabilizer/IEEEST/IeeestData.hpp | 38 +++++++++---------- .../Stabilizer/IEEEST/IeeestEnzyme.cpp | 2 +- 3 files changed, 35 insertions(+), 35 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp index f77935130..0ed5893e6 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp @@ -32,28 +32,28 @@ namespace GridKit { namespace Stabilizer { - /// Internal variables of a `Ieeest` + /// Internal variables of `Ieeest` enum class IeeestInternalVariables : size_t { - X1, ///< Notch filter state 1 - X2, ///< Notch filter state 2 - X3, ///< Notch filter state 3 - X4, ///< Notch filter state 4 - X5, ///< Lead-lag 1 state - X6, ///< Lead-lag 2 state - X7, ///< Washout state - V4, ///< Notch filter output - V5, ///< Lead-lag 1 output - V6, ///< Lead-lag 2 output - V7, ///< Unlimited stabilizer signal - VSS, ///< Limited stabilizer signal (model output) + X1, ///< \f$x_1\f$ Notch-filter signal state + X2, ///< \f$x_2\f$ First derivative of the filtered signal + X3, ///< \f$x_3\f$ Second derivative of the filtered signal + X4, ///< \f$x_4\f$ Third derivative of the filtered signal + X5, ///< \f$x_5\f$ Lead-lag 1 state + X6, ///< \f$x_6\f$ Lead-lag 2 state + X7, ///< \f$x_7\f$ Washout state + V4, ///< \f$v_4\f$ Notch-filter output + V5, ///< \f$v_5\f$ Lead-lag 1 output + V6, ///< \f$v_6\f$ Lead-lag 2 output + V7, ///< \f$v_7\f$ Unlimited stabilizer signal + VSS, ///< \f$V_{\mathrm{ss}}\f$ Limited stabilizer signal and model output MAXIMUM, }; - /// External variables of a `Ieeest` + /// External variables of an `Ieeest`. enum class IeeestExternalVariables : size_t { - U, ///< Stabilizer input signal + U, ///< \f$u\f$ Stabilizer input signal MAXIMUM, }; diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestData.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestData.hpp index 139e928af..ba8a1466b 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestData.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestData.hpp @@ -19,24 +19,24 @@ namespace GridKit */ enum class IeeestParameters { - A1, ///< Notch filter denominator coefficient - A2, ///< Notch filter denominator coefficient - A3, ///< Notch filter denominator coefficient - A4, ///< Notch filter denominator coefficient - A5, ///< Notch filter numerator coefficient - A6, ///< Notch filter numerator coefficient - T1, ///< Lead-lag 1 numerator time constant - T2, ///< Lead-lag 1 denominator time constant - T3, ///< Lead-lag 2 numerator time constant - T4, ///< Lead-lag 2 denominator time constant - T5, ///< Washout numerator time constant - T6, ///< Washout denominator time constant - Ks, ///< Stabilizer gain - Lsmin, ///< Minimum stabilizer output limit - Lsmax, ///< Maximum stabilizer output limit - Vcl, ///< Lower input cutout threshold (not modeled) - Vcu, ///< Upper input cutout threshold (not modeled) - Tdelay, ///< Input time delay (not modeled) + A1, ///< \f$A_1\f$ Notch-filter denominator coefficient + A2, ///< \f$A_2\f$ Notch-filter denominator coefficient + A3, ///< \f$A_3\f$ Notch-filter denominator coefficient + A4, ///< \f$A_4\f$ Notch-filter denominator coefficient + A5, ///< \f$A_5\f$ Notch-filter numerator coefficient + A6, ///< \f$A_6\f$ Notch-filter numerator coefficient + T1, ///< \f$T_1\f$ Lead-lag 1 numerator time constant + T2, ///< \f$T_2\f$ Lead-lag 1 denominator time constant + T3, ///< \f$T_3\f$ Lead-lag 2 numerator time constant + T4, ///< \f$T_4\f$ Lead-lag 2 denominator time constant + T5, ///< \f$T_5\f$ Washout numerator time constant + T6, ///< \f$T_6\f$ Washout denominator time constant + Ks, ///< \f$K_s\f$ Stabilizer gain + Lsmin, ///< \f$L_s^{\min}\f$ Minimum stabilizer output limit + Lsmax, ///< \f$L_s^{\max}\f$ Maximum stabilizer output limit + Vcl, ///< \f$V_{\mathrm{cl}}\f$ Unsupported lower input-cutout threshold + Vcu, ///< \f$V_{\mathrm{cu}}\f$ Unsupported upper input-cutout threshold + Tdelay, ///< \f$T_{\mathrm{delay}}\f$ Unsupported input delay }; /** @@ -70,7 +70,7 @@ namespace GridKit */ enum class IeeestMonitorableVariables { - vss, ///< Stabilizer output (limited signal) + vss, ///< \f$V_{\mathrm{ss}}\f$ Limited stabilizer signal and model output }; /** diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp index 4c1af73d9..0767a7293 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestEnzyme.cpp @@ -1,7 +1,7 @@ /** * @file IeeestEnzyme.cpp * @author Luke Lowery (lukel@tamu.edu) - * @brief Enzyme-based sparse Jacobian for IEEEST Stabilizer. + * @brief Enzyme sparse Jacobian for the IEEEST stabilizer model. */ #include From 98ddb0ae8c10f05dd418330399e2fd6e8e827089 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Mon, 17 Aug 2026 22:03:35 -0500 Subject: [PATCH 10/10] no-op add warnings as reminders of parameter changes the user may want to know --- .../Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp | 9 +++++++++ .../Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp | 9 ++++++++- 2 files changed, 17 insertions(+), 1 deletion(-) diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp index b241f2abf..c482d4090 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp @@ -247,6 +247,7 @@ namespace GridKit { Ke_eff_ = (Vrmax_ / 10.0 - static_cast(ksat)) / static_cast(efdp); + Log::warning() << "Ieeet1: Ke is zero so effective Ke is derived during initialization\n"; } else { @@ -491,6 +492,14 @@ namespace GridKit template void Ieeet1::setDerivedParameters() { + if (Tr_ < TIME_CONSTANT_MINIMUM || Ta_ < TIME_CONSTANT_MINIMUM + || Te_ < TIME_CONSTANT_MINIMUM || Tf_ < TIME_CONSTANT_MINIMUM) + { + Log::warning() << "Ieeet1: Tr, Ta, Te, and Tf below " + << TIME_CONSTANT_MINIMUM + << " s are raised to that floor\n"; + } + Tr_ = std::max(Tr_, TIME_CONSTANT_MINIMUM); Ta_ = std::max(Ta_, TIME_CONSTANT_MINIMUM); Te_ = std::max(Te_, TIME_CONSTANT_MINIMUM); diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp index 1e582930a..3c2f4c5c9 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp @@ -131,6 +131,12 @@ namespace GridKit template void Tgov1::setDerivedParams() { + if (T1_ < TIME_CONSTANT_MINIMUM || T3_ < TIME_CONSTANT_MINIMUM) + { + Log::warning() << "Tgov1: T1 and T3 below " << TIME_CONSTANT_MINIMUM + << " s are raised to that floor\n"; + } + T1_ = std::max(T1_, TIME_CONSTANT_MINIMUM); T3_ = std::max(T3_, TIME_CONSTANT_MINIMUM); @@ -287,7 +293,8 @@ namespace GridKit const RealT limit_tolerance = static_cast(4.0) * std::numeric_limits::epsilon(); if (pv0_value < Pvmin_ - limit_tolerance || pv0_value > Pvmax_ + limit_tolerance) { - Log::error() << "Tgov1: initial valve position is outside limits\n"; + Log::warning() << "Tgov1: initial valve position is outside [Pvmin, Pvmax]. " + "Check initial dispatch and valve limits\n"; return 1; }