diff --git a/GridKit/AutomaticDifferentiation/DependencyTracking/VariableOperators.hpp b/GridKit/AutomaticDifferentiation/DependencyTracking/VariableOperators.hpp index 5f2ed91ba..de7a29008 100644 --- a/GridKit/AutomaticDifferentiation/DependencyTracking/VariableOperators.hpp +++ b/GridKit/AutomaticDifferentiation/DependencyTracking/VariableOperators.hpp @@ -311,6 +311,12 @@ namespace GridKit return 1.0 / x; } + /// Derivative of log(1 + x). + inline double log1p_derivative(double x) + { + return 1.0 / (1.0 + x); + } + /// Derivative of logarithm to base 10 function: 1/(x*log(10)). inline double log10_derivative(double x) { @@ -377,6 +383,7 @@ namespace std IMPL_FUN_1(tanh, GridKit::DependencyTracking::tanh_derivative) IMPL_FUN_1(exp, GridKit::DependencyTracking::exp_derivative) IMPL_FUN_1(log, GridKit::DependencyTracking::log_derivative) + IMPL_FUN_1(log1p, GridKit::DependencyTracking::log1p_derivative) IMPL_FUN_1(log10, GridKit::DependencyTracking::log10_derivative) IMPL_FUN_1(sqrt, GridKit::DependencyTracking::sqrt_derivative) IMPL_FUN_1(abs, GridKit::DependencyTracking::abs_derivative) diff --git a/GridKit/CommonMath.hpp b/GridKit/CommonMath.hpp index 7077c2076..9bb3a5c3c 100644 --- a/GridKit/CommonMath.hpp +++ b/GridKit/CommonMath.hpp @@ -1,5 +1,6 @@ #pragma once +#include #include #include @@ -29,6 +30,200 @@ namespace GridKit return ONE / (ONE + std::exp(-MU * x)); } + /** + * @brief Smooth one-sided ramp function + * + * Smooth approximation to max(x, 0), using a stable softplus form with + * the same scale as the rest of CommonMath. + * + * @tparam ScalarT - scalar data type + * + * @param[in] x - expected to be of order 1 + * @return value of the smooth ramp function + */ + template + __attribute__((always_inline)) inline ScalarT ramp(const ScalarT x) + { + using RealT = typename GridKit::ScalarTraits::RealT; + static constexpr RealT MU = 240.0; + + ScalarT z = MU * x; + ScalarT a = std::abs(z); + return (HALF * (z + a) + std::log1p(std::exp(-a))) / MU; + } + + /** + * @brief Smooth one-sided quadratic ramp + * + * Smooth approximation to max(x, 0)^2 via a sigmoid-gated quadratic. + * Used for IEEE-style quadratic saturation curves. + * + * @note Eventually a enzyme specialization for an exact implementation + * would be nice, since the piecewise definition is C^1 continuous + * + * @tparam ScalarT - scalar data type + * + * @param[in] x - input signal + * @return value of the quadratic ramp + */ + template + __attribute__((always_inline)) inline ScalarT qramp(const ScalarT x) + { + return x * x * sigmoid(x); + } + + /** + * @brief Smooth binary maximum function + * + * Smooth approximation to max(x, y), composed from the smooth ramp + * function. + * + * @tparam LeftT - scalar type of x + * @tparam RightT - scalar type of y + * + * @param[in] x - First input signal + * @param[in] y - Second input signal + * @return Smooth maximum of x and y + * + * @note The two input types intentionally may differ. Model equations + * often compare a differentiable state or signal with a plain real + * parameter, limit, or literal bound. Keeping both template parameters + * lets the expression promote to the differentiable scalar type without + * forcing callers to cast every parameter. + */ + template + __attribute__((always_inline)) inline auto max( + const LeftT x, + const RightT y) + { + return y + ramp(x - y); + } + + /** + * @brief Smooth binary minimum function + * + * Smooth approximation to min(x, y), composed from the smooth ramp + * function. + * + * @tparam LeftT - scalar type of x + * @tparam RightT - scalar type of y + * + * @param[in] x - First input signal + * @param[in] y - Second input signal + * @return Smooth minimum of x and y + * + * @note The two input types intentionally may differ. Model equations + * often compare a differentiable state or signal with a plain real + * parameter, limit, or literal bound. Keeping both template parameters + * lets the expression promote to the differentiable scalar type without + * forcing callers to cast every parameter. + */ + template + __attribute__((always_inline)) inline auto min( + const LeftT x, + const RightT y) + { + return x - ramp(x - y); + } + + /** + * @brief Smooth clamp function + * + * Smooth approximation to min(max(x, lower), upper), composed from the + * smooth ramp function. Lower and upper bounds may be independent types + * (e.g. constant Real bounds or algebraic-variable bounds). + * + * @tparam ScalarT - scalar data type of the input signal + * @tparam LowerT - data type of the lower bound + * @tparam UpperT - data type of the upper bound + * + * @param[in] x - expected to be of order 1 + * @param[in] lower - Lower limit + * @param[in] upper - Upper limit + * @return value of the smooth clamp function + */ + template + __attribute__((always_inline)) inline auto clamp( + const ScalarT x, + const LowerT lower, + const UpperT upper) + { + assert(lower <= upper); + return lower + ramp(x - lower) - ramp(x - upper); + } + + /** + * @brief Smooth two-sided deadband function + * + * Smooth approximation to x - min(max(x, lower), upper), composed from the + * smooth ramp function. + * + * @tparam ScalarT - scalar data type + * @tparam RealT - Real data type (see GridKit::ScalarTraits::RealT) + * + * @param[in] x - Input signal + * @param[in] lower - Lower breakpoint + * @param[in] upper - Upper breakpoint + * @return Smooth deadbanded value + */ + template + __attribute__((always_inline)) inline ScalarT deadband( + const ScalarT x, + const RealT lower, + const RealT upper) + { + assert(lower <= upper); + return ramp(x - upper) - ramp(-(x - lower)); + } + + /** + * @brief Smooth slew-rate limiter + * + * Smooth approximation to min(max(f, -rate), rate). + * + * @tparam ScalarT - scalar data type + * @tparam RealT - Real data type (see GridKit::ScalarTraits::RealT) + * + * @param[in] f - Pre-limit derivative or rate signal + * @param[in] rate - Symmetric positive rate limit + * @return Slew-rate-limited value of f + */ + template + __attribute__((always_inline)) inline ScalarT slew( + const ScalarT f, + const RealT rate) + { + assert(rate >= ZERO); + return clamp(f, -rate, rate); + } + + /** + * @brief Smooth linear segment contribution + * + * Smooth approximation to a linear segment contribution that is zero below + * lower, linear over [lower, upper], and saturated at height above upper. + * Callers should supply lower < upper; height may be positive or negative. + * + * @tparam ScalarT - scalar data type + * @tparam RealT - Real data type (see GridKit::ScalarTraits::RealT) + * + * @param[in] x - Input signal + * @param[in] lower - Lower breakpoint + * @param[in] upper - Upper breakpoint + * @param[in] height - Saturated value above the upper breakpoint + * @return Smooth linear segment contribution + */ + template + __attribute__((always_inline)) inline ScalarT linseg( + const ScalarT x, + const RealT lower, + const RealT upper, + const RealT height) + { + assert(lower < upper); + return height / (upper - lower) * (ramp(x - lower) - ramp(x - upper)); + } + /** * @brief Derivative of the scaled sigmoid activation function * (i.e., approximation to the delta dirac function) @@ -46,86 +241,138 @@ namespace GridKit } /** - * @brief Low indicator function for regulator limits + * @brief Smooth above-limit indicator * * @tparam ScalarT - Scalar data type * @tparam RealT - Real data type (see GridKit::ScalarTraits::RealT) * - * @param[in] limit_min - Minimum limit * @param[in] x - State variable - * @return Scalar value indicating limit activation + * @param[in] limit_min - Minimum limit + * @return Smooth indicator that x is above limit_min */ template - __attribute__((always_inline)) inline ScalarT indicator_low( - const RealT limit_min, - const ScalarT x) + __attribute__((always_inline)) inline ScalarT above( + const ScalarT x, + const RealT limit_min) { return sigmoid(x - limit_min); } /** - * @brief High indicator function for regulator limits + * @brief Smooth below-limit indicator * * @tparam ScalarT - Scalar data type * @tparam RealT - Real data type (see GridKit::ScalarTraits::RealT) * - * @param[in] limit_max - Maximum limit * @param[in] x - State variable - * @return Scalar value indicating limit activation + * @param[in] limit_max - Maximum limit + * @return Smooth indicator that x is below limit_max */ template - __attribute__((always_inline)) inline ScalarT indicator_high( - const RealT limit_max, - const ScalarT x) + __attribute__((always_inline)) inline ScalarT below( + const ScalarT x, + const RealT limit_max) { return sigmoid(limit_max - x); } /** - * @brief Zero indicator function for regulator limits + * @brief Smooth inside-limits indicator * * @tparam ScalarT - Scalar data type * @tparam RealT - Real data type (see GridKit::ScalarTraits::RealT) * - * @param[in] limit_max - Maximum limit * @param[in] x - State variable - * @return Scalar value indicating limit activation + * @param[in] limit_min - Minimum limit + * @param[in] limit_max - Maximum limit + * @return Smooth indicator that x is inside [limit_min, limit_max] */ template - __attribute__((always_inline)) inline ScalarT indicator_zero( + __attribute__((always_inline)) inline ScalarT inside( + const ScalarT x, const RealT limit_min, - const RealT limit_max, - const ScalarT x) + const RealT limit_max) { - return indicator_low(limit_min, x) + indicator_high(limit_max, x) - ONE; + assert(limit_min <= limit_max); + return above(x, limit_min) + below(x, limit_max) - ONE; } /** - * @brief Smooth anti-windup indicator for a limited state variable + * @brief Smooth outside-limits indicator * * @tparam ScalarT - Scalar data type * @tparam RealT - Real data type (see GridKit::ScalarTraits::RealT) * + * @param[in] x - State variable * @param[in] limit_min - Minimum limit * @param[in] limit_max - Maximum limit + * @return Smooth indicator that x is outside [limit_min, limit_max] + */ + template + __attribute__((always_inline)) inline ScalarT outside( + const ScalarT x, + const RealT limit_min, + const RealT limit_max) + { + assert(limit_min <= limit_max); + return below(x, limit_min) + above(x, limit_max); + } + + /** + * @brief Smooth anti-windup indicator for a limited state variable + * + * @tparam ScalarT - Scalar data type + * @tparam RealT - Real data type (see GridKit::ScalarTraits::RealT) + * * @param[in] x - State variable * @param[in] f - Pre-limit derivative of the state variable + * @param[in] limit_min - Minimum limit + * @param[in] limit_max - Maximum limit * @return Scalar value in [0, 1]: 1 when dynamics should pass through, * 0 when integration should be blocked. */ template __attribute__((always_inline)) inline ScalarT indicator( - const RealT limit_min, - const RealT limit_max, const ScalarT x, - const ScalarT f) + const ScalarT f, + const RealT limit_min, + const RealT limit_max) { - ScalarT above_min = indicator_low(limit_min, x); - ScalarT below_max = indicator_high(limit_max, x); + assert(limit_min <= limit_max); + + ScalarT above_min = above(x, limit_min); + ScalarT below_max = below(x, limit_max); return above_min * below_max + // (ONE - below_max) * sigmoid(-f) + // (ONE - above_min) * sigmoid(f); } + + /** + * @brief Smooth anti-windup limited derivative + * + * Applies the smooth anti-windup indicator gate to a pre-limit derivative. + * The returned value approximates the conditional-integration rule that + * passes interior dynamics, passes restoring motion from saturated limits, + * and blocks motion that would push further into saturation. + * + * @tparam ScalarT - Scalar data type + * @tparam RealT - Real data type (see GridKit::ScalarTraits::RealT) + * + * @param[in] x - Limited state or limited output signal + * @param[in] f - Pre-limit derivative + * @param[in] limit_min - Minimum limit + * @param[in] limit_max - Maximum limit + * @return Smooth anti-windup limited derivative + */ + template + __attribute__((always_inline)) inline ScalarT antiwindup( + const ScalarT x, + const ScalarT f, + const RealT limit_min, + const RealT limit_max) + { + return indicator(x, f, limit_min, limit_max) * f; + } } // namespace Math } // namespace GridKit diff --git a/GridKit/CommonMath.md b/GridKit/CommonMath.md index c4369c303..a654487a2 100644 --- a/GridKit/CommonMath.md +++ b/GridKit/CommonMath.md @@ -1,62 +1,171 @@ # CommonMath -Numerical utilities in [CommonMath.hpp](CommonMath.hpp): smooth, autodiff-friendly replacements for piecewise functions used across GridKit component models. +Smooth, autodiff-friendly replacements for piecewise functions used across GridKit component models. See [CommonMath.hpp](CommonMath.hpp) for implementation details. -## Sigmoid +## Primitives -Smooth approximation to the step function. +### Descriptions -```math -\sigma(x) = \dfrac{1}{1+\exp(-\alpha x)} -``` +| Name | Description | Usage | +|------|-------------|-------| +| `sigmoid` | Step function | `IEEET1` | +| `ramp` | Smooth one-sided ramp | `REGCA` | +| `qramp` | Exact one-sided quadratic ramp | `IEEET1` | -The scale $\alpha$ (currently $240$) is chosen large enough that $\sigma$ behaves as a step on inputs of order 1. +### Exact Equations -## Limit Indicators +```math +\begin{aligned} +\sigma(x) +&= +\begin{cases} +0 & x \le 0 \\ +1 & x > 0 +\end{cases} +\\ +\rho(x) &= x\,\sigma(x) \\ +q(x) &= x^2\,\sigma(x) +\end{aligned} +``` -For a state variable $x$ with limits $(x_{\min}, x_{\max})$: +### Smooth Equations ```math \begin{aligned} - \phi_L(x) &= \sigma(x - x_{\min}) \\ - \phi_U(x) &= \sigma(x_{\max} - x) \\ - \phi_0(x) &= \phi_L + \phi_U - 1 +\sigma(x) &= \dfrac{1}{1+e^{-\mu x}} \\ +\rho(x) &= \dfrac{(\mu x+\lvert\mu x\rvert)/2+\log(1+e^{-\lvert\mu x\rvert})}{\mu} \\ +q(x) &= x^2\,\sigma(x) \end{aligned} ``` -$\phi_L$ and $\phi_U$ are soft indicators for "above lower limit" and "below upper limit"; their product (or $\phi_0$) is an interior indicator. +The scale $\mu=4\cdot f_{\text{sync}}=240$ is chosen so $\sigma$ behaves like a step on inputs of order 1 while keeping derivatives finite. As $\mu \to \infty$, these functions approach their exact targets. (*Note*: the implementation of the quadratic ramp `q(x)` could be optimized with Enzyme features down the road). -## Anti-Windup Indicator +## Derived Functions -Component models with a limited state $x \in (x_{\min}, x_{\max})$ and pre-limit derivative $f$ express the gated dynamics piecewise as +### Descriptions -```math -\dot x = - \begin{cases} - f - & \text{if } (x_{\min} < x < x_{\max}) & \lor \\ - & \quad (x \leq x_{\min} \land f > 0) & \lor \\ - & \quad (x \geq x_{\max} \land f < 0) \\ - 0 & \text{else} - \end{cases} -``` +| Name | Description | Usage | +|------|-------------|-------| +| `max` | Smooth binary maximum | `REECA` | +| `min` | Smooth binary minimum | `REECA` | +| `clamp` | Bounded saturation | `IEEEST` | +| `deadband` | Signed two-sided deadband | `REECA` | +| `slew` | Symmetric slew-rate limiter | - | +| `linseg` | Saturated linear segment contribution | `REGCA`, `REECA` | +| `above` | Above-lower-limit indicator | - | +| `below` | Below-upper-limit indicator | - | +| `inside` | Interior pulse indicator | `IEEEST` | +| `outside` | Outside-band indicator | `REECA` | +| `antiwindup` | Anti-windup limited derivative | `IEEET1`, `TGOV1`, `SEXS-PTI`, `REECA` | -In simulation this is replaced with the smooth approximation $\dot x = \phi(x, f) \cdot f$, where +### Exact Equations ```math -\phi(x, f) = \phi_L \phi_U + (1 - \phi_U)\,\sigma(-f) + (1 - \phi_L)\,\sigma(f). +\begin{aligned} +\text{max}(x,y) +&= +\begin{cases} +x & x > y \\ +y & x \le y +\end{cases} +\\ +\text{min}(x,y) +&= +\begin{cases} +x & x < y \\ +y & x \ge y +\end{cases} +\\ +\text{clamp}(x;\ell,u) +&= +\begin{cases} +\ell & x < \ell \\ +x & \ell \le x \le u \\ +u & x > u +\end{cases} +\\ +\text{deadband}(x;\ell,u) +&= +\begin{cases} +x-\ell & x < \ell \\ +0 & \ell \le x \le u \\ +x-u & x > u +\end{cases} +\\ +\text{slew}(f;r) +&= +\begin{cases} +-r & f < -r \\ +f & -r \le f \le r \\ +r & f > r +\end{cases} +\\ +\text{linseg}(x;a,b,h) +&= +\begin{cases} +0 & x < a \\ +\dfrac{h}{b-a}(x-a) & a \le x \le b \\ +h & x > b +\end{cases} +\\ +\text{above}(x;\ell) +&= +\begin{cases} +0 & x \le \ell \\ +1 & x > \ell +\end{cases} +\\ +\text{below}(x;u) +&= +\begin{cases} +1 & x < u \\ +0 & x \ge u +\end{cases} +\\ +\text{inside}(x;\ell,u) +&= +\begin{cases} +1 & \ell < x < u \\ +0 & \text{else} +\end{cases} +\\ +\text{outside}(x;\ell,u) +&= +\begin{cases} +1 & x < \ell \lor x > u \\ +0 & \text{else} +\end{cases} +\\ +\text{antiwindup}(x,f;\ell,u) +&= +\begin{cases} +f & \ell < x < u \\ +f & x \le \ell \land f > 0 \\ +f & x \ge u \land f < 0 \\ +0 & \text{otherwise} +\end{cases} +\end{aligned} ``` -The first term passes the dynamics through in the interior. The second re-admits them when $x$ is above $x_{\max}$ and $f$ is pulling back down; the third re-admits them when $x$ is below $x_{\min}$ and $f$ is driving back up. When $x$ is outside a limit *and* $f$ is still driving further beyond it, $\phi$ vanishes smoothly and blocks the windup. - -## Callers - -Anti-windup indicator (`Math::indicator`): - -- [IEEET1](Model/PhasorDynamics/Exciter/IEEET1/README.md): gates $\dot V_R$ on $V_R \in (V_{rmin}, V_{rmax})$ -- [TGOV1](Model/PhasorDynamics/Governor/Tgov1/README.md): gates $\dot P_v$ on $P_v \in (P_{vmin}, P_{vmax})$ -- [SEXS-PTI](Model/PhasorDynamics/Exciter/SEXS-PTI/README.md): gates $\dot E_{fd}$ on $E_{fd} \in (E_{fd,\min}, E_{fd,\max})$ +### Smooth Equations -Interior indicator (`Math::indicator_zero`): + -- [IEEEST](Model/PhasorDynamics/Stabilizer/IEEEST/README.md): clips the stabilizer output $v_7$ to $[L_{s\min}, L_{s\max}]$, passing $v_7$ through in the interior and saturating to the limits outside +```math +\begin{aligned} +\text{max}(x,y) &= y+\rho(x-y) \\ +\text{min}(x,y) &= x-\rho(x-y) \\ +\text{clamp}(x;\ell,u) &= \ell+\rho(x-\ell)-\rho(x-u) \\ +\text{deadband}(x;\ell,u) &= \rho(x-u)-\rho(\ell-x) \\ +\text{slew}(f;r) &= -r+\rho(f+r)-\rho(f-r) \\ +\text{linseg}(x;a,b,h) &= \dfrac{h}{b-a}\left[\rho(x-a)-\rho(x-b)\right] \\ +\text{above}(x;\ell) &= \sigma(x-\ell) \\ +\text{below}(x;u) &= \sigma(u-x) \\ +\text{inside}(x;\ell,u) &= \sigma(x-\ell)+\sigma(u-x)-1 \\ +\text{outside}(x;\ell,u) &= \sigma(\ell-x)+\sigma(x-u) \\ +\phi_L &= \text{above}(x;\ell) \\ +\phi_U &= \text{below}(x;u) \\ +\phi(x,f) &= \phi_L\phi_U+(1-\phi_U)\sigma(-f)+(1-\phi_L)\sigma(f) \\ +\text{antiwindup}(x,f;\ell,u) &= \phi(x,f)f +\end{aligned} +``` diff --git a/GridKit/Model/PhasorDynamics/Converter/README.md b/GridKit/Model/PhasorDynamics/Converter/README.md new file mode 100644 index 000000000..ad38ba19a --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/README.md @@ -0,0 +1,14 @@ +# **Converter Models** + +## Introduction + +Converter models represent inverter-coupled resources in the phasor dynamics model. They provide the network interface between renewable-energy control +models and the bus equations, typically through commanded active and reactive current components. + +## Types + +The GridKit converter documentation includes: + +- Renewable Energy Generator/Converter Model REGCA (See [REGCA](REGCA/README.md)) +- Renewable Energy Generator/Converter Model REGCB (See [REGCB](REGCB/README.md)) +- Renewable Energy Electrical Control Model REECA (See [REECA](REECA/README.md)) diff --git a/GridKit/Model/PhasorDynamics/Converter/REECA/README.md b/GridKit/Model/PhasorDynamics/Converter/REECA/README.md new file mode 100644 index 000000000..6e0d0b83e --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REECA/README.md @@ -0,0 +1,426 @@ +# **Renewable Energy Electrical Control Model (REECA)** + +REECA is a WECC renewable energy electrical control model for inverter-coupled resources. In GridKit it is represented as a signal-control model that computes active- and reactive-current commands. + +Notes: +- Internal electrical quantities and current commands are on model base unless otherwise stated. +- Optional signal inputs default to their documented constant values when omitted. +- Timer-based post-dip reactive-current injection hold and active-current limit hold are not modeled in this version; $T_{\mathrm{hld}}$ and $T_{\mathrm{hld2}}$ must be zero. + +## Block Diagram + +Standard REECA block diagram. + +
+ + + Figure 1: REECA block diagram. Figure courtesy of [PowerWorld](https://www.powerworld.com/WebHelp/) +
+ +## Model Parameters + +Symbol | Units | Description | Typical Value | Note +------------------------------------|----------|---------------------------------------------------------|---------------|------ +$S^{\mathrm{base}}$ | [MVA] | REECA model power base | TBD | Block name: `MVABase` +$s_{\mathrm{pf}}$ | [binary] | Power-factor control flag | TBD | Block name: `PfFlag`; 1 = power-factor control, 0 = Q control +$s_V$ | [binary] | Voltage-control mode flag | TBD | Block name: `VFlag`; 1 = Q control, 0 = voltage control +$s_Q$ | [binary] | Reactive-power control flag | TBD | Block name: `QFlag`; 1 = voltage/Q control, 0 = constant pf or Q control +$s_P$ | [binary] | Active-power reference speed-multiplier flag | TBD | Block name: `Pflag`; 1 = multiply by generator speed +$s_{PQ}$ | [binary] | P/Q priority flag for converter current limit | TBD | Block name: `Pqflag`; 0 = Q priority, 1 = P priority +$T_{\mathrm{rv}}$ | [sec] | Voltage-measurement filter time constant | TBD | Block name: `Trv`; if zero, $V_{\mathrm{meas}}$ is algebraic +$T_{\mathrm{p}}$ | [sec] | Electrical-power measurement filter time constant | TBD | Block name: `Tp`; if zero, $P_{\mathrm{meas}}$ is algebraic +$V_{\mathrm{ref0}}$ | [p.u.] | Outer-loop voltage reference | TBD | Block name: `Vref0`; initialized to terminal voltage if omitted +$V_{\mathrm{dip}}$ | [p.u.] | Low-voltage threshold for reactive-current injection logic | TBD | Block name: `Vdip` +$V_{\mathrm{up}}$ | [p.u.] | High-voltage threshold for reactive-current injection logic | TBD | Block name: `Vup` +$D_{\mathrm{bd1}}$ | [p.u.] | Overvoltage deadband for voltage-error response | TBD | Block name: `dbd1` +$D_{\mathrm{bd2}}$ | [p.u.] | Undervoltage deadband for voltage-error response | TBD | Block name: `dbd2` +$K_{\mathrm{qv}}$ | [p.u.] | Reactive-current injection gain during voltage dip/overvoltage logic | TBD | Block name: `kqv` +$I_{\mathrm{qinj}}^{\min}$ | [p.u.] | Minimum reactive-current injection limit | TBD | Block name: `Iql1` +$I_{\mathrm{qinj}}^{\max}$ | [p.u.] | Maximum reactive-current injection limit | TBD | Block name: `Iqh1` +$I_{\mathrm{qinj}}^{\mathrm{frz}}$ | [p.u.] | Held reactive-current injection value after voltage dip | TBD | Block name: `Iqfrz`; unused when $T_{\mathrm{hld}} = 0$ +$T_{\mathrm{hld}}$ | [sec] | Reactive-current injection hold time after voltage dip clears | TBD | Block name: `Thld`; required to be zero in this version +$Q^{\max}$ | [p.u.] | Maximum reactive-power control limit | TBD | Block name: `Qmax` +$Q^{\min}$ | [p.u.] | Minimum reactive-power control limit | TBD | Block name: `Qmin` +$K_{\mathrm{qp}}$ | [p.u.] | Reactive-power control proportional gain | TBD | Block name: `Kqp` +$K_{\mathrm{qi}}$ | [p.u./s] | Reactive-power control integral gain | TBD | Block name: `Kqi` +$V^{\max}$ | [p.u.] | Maximum voltage-control limit | TBD | Block name: `Vmax` +$V^{\min}$ | [p.u.] | Minimum voltage-control limit | TBD | Block name: `Vmin` +$V_{\mathrm{ref1}}$ | [p.u.] | Inner-loop voltage-control reference/bias | 0 | Block name: `Vref1` +$K_{\mathrm{vp}}$ | [p.u.] | Voltage-control proportional gain | TBD | Block name: `Kvp` +$K_{\mathrm{vi}}$ | [p.u./s] | Voltage-control integral gain | TBD | Block name: `Kvi` +$T_{\mathrm{iq}}$ | [sec] | Reactive-current command lag time constant | TBD | Block name: `Tiq` +$T_{\mathrm{pord}}$ | [sec] | Active-power order filter time constant | TBD | Block name: `Tpord` +$R_P^{\max}$ | [p.u./s] | Positive active-power order ramp-rate limit | TBD | Block name: `dPmax` +$R_P^{\min}$ | [p.u./s] | Negative active-power order ramp-rate limit | TBD | Block name: `dPmin` +$P^{\max}$ | [p.u.] | Maximum active-power order limit | TBD | Block name: `Pmax` +$P^{\min}$ | [p.u.] | Minimum active-power order limit | TBD | Block name: `Pmin` +$I^{\max}$ | [p.u.] | Maximum total converter current | TBD | Block name: `Imax` +$V_{\mathrm{q},1}$ | [p.u.] | VDL1 voltage point 1 | TBD | Block name: `vq1` +$I_{\mathrm{q},1}^{\max}$ | [p.u.] | VDL1 reactive-current limit point 1 | TBD | Block name: `lq1` +$V_{\mathrm{q},2}$ | [p.u.] | VDL1 voltage point 2 | TBD | Block name: `vq2` +$I_{\mathrm{q},2}^{\max}$ | [p.u.] | VDL1 reactive-current limit point 2 | TBD | Block name: `lq2` +$V_{\mathrm{q},3}$ | [p.u.] | VDL1 voltage point 3 | TBD | Block name: `vq3` +$I_{\mathrm{q},3}^{\max}$ | [p.u.] | VDL1 reactive-current limit point 3 | TBD | Block name: `lq3` +$V_{\mathrm{q},4}$ | [p.u.] | VDL1 voltage point 4 | TBD | Block name: `vq4` +$I_{\mathrm{q},4}^{\max}$ | [p.u.] | VDL1 reactive-current limit point 4 | TBD | Block name: `lq4` +$V_{\mathrm{p},1}$ | [p.u.] | VDL2 voltage point 1 | TBD | Block name: `vp1` +$I_{\mathrm{p},1}^{\max}$ | [p.u.] | VDL2 active-current limit point 1 | TBD | Block name: `lp1` +$V_{\mathrm{p},2}$ | [p.u.] | VDL2 voltage point 2 | TBD | Block name: `vp2` +$I_{\mathrm{p},2}^{\max}$ | [p.u.] | VDL2 active-current limit point 2 | TBD | Block name: `lp2` +$V_{\mathrm{p},3}$ | [p.u.] | VDL2 voltage point 3 | TBD | Block name: `vp3` +$I_{\mathrm{p},3}^{\max}$ | [p.u.] | VDL2 active-current limit point 3 | TBD | Block name: `lp3` +$V_{\mathrm{p},4}$ | [p.u.] | VDL2 voltage point 4 | TBD | Block name: `vp4` +$I_{\mathrm{p},4}^{\max}$ | [p.u.] | VDL2 active-current limit point 4 | TBD | Block name: `lp4` +$T_{\mathrm{hld2}}$ | [sec] | Active-current limit hold time after voltage dip clears | TBD | Block name: `Thld2`; required to be zero in this version + +### Parameter Validation + +Implementations should reject invalid REECA parameter sets. If source data preprocessing adjusts active-power ramp-rate or order limits, apply these checks to the effective values used by the equations. + +The required checks are: + +```math +\begin{aligned} + &S^{\mathrm{base}} > 0 \\ + &s_{\mathrm{pf}}, s_V, s_Q, s_P, s_{PQ} \in \{0,1\} \\ + &T_{\mathrm{rv}}, T_{\mathrm{p}} \ge 0 \\ + &0 \le V_{\mathrm{dip}} < V_{\mathrm{up}} \\ + &D_{\mathrm{bd1}} \le 0 \le D_{\mathrm{bd2}} \\ + &I_{\mathrm{qinj}}^{\min} \le I_{\mathrm{qinj}}^{\max} \\ + &T_{\mathrm{hld}} = T_{\mathrm{hld2}} = 0 \\ + &Q^{\min} \le Q^{\max} \\ + &V^{\min} \le V^{\max} \\ + &T_{\mathrm{iq}}, T_{\mathrm{pord}} > 0 \\ + &R_P^{\min} < 0 < R_P^{\max} \\ + &P^{\min} \le P^{\max} \\ + &I^{\max} \ge 0 \\ + &0 \le V_{\mathrm{q},1} < V_{\mathrm{q},2} < V_{\mathrm{q},3} < V_{\mathrm{q},4} \\ + &I_{\mathrm{q},k}^{\max} \ge 0\ \text{for } k=1,\ldots,4 \\ + &0 \le V_{\mathrm{p},1} < V_{\mathrm{p},2} < V_{\mathrm{p},3} < V_{\mathrm{p},4} \\ + &I_{\mathrm{p},k}^{\max} \ge 0\ \text{for } k=1,\ldots,4 +\end{aligned} +``` + +### Model Derived Parameters + +The off-mode flag complements are: + +```math +\begin{aligned} + s_{\mathrm{pf}}^{\mathrm{off}} &= 1 - s_{\mathrm{pf}} \\ + s_V^{\mathrm{off}} &= 1 - s_V \\ + s_Q^{\mathrm{off}} &= 1 - s_Q \\ + s_{PQ}^{\mathrm{off}} &= 1 - s_{PQ} +\end{aligned} +``` + +The VDL functions use GridKit's smooth [Linear Segment](../../../../CommonMath.md#derived-functions) helper and provide flat extrapolation outside the first and fourth voltage points: + +```math +\begin{aligned} + g_q(x) &= + I_{\mathrm{q},1}^{\max} + + \sum_{k=1}^{3} + \text{linseg}\!\left( + x;\, + V_{\mathrm{q},k},\, + V_{\mathrm{q},k+1},\, + I_{\mathrm{q},k+1}^{\max} - I_{\mathrm{q},k}^{\max} + \right) \\ + g_p(x) &= + I_{\mathrm{p},1}^{\max} + + \sum_{k=1}^{3} + \text{linseg}\!\left( + x;\, + V_{\mathrm{p},k},\, + V_{\mathrm{p},k+1},\, + I_{\mathrm{p},k+1}^{\max} - I_{\mathrm{p},k}^{\max} + \right) +\end{aligned} +``` + +## Model Variables + +### Internal Variables + +#### Differential + +Symbol | Units | Description | Note +------------------------|--------|-------------------------------------|------ +$V_{\mathrm{meas}}$ | [p.u.] | Filtered terminal voltage | State 1 in Fig. 1; source label: `Vmeas`; algebraic when $T_{\mathrm{rv}} = 0$ +$P_{\mathrm{meas}}$ | [p.u.] | Filtered electrical power | State 2 in Fig. 1; source label: `Pmeas`; algebraic when $T_{\mathrm{p}} = 0$ +$x_{\mathrm{PIQ}}$ | [p.u.] | Reactive-power PI controller state | State 3 in Fig. 1; source label: `PIQ` +$x_{\mathrm{PIV}}$ | [p.u.] | Voltage PI controller state | State 4 in Fig. 1; source label: `PIV` +$Q_V$ | [p.u.] | Reactive-current command lag state | State 5 in Fig. 1; source label: `Q_V` +$P_{\mathrm{ord}}$ | [p.u.] | Filtered active-power order | State 6 in Fig. 1; source label: `Pord` + +#### Algebraic + +Symbol | Units | Description | Note +--------------------------------|--------|-------------------------------------|------ +$V_T$ | [p.u.] | Terminal voltage magnitude | +$V_{\mathrm{meas}}^{\mathrm{safe}}$ | [p.u.] | Safe filtered terminal voltage for divider blocks | Lower bounded by 0.01 +$s_{\mathrm{dip}}$ | [binary] | Voltage-dip/overvoltage freeze indicator | 1 when outside voltage thresholds +$V_{\mathrm{err}}$ | [p.u.] | Deadbanded voltage error | Defined by CommonMath `deadband` +$I_{\mathrm{qv}}$ | [p.u.] | Reactive-current injection candidate | Converter base +$Q_{\mathrm{ref}}$ | [p.u.] | Selected reactive-power reference | From power-factor or external reactive-power command +$e_Q$ | [p.u.] | Reactive-power control error | Limited $Q_{\mathrm{ref}}$ minus $Q_{\mathrm{gen}}$ +$V_{\mathrm{PIQ}}$ | [p.u.] | Reactive-power control PI output | Limited by $V^{\min}$ and $V^{\max}$ +$e_{\mathrm{PIV}}$ | [p.u.] | Voltage-control PI error | Selected voltage-control signal minus $V_{\mathrm{meas}}$ +$f_{\mathrm{pord}}$ | [p.u./s] | Active-power order derivative before ramp-rate limiting | Feeds $r_{\mathrm{pord}}$ +$r_{\mathrm{pord}}$ | [p.u./s] | Ramp-rate-limited active-power order derivative | Feeds $P_{\mathrm{ord}}$ anti-windup +$I_{\mathrm{q}}^{\mathrm{circ}}$ | [p.u.] | Reactive-current limit from converter current circle | Converter base; nonnegative algebraic branch +$I_{\mathrm{p}}^{\mathrm{circ}}$ | [p.u.] | Active-current limit from converter current circle | Converter base; nonnegative algebraic branch +$I_{\mathrm{q}}^{\max}$ | [p.u.] | Final reactive-current upper limit | Converter base; updated by VDL1 and current-limit logic +$I_{\mathrm{p}}^{\max}$ | [p.u.] | Final active-current upper limit | Converter base; updated by VDL2 and current-limit logic +$I_{\mathrm{qbase}}$ | [p.u.] | Base reactive-current command | Converter base; before $s_Q$ selection and reactive-current injection +$I_{\mathrm{q}}^{\mathrm{raw}}$ | [p.u.] | Raw reactive-current command before final limit | Converter base +$I_{\mathrm{q}}^{\mathrm{cmd}}$ | [p.u.] | Reactive-current command output | Converter base +$I_{\mathrm{p}}^{\mathrm{cmd}}$ | [p.u.] | Active-current command output | Converter base + +### External Variables + +#### Differential + +Symbol | Units | Description | Note +-----------|--------|-------------------------|------ +$\omega$ | [p.u.] | Generator speed deviation | Optional, defaults to zero; source diagram $\omega_g = 1 + \omega$ + +#### Algebraic + +Symbol | Units | Description | Note +-------------------------------------|--------|-----------------------------------|------ +$V_{\mathrm{r}}$ | [p.u.] | Terminal voltage, real component | Owned by bus object +$V_{\mathrm{i}}$ | [p.u.] | Terminal voltage, imaginary component | Owned by bus object +$P_e$ | [p.u.] | Electrical active power | Source label: `Pe` +$Q_{\mathrm{gen}}$ | [p.u.] | Reactive-power feedback | Source label: `Qgen` +$Q_{\mathrm{ext}}$ | [p.u.] | External reactive-power command | Optional, defaults to initialized constant +$\phi_{\mathrm{pf}}^{\mathrm{ref}}$ | [rad] | Power-factor angle reference | Source label: `pfaref`; used through tangent block +$P_{\mathrm{ref}}$ | [p.u.] | External active-power reference | Optional, defaults to initialized constant + +## Model Equations + +For readability, define: + +```math +\begin{aligned} + f_{\mathrm{PIQ}} &= K_{\mathrm{qi}} e_Q \\ + f_{\mathrm{PIV}} &= K_{\mathrm{vi}} e_{\mathrm{PIV}} +\end{aligned} +``` + +### Differential Equations + +The state-equation residuals use compact limiter notation where applicable. The measurement filters are written in descriptor form: if $T_{\mathrm{rv}} = 0$ or $T_{\mathrm{p}} = 0$, the corresponding variable should be tagged algebraic. The $Q_V$ equation also uses $T_{\mathrm{iq}}$ as a derivative coefficient, but $T_{\mathrm{iq}} > 0$ remains required because the freeze multiplier makes the zero-time case structurally different. + +```math +\begin{aligned} + 0 &= -T_{\mathrm{rv}}\dot V_{\mathrm{meas}} - V_{\mathrm{meas}} + V_T \\ + 0 &= -T_{\mathrm{p}}\dot P_{\mathrm{meas}} - P_{\mathrm{meas}} + P_e \\ + 0 &= + -\dot x_{\mathrm{PIQ}} + + (1 - s_{\mathrm{dip}}) + \text{antiwindup}\!\left( + V_{\mathrm{PIQ}}, + f_{\mathrm{PIQ}}, + V^{\min}, + V^{\max} + \right) \\ + 0 &= + -\dot x_{\mathrm{PIV}} + + (1 - s_{\mathrm{dip}}) + \text{antiwindup}\!\left( + I_{\mathrm{qbase}}, + f_{\mathrm{PIV}}, + -I_{\mathrm{q}}^{\max}, + I_{\mathrm{q}}^{\max} + \right) \\ + 0 &= + -T_{\mathrm{iq}}\dot Q_V + - (1 - s_{\mathrm{dip}})Q_V + + (1 - s_{\mathrm{dip}})Q_{\mathrm{ref}}/V_{\mathrm{meas}}^{\mathrm{safe}} \\ + 0 &= + -\dot P_{\mathrm{ord}} + + (1 - s_{\mathrm{dip}}) + \text{antiwindup}\!\left( + P_{\mathrm{ord}}, + r_{\mathrm{pord}}, + P^{\min}, + P^{\max} + \right) +\end{aligned} +``` + +CommonMath defines the [Anti-Windup](../../../../CommonMath.md#anti-windup-indicator) target and smooth approximation. + +### Algebraic Equations + +The algebraic targets use CommonMath helper notation where applicable: + +```math +\begin{aligned} + 0 &= -V_T^2 + V_\mathrm r^2 + V_\mathrm i^2 \\ + 0 &= -V_\mathrm{meas}^\mathrm{safe} + \max(V_\mathrm{meas}, 0.01) \\ + 0 &= -s_\mathrm{dip} + \text{outside}(V_T, V_\mathrm{dip}, V_\mathrm{up}) \\ + 0 &= -V_\mathrm{err} + \text{deadband}(V_\mathrm{ref0} - V_\mathrm{meas}, D_\mathrm{bd1}, D_\mathrm{bd2}) \\ + 0 &= -I_\mathrm{qv} + \text{clamp}(K_\mathrm{qv} V_\mathrm{err}, I_\mathrm{qinj}^{\min}, I_\mathrm{qinj}^{\max}) \\ + 0 &= -Q_\mathrm{ref} + + s_\mathrm{pf} P_\mathrm{meas}\tan(\phi_\mathrm{pf}^\mathrm{ref}) + + s_\mathrm{pf}^\mathrm{off} Q_\mathrm{ext} \\ + 0 &= -e_Q + \text{clamp}(Q_\mathrm{ref}, Q^{\min}, Q^{\max}) - Q_\mathrm{gen} \\ + 0 &= -V_\mathrm{PIQ} + \text{clamp}(K_\mathrm{qp} e_Q + x_\mathrm{PIQ}, V^{\min}, V^{\max}) \\ + 0 &= -e_\mathrm{PIV} + s_V V_\mathrm{PIQ} + s_V^\mathrm{off}(Q_\mathrm{ref} + V_\mathrm{ref1}) - V_\mathrm{meas} \\ + 0 &= -T_\mathrm{pord} f_\mathrm{pord} + (1 + s_P\omega)P_\mathrm{ref} - P_\mathrm{ord} \\ + 0 &= -r_\mathrm{pord} + \text{clamp}(f_\mathrm{pord}, R_P^{\min}, R_P^{\max}) +\end{aligned} +``` + +```math +\begin{aligned} + 0 &= -{I_\mathrm{q}^\mathrm{circ}}^2 + (I^{\max})^2 - s_{PQ}(I_\mathrm{p}^\mathrm{cmd})^2 \\ + 0 &= -{I_\mathrm{p}^\mathrm{circ}}^2 + (I^{\max})^2 - s_{PQ}^\mathrm{off}(I_\mathrm{q}^\mathrm{cmd})^2 \\ + 0 &= -I_\mathrm{q}^{\max} + \text{min}(g_q(V_\mathrm{meas}), I_\mathrm{q}^\mathrm{circ}) \\ + 0 &= -I_\mathrm{p}^{\max} + \text{min}(g_p(V_\mathrm{meas}), I_\mathrm{p}^\mathrm{circ}) \\ + 0 &= -I_\mathrm{qbase} + \text{clamp}(K_\mathrm{vp} e_\mathrm{PIV} + x_\mathrm{PIV}, -I_\mathrm{q}^{\max}, I_\mathrm{q}^{\max}) \\ + 0 &= -I_\mathrm{q}^\mathrm{raw} + s_Q I_\mathrm{qbase} + s_Q^\mathrm{off} Q_V + s_\mathrm{dip} I_\mathrm{qv} \\ + 0 &= -I_\mathrm{q}^\mathrm{cmd} + \text{clamp}(I_\mathrm{q}^\mathrm{raw}, -I_\mathrm{q}^{\max}, I_\mathrm{q}^{\max}) \\ + 0 &= -I_\mathrm{p}^\mathrm{cmd} + \text{clamp}(P_\mathrm{ord}/V_\mathrm{meas}^\mathrm{safe}, 0, I_\mathrm{p}^{\max}) +\end{aligned} +``` + +The $V_T$, $I_{\mathrm{q}}^{\mathrm{circ}}$, and $I_{\mathrm{p}}^{\mathrm{circ}}$ variables use nonnegative branches of squared algebraic residuals. This preserves the $s_{PQ}=0$ Q-priority and $s_{PQ}=1$ P-priority current-circle behavior without explicit square roots; a consistent solution should satisfy the nonnegative branch and nonnegative radicands. + +CommonMath defines the helper targets and smooth approximations for [min, max, clamp, deadband, and outside](../../../../CommonMath.md#derived-functions). + +## Initialization + +Initialization is performed by evaluating the steady-state residuals in dependency order. Let subscript $0$ denote initial values and set all internal derivatives to zero. If optional signals are not connected, use steady-state constants: + +```math +\begin{aligned} + V_{T,0} &= \sqrt{V_{\mathrm{r},0}^2 + V_{\mathrm{i},0}^2} \\ + \omega_0 &= 0 \\ + Q_{\mathrm{ext},0} &= Q_{\mathrm{gen},0} \\ + P_{\mathrm{ref},0} &= \dfrac{P_{e,0}}{1+s_P\omega_0} +\end{aligned} +``` + +Connected optional signals use their supplied initial values; if only some are omitted, compute the omitted constants with the connected initial values. Inconsistent supplied commands require a residual solve or initialization rejection. + +If $V_{\mathrm{ref0}}$ is omitted, set $V_{\mathrm{ref0}} = V_{T,0}$. Initialize the measurement variables from the descriptor-form filter residuals: + +```math +\begin{aligned} + V_{\mathrm{meas},0} &= V_{T,0} \\ + P_{\mathrm{meas},0} &= P_{e,0} +\end{aligned} +``` + +When $T_{\mathrm{rv}} = 0$ or $T_{\mathrm{p}} = 0$, the corresponding relation is an algebraic residual rather than a differential-state initial condition. + +Then evaluate the upstream algebraic chain: + +```math +\begin{aligned} + V_{\mathrm{meas},0}^{\mathrm{safe}} &= \text{max}(V_{\mathrm{meas},0}, 0.01) \\ + s_{\mathrm{dip},0} &= \text{outside}(V_{T,0}, V_{\mathrm{dip}}, V_{\mathrm{up}}) \\ + V_{\mathrm{err},0} &= \text{deadband}(V_{\mathrm{ref0}} - V_{\mathrm{meas},0}, D_{\mathrm{bd1}}, D_{\mathrm{bd2}}) \\ + I_{\mathrm{qv},0} &= \text{clamp}(K_{\mathrm{qv}} V_{\mathrm{err},0}, I_{\mathrm{qinj}}^{\min}, I_{\mathrm{qinj}}^{\max}) \\ + Q_{\mathrm{ref},0} &= s_{\mathrm{pf}} P_{\mathrm{meas},0}\tan(\phi_{\mathrm{pf},0}^{\mathrm{ref}}) + s_{\mathrm{pf}}^{\mathrm{off}} Q_{\mathrm{ext},0} \\ + e_{Q,0} &= \text{clamp}(Q_{\mathrm{ref},0}, Q^{\min}, Q^{\max}) - Q_{\mathrm{gen},0} \\ + Q_{V,0} &= \dfrac{Q_{\mathrm{ref},0}}{V_{\mathrm{meas},0}^{\mathrm{safe}}} \\ + P_{\mathrm{ord},0} &= (1+s_P\omega_0)P_{\mathrm{ref},0} +\end{aligned} +``` + +Initialize the reactive-power PI output so its residual and zero-derivative anti-windup condition hold: + +```math +\begin{aligned} + V_{\mathrm{PIQ},0} &= \text{clamp}(K_{\mathrm{qp}} e_{Q,0} + x_{\mathrm{PIQ},0}, V^{\min}, V^{\max}) \\ + e_{\mathrm{PIV},0} &= s_V V_{\mathrm{PIQ},0} + s_V^{\mathrm{off}}(Q_{\mathrm{ref},0} + V_{\mathrm{ref1}}) - V_{\mathrm{meas},0} +\end{aligned} +``` + +For an unsaturated zero-derivative start, require $e_{Q,0}=0$ for $x_{\mathrm{PIQ}}$ and choose or verify $e_{\mathrm{PIV},0}=0$ for $x_{\mathrm{PIV}}$. Then $x_{\mathrm{PIQ},0}=V_{\mathrm{PIQ},0}-K_{\mathrm{qp}}e_{Q,0}$; when $s_V=1$, set $V_{\mathrm{PIQ},0}=V_{\mathrm{meas},0}$, and when $s_V=0$, the supplied $Q_{\mathrm{ref},0}+V_{\mathrm{ref1}}$ must equal $V_{\mathrm{meas},0}$. Saturated initial conditions should be solved against the anti-windup residuals, not forced by this unsaturated formula. + +Finish by evaluating $g_q(V_{\mathrm{meas},0})$, $g_p(V_{\mathrm{meas},0})$, and the current-limit and current-command algebraic residuals in priority order. At the command steps, use the power-flow current targets before final limiting: + +```math +\begin{aligned} + I_{\mathrm{qbase},0}^{\star} &= \dfrac{Q_{\mathrm{gen},0}}{V_{\mathrm{meas},0}^{\mathrm{safe}}} \\ + I_{\mathrm{p},0}^{\star} &= \dfrac{P_{\mathrm{ord},0}}{V_{\mathrm{meas},0}^{\mathrm{safe}}} +\end{aligned} +``` + +For $s_{PQ}=0$, use: + +```math +\begin{aligned} +I_{\mathrm{q},0}^{\mathrm{circ}} +\rightarrow I_{\mathrm{q},0}^{\max} +\rightarrow I_{\mathrm{qbase},0} +\rightarrow I_{\mathrm{q},0}^{\mathrm{raw}} +\rightarrow I_{\mathrm{q},0}^{\mathrm{cmd}} +\rightarrow I_{\mathrm{p},0}^{\mathrm{circ}} +\rightarrow I_{\mathrm{p},0}^{\max} +\rightarrow I_{\mathrm{p},0}^{\mathrm{cmd}} +\end{aligned} +``` + +For $s_{PQ}=1$, use: + +```math +\begin{aligned} +I_{\mathrm{p},0}^{\mathrm{circ}} +\rightarrow I_{\mathrm{p},0}^{\max} +\rightarrow I_{\mathrm{p},0}^{\mathrm{cmd}} +\rightarrow I_{\mathrm{q},0}^{\mathrm{circ}} +\rightarrow I_{\mathrm{q},0}^{\max} +\rightarrow I_{\mathrm{qbase},0} +\rightarrow I_{\mathrm{q},0}^{\mathrm{raw}} +\rightarrow I_{\mathrm{q},0}^{\mathrm{cmd}} +\end{aligned} +``` + +After $I_{\mathrm{q},0}^{\max}$ and $I_{\mathrm{qbase},0}$ are known, initialize the voltage PI state from its output residual; the unsaturated zero-derivative start also requires the $e_{\mathrm{PIV},0}=0$ condition above: + +```math +x_{\mathrm{PIV},0} = I_{\mathrm{qbase},0} - K_{\mathrm{vp}} e_{\mathrm{PIV},0} +``` + +The current-circle variables use the nonnegative branch of the squared algebraic residuals; initialization must reject negative radicands. A standard steady-state initialization assumes $s_{\mathrm{dip},0}=0$. If initialized during voltage-dip or overvoltage logic, $Q_V$, $P_{\mathrm{ord}}$, and the PI histories are not uniquely determined without the unsupported hold-timer histories, so the implementation should solve a saturation-consistent state or reject the start. + +## Model Outputs + +Output | Units | Description | Note +----------------|--------|-------------------------------------|------ +`iqcmd` | [p.u.] | Reactive-current command output | Converter base +`ipcmd` | [p.u.] | Active-current command output | Converter base +`vmeas` | [p.u.] | Filtered terminal voltage | +`pmeas` | [p.u.] | Filtered electrical power | +`piq` | [p.u.] | Reactive-power PI controller state | +`piv` | [p.u.] | Voltage PI controller state | +`qv` | [p.u.] | Reactive-current command lag state | +`pord` | [p.u.] | Filtered active-power order | +`qref` | [p.u.] | Selected reactive-power reference | +`sdip` | [binary] | Voltage-dip/overvoltage freeze indicator | +`iqmax` | [p.u.] | Final reactive-current upper limit | Converter base +`ipmax` | [p.u.] | Final active-current upper limit | Converter base +`iqv` | [p.u.] | Reactive-current injection candidate | Converter base +`vqctrl` | [p.u.] | Reactive-power control PI output | +`iqbase` | [p.u.] | Base reactive-current command | Converter base + +## Outstanding + +Nonzero $T_{\mathrm{hld}}$ and $T_{\mathrm{hld2}}$ require timer/history-state support and are not modeled yet. With the required zero values, $I_{\mathrm{qinj}}^{\mathrm{frz}}$ is unreachable and $I_{\mathrm{p}}^{\max}$ is recalculated from VDL2 and current-circle logic at each residual evaluation instead of held after voltage recovery. + +A future smooth approximation of the held reactive-current path could introduce a continuous gate $h_q$: + +```math +\dot h_q = + \dfrac{1}{T_{\mathrm{rise}}} s_{\mathrm{dip}}(1-h_q) + - \dfrac{1}{T_{\mathrm{hld}}} (1-s_{\mathrm{dip}})h_q +``` + +That approximation is a modeling choice and is not part of the present equations. diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md b/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md new file mode 100644 index 000000000..ea352c520 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md @@ -0,0 +1,311 @@ +# **Renewable Energy Generator/Converter Model (REGCA)** + +REGCA is a first-generation WECC renewable generator/converter model for +inverter-coupled resources. In GridKit it is represented as a controlled +current source at the network interface. + +Notes: +- LVACM uses the unfiltered terminal voltage $V_T$; LVPL uses the filtered voltage $V_M$. +- Internal currents are on converter base; bus injections are converted to system base in the network interface. +- HVRCM is represented by internal algebraic current $I_{\mathrm{q}}^{\mathrm{extra}}$. + +## Block Diagram + +Standard REGCA converter-interface model. + +
+ + + Figure 1: Generator/Converter REGCA model. Figure courtesy of [PowerWorld](https://www.powerworld.com/WebHelp/) +
+ +## Model Parameters + +Symbol | Units | Description | Typical Value | Note +---------------------------------|----------|-------------------------------------------------------|---------------|------ +$P_{\mathrm{0}}$ | [p.u.] | Initial active power injection | | On system base +$Q_{\mathrm{0}}$ | [p.u.] | Initial reactive power injection | | On system base +$S^{\mathrm{conv}}$ | [MVA] | Converter/model power base | TBD | +$T_{\mathrm{g}}$ | [sec] | Converter current-control lag time constant | TBD | +$T_M$ | [sec] | Terminal voltage sensor time constant | TBD | Block name: `Tfltr` +$R_{\mathrm{q}}^{\max}$ | [p.u./s] | Reactive-current recovery positive rate limit | TBD | Block name: `Iqrmax` +$R_{\mathrm{q}}^{\min}$ | [p.u./s] | Reactive-current recovery negative rate limit | TBD | Block name: `Iqrmin` +$R_{\mathrm{p}}^{\max}$ | [p.u./s] | Active-current magnitude recovery rate limit | TBD | Block name: `rrpwr` +$s_L$ | [binary] | LVPL switch | TBD | Block name: `LPVLSW` +$I_{L1}$ | [p.u.] | LVPL upper-current ceiling | TBD | Block name: `Lvpl1` +$V_{L0}$ | [p.u.] | LVPL zero-crossing voltage | TBD | Block name: `xerox` +$V_{L1}$ | [p.u.] | LVPL upper breakpoint voltage | TBD | Block name: `brkpt` +$V_{A0}$ | [p.u.] | LVACM lower breakpoint voltage | TBD | Block name: `Lvpnt0` +$V_{A1}$ | [p.u.] | LVACM upper breakpoint voltage | TBD | Block name: `Lvpnt1` +$V_{\mathrm{hv}}^{\max}$ | [p.u.] | Terminal-voltage ceiling for HV reactive management | TBD | Block name: `Vlim` + +### Parameter Validation + +Implementations should reject or report invalid parameter sets: + +```math +\begin{aligned} + S^{\mathrm{conv}} &> 0 & + T_{\mathrm{g}} &> 0 & + T_M &> 0 \\ + R_{\mathrm{p}}^{\max} &> 0 & + R_{\mathrm{q}}^{\min} &< 0 < R_{\mathrm{q}}^{\max} & + s_L &\in \{0,1\} \\ + I_{L1} &\ge 0 & + 0 &\le V_{L0} < V_{L1} & + 0 &\le V_{A0} < V_{A1} \\ + V_{\mathrm{hv}}^{\max} &> 0 +\end{aligned} +``` + +### Model Derived Parameters + +The smooth active-current bound equations use $M_{\mathrm{p}}$, a numerical +relaxation for inactive $\pm\infty$ rate bounds: + +```math +M_{\mathrm{p}} = 100 R_{\mathrm{p}}^{\max} +``` + +$M_{\mathrm{p}}$ is not a physical REGCA parameter; it should be large enough +that inactive bounds do not bind expected $f_{\mathrm{p}}$ values while staying +moderate enough to keep the smooth clamp well conditioned. + +## Model Variables + +### Internal Variables + +#### Differential + +Symbol | Units | Description | Note +----------------------|--------|---------------------------|------ +$V_M$ | [p.u.] | Filtered terminal voltage | State 3 in Fig. 1 +$I_{\mathrm{q}}$ | [p.u.] | Reactive-current state | State 1 in Fig. 1 before the `-1` block; converter base +$I_{\mathrm{p}}$ | [p.u.] | Active-current state | State 2 in Fig. 1; converter base + +#### Algebraic + +Symbol | Units | Description | Note +---------------------------|--------|-------------------------------------------------------------|------ +$V_T$ | [p.u.] | Terminal voltage magnitude | +$I_{\mathrm{i}}$ | [p.u.] | Injected current, imaginary component on network reference frame | Converter base +$I_{\mathrm{q}}^{\mathrm{extra}}$ | [p.u.] | Extra inductive current from high-voltage reactive current management | Converter base +$I_L$ | [p.u.] | LVPL upper-limit current curve | Function of $V_M$ +$I_{\mathrm{r}}$ | [p.u.] | Injected current, real component on network reference frame | Converter base +$\ell_{\mathrm{p}}$ | [p.u./s] | Smooth active-current lower rate bound | Equivalent to diagram `Rdown` +$u_{\mathrm{p}}$ | [p.u./s] | Smooth active-current upper rate bound | Effective `Rup`; includes LVPL anti-windup when $s_L=1$ + +### External Variables + +#### Differential +None. + +#### Algebraic + +Symbol | Units | Description | Note +--------------------------------|--------|------------------------------------------------------------------|------ +$V_{\mathrm{r}}$ | [p.u.] | Terminal voltage, real component on network reference frame | Owned by bus object +$V_{\mathrm{i}}$ | [p.u.] | Terminal voltage, imaginary component on network reference frame | Owned by bus object +$I_{\mathrm{q}}^{\mathrm{cmd}}$ | [p.u.] | Reactive-current command | Converter base; owned by REEC, constant if no REEC is connected +$I_{\mathrm{p}}^{\mathrm{cmd}}$ | [p.u.] | Active-current command | Converter base; owned by REEC, constant if no REEC is connected + +## Model Equations + +Define the pre-limit current derivatives: + +```math +\begin{aligned} + f_{\mathrm{q}} &= \dfrac{1}{T_{\mathrm{g}}}(I_{\mathrm{q}}^{\mathrm{cmd}} - I_{\mathrm{q}}) \\ + f_{\mathrm{p}} &= \dfrac{1}{T_{\mathrm{g}}}(I_{\mathrm{p}}^{\mathrm{cmd}} - I_{\mathrm{p}}) +\end{aligned} +``` + +### Differential Equations + +The exact state equations are + +```math +\begin{aligned} + 0 &= -T_M \dot V_M - V_M + V_T \\ + 0 &= -\dot I_{\mathrm{q}} + + \begin{cases} + \min(f_{\mathrm{q}}, R_{\mathrm{q}}^{\max}) & Q_{\mathrm{0}} > 0 \\ + \max(f_{\mathrm{q}}, R_{\mathrm{q}}^{\min}) & Q_{\mathrm{0}} \le 0 + \end{cases} \\ + 0 &= -\dot I_{\mathrm{p}} + \text{clamp}(f_{\mathrm{p}}, \ell_{\mathrm{p}}, u_{\mathrm{p}}) +\end{aligned} +``` + +The implemented smooth state equations are + +```math +\begin{aligned} + 0 &= -T_M \dot V_M - V_M + V_T \\ + 0 &= -\dot I_{\mathrm{q}} + + \begin{cases} + f_{\mathrm{q}} - \rho(f_{\mathrm{q}} - R_{\mathrm{q}}^{\max}) + & Q_{\mathrm{0}} > 0 \\ + f_{\mathrm{q}} + \rho(R_{\mathrm{q}}^{\min} - f_{\mathrm{q}}) + & Q_{\mathrm{0}} \le 0 + \end{cases} \\ + 0 &= -\dot I_{\mathrm{p}} + + \ell_{\mathrm{p}} + + \rho(f_{\mathrm{p}} - \ell_{\mathrm{p}}) + - \rho(f_{\mathrm{p}} - u_{\mathrm{p}}) +\end{aligned} +``` + +Here $\rho$ is GridKit's smooth ramp function. The $I_{\mathrm{q}}$ branch is +selected by initial reactive power $Q_{\mathrm{0}}$. The $I_{\mathrm{p}}$ +equation is the smooth clamp of $f_{\mathrm{p}}$ between the algebraic bounds +$\ell_{\mathrm{p}}$ and $u_{\mathrm{p}}$. + +### Algebraic Equations + +The piecewise definitions in this section switch on continuous states, unlike +the $I_{\mathrm{q}}$ differential branch selected by initial conditions. The +exact algebraic targets are: + +```math +\begin{aligned} + 0 &= -V_T^2 + V_{\mathrm{r}}^2 + V_{\mathrm{i}}^2 \\ + I_{\mathrm{i}} &= -I_{\mathrm{q}} + I_{\mathrm{q}}^{\mathrm{extra}} \\ + 0 &= + \begin{cases} + I_{\mathrm{q}}^{\mathrm{extra}} + & V_T < V_{\mathrm{hv}}^{\max} \\ + V_T - V_{\mathrm{hv}}^{\max} + & I_{\mathrm{q}}^{\mathrm{extra}} > 0 + \end{cases} \\ + I_L &= I_{L1} + \begin{cases} + 0 & V_M \le V_{L0} \\ + \dfrac{V_M - V_{L0}}{V_{L1} - V_{L0}} & V_{L0} < V_M < V_{L1} \\ + 1 & V_M \ge V_{L1} + \end{cases} \\ + I_{\mathrm{r}} &= I_{\mathrm{p}} + \begin{cases} + 0 & V_T \le V_{A0} \\ + \dfrac{V_T - V_{A0}}{V_{A1} - V_{A0}} & V_{A0} < V_T < V_{A1} \\ + 1 & V_T \ge V_{A1} + \end{cases} \\ + \ell_{\mathrm{p}} &= + \begin{cases} + -R_{\mathrm{p}}^{\max} & I_{\mathrm{p}} \le 0 \\ + -\infty & I_{\mathrm{p}} > 0 + \end{cases} \\ + u_{\mathrm{p}} &= + \begin{cases} + R_{\mathrm{p}}^{\max} & I_{\mathrm{p}} \ge 0 \ \land\ (s_L = 0 \lor I_{\mathrm{p}} < I_L) \\ + 0 & s_L = 1 \ \land\ I_{\mathrm{p}} \ge I_L \\ + \infty & I_{\mathrm{p}} < 0 + \end{cases} +\end{aligned} +``` + +The implemented algebraic residuals use smooth $\text{linseg}$, +$\rho$, and $\sigma$ operators: + +```math +\begin{aligned} + 0 &= -V_T^2 + V_{\mathrm{r}}^2 + V_{\mathrm{i}}^2 \\ + 0 &= -I_{\mathrm{i}} - I_{\mathrm{q}} + I_{\mathrm{q}}^{\mathrm{extra}} \\ + 0 &= -I_{\mathrm{q}}^{\mathrm{extra}} + + \rho\!\left(I_{\mathrm{q}}^{\mathrm{extra}} - (V_{\mathrm{hv}}^{\max} - V_T)\right) \\ + 0 &= -I_L + + \text{linseg}(V_M;\ V_{L0},\ V_{L1},\ I_{L1}) \\ + 0 &= -I_{\mathrm{r}} + + I_{\mathrm{p}}\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1) \\ + 0 &= -\ell_{\mathrm{p}} + - R_{\mathrm{p}}^{\max} + - (M_{\mathrm{p}} - R_{\mathrm{p}}^{\max})\sigma(I_{\mathrm{p}}) \\ + 0 &= -u_{\mathrm{p}} + + \begin{cases} + M_{\mathrm{p}}(1-\sigma(I_{\mathrm{p}})) + + R_{\mathrm{p}}^{\max}\sigma(I_{\mathrm{p}}) + \sigma(I_L - I_{\mathrm{p}}) + & s_L = 1 \\ + M_{\mathrm{p}}(1-\sigma(I_{\mathrm{p}})) + + R_{\mathrm{p}}^{\max}\sigma(I_{\mathrm{p}}) + & s_L = 0 + \end{cases} +\end{aligned} +``` + +The $V_T$ residual is kept in squared form for smoothness at the origin. + +## Network Interface + +The bus receives system-base current injections converted from converter-base +REGCA currents: + +```math +\begin{aligned} + I_{\mathrm{r}}^{\mathrm{inj}} &:= I_{\mathrm{r}}\dfrac{S^{\mathrm{conv}}}{S^{\mathrm{sys}}} \\ + I_{\mathrm{i}}^{\mathrm{inj}} &:= I_{\mathrm{i}}\dfrac{S^{\mathrm{conv}}}{S^{\mathrm{sys}}} +\end{aligned} +``` + +Positive current injection is into the bus. + +## Initialization + +Given initialized bus voltage $V_{\mathrm{r}}, V_{\mathrm{i}}$, compute the +steady-state initial values: + +```math +\begin{aligned} + V_T &= \sqrt{V_\mathrm{r}^2 + V_\mathrm{i}^2} \\ + I_\mathrm{r0} &= \dfrac{P_0 V_\mathrm{r} + Q_0 V_\mathrm{i}}{V_T^2} + \dfrac{S^\mathrm{sys}}{S^\mathrm{conv}} \\ + I_\mathrm{i0} &= \dfrac{P_0 V_\mathrm{i} - Q_0 V_\mathrm{r}}{V_T^2} + \dfrac{S^\mathrm{sys}}{S^\mathrm{conv}} \\ + V_{M0} &= V_T \\ + I_{L0} &= \text{linseg}(V_T;\ V_{L0},\ V_{L1},\ I_{L1}) \\ + I_\mathrm{p0} &= \dfrac{I_\mathrm{r0}} + {\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1)} \\ + \ell_\mathrm{p0} &= -R_\mathrm{p}^{\max} + - (M_\mathrm{p} - R_\mathrm{p}^{\max})\sigma(I_\mathrm{p0}) \\ + u_\mathrm{p0} &= + \begin{cases} + M_\mathrm{p}(1-\sigma(I_\mathrm{p0})) + + R_\mathrm{p}^{\max}\sigma(I_\mathrm{p0}) + \sigma(I_{L0} - I_\mathrm{p0}) + & s_L = 1 \\ + M_\mathrm{p}(1-\sigma(I_\mathrm{p0})) + + R_\mathrm{p}^{\max}\sigma(I_\mathrm{p0}) + & s_L = 0 + \end{cases} \\ + I_\mathrm{q0}^\mathrm{extra} &= 0 \\ + I_\mathrm{q0} &= -I_\mathrm{i0} + I_\mathrm{q0}^\mathrm{extra} \\ + I_\mathrm{p0}^\mathrm{cmd} &= I_\mathrm{p0} \\ + I_\mathrm{q0}^\mathrm{cmd} &= I_\mathrm{q0} +\end{aligned} +``` + +For normal power-flow starts, $V_T > V_{A1}$, so +$\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1) = 1$ and the +$I_{\mathrm{p0}}$ formula is well defined. + +Initialization should verify: +- $V_T \le V_{\mathrm{hv}}^{\max}$. If $V_T \ge V_{\mathrm{hv}}^{\max}$, + $I_{\mathrm{q0}}^{\mathrm{extra}} = 0$ may not satisfy the HVRCM algebraic + condition, and a nonzero value should be solved or the initialization rejected. +- $\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1) > 0$ when + $I_{\mathrm{r0}} \ne 0$. If the LVACM gain is zero, no finite + $I_{\mathrm{p0}}$ can reproduce nonzero initial active current. + +All internal derivatives initialize to zero. + +## Model Outputs + +Real and imaginary injected currents, $I_{\mathrm{r}}$ and +$I_{\mathrm{i}}$, are converter-base algebraic variables. System-base power +outputs use the bus-facing currents: +```math +\begin{aligned} + P &= V_{\mathrm{r}} I_{\mathrm{r}}^{\mathrm{inj}} + V_{\mathrm{i}} I_{\mathrm{i}}^{\mathrm{inj}} \\ + Q &= V_{\mathrm{i}} I_{\mathrm{r}}^{\mathrm{inj}} - V_{\mathrm{r}} I_{\mathrm{i}}^{\mathrm{inj}} +\end{aligned} +``` +Power outputs are positive leaving the converter and entering the bus. diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCB/README.md b/GridKit/Model/PhasorDynamics/Converter/REGCB/README.md new file mode 100644 index 000000000..f10fe2218 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCB/README.md @@ -0,0 +1,117 @@ +# **Renewable Energy Generator/Converter Model (REGCB)** + +REGCB is a WECC renewable energy generator/converter model for inverter-coupled +resources. This document is a skeleton for the model specification; parameters, +equations, initialization details, and default values must be validated against +the REGCB source standard before implementation. + +## Block Diagram + +Standard model diagram for the REGCB converter interface. + +
+ + + Figure 1: Generator/Converter REGCB model. Figure courtesy of [PowerWorld](https://www.powerworld.com/WebHelp/) +
+ +Detailed REGCB parameters, variables, equations, initialization details, and +outputs will be added after validation against the REGCB source standard. + + diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp index e2fa7dab5..f535862f0 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp @@ -231,12 +231,11 @@ namespace GridKit ScalarT vimag = bus_->Vi(); ScalarT Ec = std::sqrt(vreal * vreal + vimag * vimag); - ScalarT efdp = efd0 / (ONE + omega * Ispdlim_); - ScalarT efd_sat = (efdp - SA_) * Math::sigmoid(efdp - SA_); - ScalarT ksat = SB_ * efd_sat * efd_sat; - ScalarT ve = ksat * efdp; - ScalarT vr = Ke_ * efdp + ve; - ScalarT vtr = vr / Ka_; + 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_; ScalarT vf{0}; ScalarT vfx = (Kf_ / Tf_) * efdp; @@ -325,7 +324,7 @@ namespace GridKit // The 'pre-limit' derivative of Vr. ScalarT func = (-vr + Ka_ * vtr) / Ta_; ScalarT func_normalized = func / static_cast(500.0); // TODO This is arbitrary, need more general conditioning method that is fast - ScalarT vr_ind = Math::indicator(Vrmin_, Vrmax_, vr, func_normalized); + ScalarT vr_ind = Math::indicator(vr, func_normalized, Vrmin_, Vrmax_); // Internal Differential Equations f[0] = -vts_dot + (Ec - vts) / Tr_; @@ -339,8 +338,7 @@ namespace GridKit f[6] = -ve + ksat * efdp; f[7] = -efd + efdp + omega * efdp * Ispdlim_; - ScalarT efd_sat = (efdp - SA_) * (Math::sigmoid(efdp - SA_)); - f[8] = -ksat + SB_ * efd_sat * efd_sat; + f[8] = -ksat + SB_ * Math::qramp(efdp - SA_); return 0; } @@ -471,7 +469,7 @@ namespace GridKit monitor_->set(Variable::efd, [this] { return efd_signal_->read(); }); monitor_->set(Variable::ksat, [this] - { return SB_ * (((y_[2] - SA_) * Math::sigmoid(y_[2] - SA_)) * ((y_[2] - SA_) * Math::sigmoid(y_[2] - SA_))); }); + { return SB_ * Math::qramp(y_[2] - SA_); }); } } // namespace Exciter } // namespace PhasorDynamics diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md index 7b2c6e358..9bb2ef5fb 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md @@ -4,10 +4,15 @@ Standard model of the IEEET1 Exciter. +Notes: +- $E_C$ should be an external signal from the generator or machine, not computed through an exciter-owned `Bus` reference. +- The current implementation uses its `Bus` reference as a proxy for $E_C$. +- This direct coupling affects numerical conditioning; production models typically use a decoupling reactance for the exciter-current path that forms $E_C$. +
- - + + Figure 1: Exciter IEEET1 model. Figure courtesy of [PowerWorld](https://www.powerworld.com/WebHelp/)
@@ -17,19 +22,19 @@ Standard model of the IEEET1 Exciter. Symbol | Units | Description | Typical Value | Note ------------|--------|--------------------------------------|---------| ------ $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 | -$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 | -$V_{rmin}$ | [p.u.] | Lower limit to voltage regulation | -1 | -$V_{rmax}$ | [p.u.] | Upper limit to voltage regulation | 1 | -$E_1$ | [p.u.] | Saturation Parameter | 2.8 | -$E_2$ | [p.u.] | Saturation Parameter | 3.73 | -$S_{e1}$ | [p.u.] | Saturation Parameter | 0.04 | -$S_{e2}$ | [p.u.] | Saturation Parameter | 0.33 | -$I_{spdlm}$ | [binary] | Speed Limit flag indicator | 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 | +$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 | +$V_R^{\min}$ | [p.u.] | Lower limit to voltage regulation | -1 | +$V_R^{\max}$ | [p.u.] | Upper limit to voltage regulation | 1 | +$E_1$ | [p.u.] | Saturation Parameter | 2.8 | +$E_2$ | [p.u.] | Saturation Parameter | 3.73 | +$S_1$ | [p.u.] | Saturation Parameter | 0.04 | +$S_2$ | [p.u.] | Saturation Parameter | 0.33 | +$I_\text{spdlm}$ | [binary] | Speed limit flag indicator | 0 | ### Model Derived Parameters @@ -37,8 +42,8 @@ The relationship of the derived parameters is defined by the following quadratic the expected saturation near the operating region. ``` math \begin{aligned} - S_{e1} &= S_B(E_1-S_A)^2 \\ - S_{e2} &= S_B(E_2-S_A)^2 \\ + S_1 &= S_B(E_1-S_A)^2 \\ + S_2 &= S_B(E_2-S_A)^2 \\ \end{aligned} ``` @@ -47,23 +52,23 @@ Generally, this system has two solutions. The non-extraneous solution is as foll \begin{aligned} C &= \sqrt{ \dfrac - {S_{e2}} - {S_{e1}} - } + {S_2} + {S_1} + } \\ - S_A &= + S_A &= \dfrac {C E_1 - E_2} - {C - 1} + {C - 1} \\ - S_B &= + S_B &= \dfrac - {S_{e1}} + {S_1} {(E_1-S_A)^2} \end{aligned} ``` -## Model Variables +## Model Variables ### Internal Variables @@ -72,9 +77,9 @@ Generally, this system has two solutions. The non-extraneous solution is as foll Symbol | Units | Description | Note ----------|--------|-----------------------------------|------- $V_{ts}$ | [p.u.] | Sensed terminal voltage | -$V_{R}$ | [p.u.] | Voltage regulator | -$E_{fd}'$ | [p.u.] | Field-current pre-speed multiplier| -$V_{fx}$ | [p.u.] | Exciter feedback internal state | +$V_R$ | [p.u.] | Voltage regulator | +$E_{fd}'$ | [p.u.] | Field-current pre-speed multiplier| +$V_{fx}$ | [p.u.] | Exciter feedback internal state | #### Algebraic @@ -83,10 +88,10 @@ $V_{fx}$ | [p.u.] | Exciter feedback internal state | Symbol | Units | Description | Note ----------------|--------|-----------------------------------|------- $V_{tr}$ | [p.u.] | Terminal Voltage Error | -$V_{f}$ | [p.u.] | Feedback Voltage | -$V_{E}$ | [p.u.] | Excitation control voltage | +$V_f$ | [p.u.] | Feedback Voltage | +$V_E$ | [p.u.] | Excitation control voltage | $E_{fd}$ | [p.u.] | Field winding voltage | -$k_{sat}$ | [p.u.] | Saturation variable | +$k_\text{sat}$ | [p.u.] | Saturation variable | ### External Variables @@ -100,11 +105,11 @@ $\omega$ | [p.u.] | Machine Speed Deviation | Read from a Mac Symbol | Units | Description | Note ----------------|--------|-----------------------------------|------- -$E_{C}$ | [p.u.] | Compensated machine terminal voltage magnitude | -$V_{ref}$ | [p.u.] | Reference terminal voltage | +$E_C$ | [p.u.] | Compensated machine terminal voltage magnitude | +$V_\text{ref}$ | [p.u.] | Reference terminal voltage | $V_{UEL}$ | [p.u.] | Input from under excitation limiter | $V_{OEL}$ | [p.u.] | Input from over excitation limiter | -$V_{S}$ | [p.u.] | Input from stabilizer controller | +$V_S$ | [p.u.] | Input from stabilizer controller | ## Model Equations @@ -114,76 +119,70 @@ $V_{S}$ | [p.u.] | Input from stabilizer controller | The IEEET1 differential equations, as derived from the model diagram. Define the pre-limit derivative of $V_R$ ```math -f = \dfrac{1}{T_A}\left[-V_R + K_a V_{tr}\right] +f = \dfrac{1}{T_A}\left[-V_R + K_A V_{tr}\right] ``` -so that $\dot V_R$ can be written in piecewise form compactly. +so that $\dot V_R$ is the anti-windup limited derivative. ```math \begin{aligned} - \dot{V}_{ts} &= \dfrac{1}{T_R}(E_C-V_{ts}) \\ - \dot{V}_{R} &= - \begin{cases} - f - & \text{if } (V_{rmin} < V_R < V_{rmax}) & \lor \\ - & \quad (V_R \leq V_{rmin} \land f>0) & \lor \\ - & \quad(V_R \geq V_{rmax} \land f<0) \\ - 0 & \text{else } \\ - \end{cases} \\ - \dot{E}_{fd}' &= \dfrac{1}{T_E}(V_R-V_E-K_E E_{fd}') \\ - \dot{V}_{fx} &= \dfrac{1}{T_F}V_f \\ + \dot V_{ts} &= \dfrac{1}{T_R}(E_C-V_{ts}) \\ + \dot V_R &= \text{antiwindup} + \left(V_R, f;\ V_R^{\min}, V_R^{\max}\right) \\ + \dot E_{fd}' &= \dfrac{1}{T_E}(V_R-V_E-K_E E_{fd}') \\ + \dot V_{fx} &= \dfrac{1}{T_F}V_f \\ \end{aligned} ``` -In simulation the piecewise form above is replaced with a smooth approximation where $\phi$ is GridKit's smooth anti-windup indicator. See [CommonMath: Anti-Windup Indicator](../../../../CommonMath.md#anti-windup-indicator) for its definition, behavior, and design rationale. +CommonMath defines the [Anti-Windup](../../../../CommonMath.md#anti-windup-indicator) target and smooth approximation. ### Algebraic Equations The algebraic equations of the exciter. ```math \begin{aligned} - V_{tr} &= V_{ref} +V_{UEL} + V_{OEL} + V_S - V_f- V_{ts}\\ - V_{f} &= \dfrac{E_{fd}' K_F}{T_F} - V_{fx}\\ - V_{E} &= k_{sat}\cdot E_{fd}' \\ + V_{tr} &= V_\text{ref} +V_{UEL} + V_{OEL} + V_S - V_f- V_{ts}\\ + V_f &= \dfrac{E_{fd}' K_F}{T_F} - V_{fx}\\ + V_E &= k_\text{sat}\cdot E_{fd}' \\ E_{fd}&= \begin{cases} - (1+\omega)E_{fd}' & \text{if } I_{spdlm}=1\\ + (1+\omega)E_{fd}' & \text{if } I_\text{spdlm}=1\\ E_{fd}' & \text{else} \\ \end{cases}\\ - k_{sat}&= \begin{cases} + k_\text{sat}&= \begin{cases} S_B(E_{fd}' -S_A)^2 & \text{if } E_{fd}' >S_A\\ 0 & \text{else } \\ \end{cases} \\ \end{aligned} ``` -#### Smooth Piecewise Approximation (Algebraic) +#### Smooth Piecewise Approximation (Algebraic) For the algebraic piecewise functions (non-flags), this implementation is straightforward when the approximation above is used. +Here $q$ is GridKit's [Quadratic Ramp](../../../../CommonMath.md#primitives). ```math \begin{aligned} E_{fd} - &=(1 + \omega I_{spdlm})E_{fd}' \\ - k_{sat} - &=S_B\left[(E_{fd}' -S_A) \cdot \sigma (E_{fd}' -S_A)\right]^2 + &=(1 + \omega I_\text{spdlm})E_{fd}' \\ + k_\text{sat} + &=S_B\, q(E_{fd}' -S_A) \end{aligned} ``` -The approximation approaches an exact solution as $\alpha\to\infty$. ## Initialization -The machine initializes $E_{fd}$ first. IEEET1 reads that value as $E_{fd,0}$, along with any attached $\omega$ and $V_S$, and solves the steady-state algebraic chain so all residuals vanish with $\dot y = 0$. There is no compensation impedance, so $E_C$ is taken as the terminal voltage magnitude. Saturation and the speed-limit flag are included directly; $V_{ref}$ is set to close the $V_{tr}$ equation with the current auxiliary inputs. +The machine initializes $E_{fd}$ first. IEEET1 reads that value as $E_{fd,0}$, along with any attached $\omega$ and $V_S$, and solves the steady-state algebraic chain so all residuals vanish with $\dot y = 0$. There is no compensation impedance, so $E_C$ is taken as the terminal voltage magnitude. Saturation and the speed-limit flag are included directly; $V_\text{ref}$ is set to close the $V_{tr}$ equation with the current auxiliary inputs. ```math \begin{aligned} E_C &= \sqrt{V_r^2 + V_i^2} \\ - E_{fd}' &= \dfrac{E_{fd,0}}{1 + I_{spdlm}\,\omega} \\ - k_{sat} &= S_B\big[(E_{fd}' - S_A)\,\sigma(E_{fd}' - S_A)\big]^2 \\ - V_E &= k_{sat}\, E_{fd}' \\ + E_{fd}' &= \dfrac{E_{fd,0}}{1 + I_\text{spdlm}\,\omega} \\ + k_\text{sat} &= S_B\, q(E_{fd}' - S_A) \\ + V_E &= k_\text{sat}\, E_{fd}' \\ V_R &= K_E\, E_{fd}' + V_E \\ - V_{tr} &= \dfrac{V_R}{K_a} \\ + V_{tr} &= \dfrac{V_R}{K_A} \\ V_{fx} &= \dfrac{K_F}{T_F}\, E_{fd}' \\ V_{ts} &= E_C \\ V_f &= 0 \\ - V_{ref} &= E_C + V_{tr} - V_{UEL} - V_{OEL} - V_S + V_\text{ref} &= E_C + V_{tr} - V_{UEL} - V_{OEL} - V_S \end{aligned} ``` @@ -192,4 +191,4 @@ All internal derivatives initialize to zero. The field voltage, $E_{fd}$, is an internal model variable. -The magnetic saturation coefficient $k_{sat}$ is calculated from $E_{fd}$ using the the smooth piecewise version (above). +The magnetic saturation coefficient $k_\text{sat}$ is calculated from $E_{fd}$ using the smooth piecewise version above. diff --git a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/README.md b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/README.md index 5bdddef7c..5ee75a980 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/README.md +++ b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/README.md @@ -18,8 +18,8 @@ $T_A$ | [sec] | Numerator time constant of lag-lead block | $T_B$ | [sec] | Denominator time constant of lag-lead block | | $T_E$ | [sec] | Exciter field time constant | | $K$ | [p.u.] | Voltage regulator gain | | -$E_{fd,\max}$ | [p.u.] | Maximum excitation output | | -$E_{fd,\min}$ | [p.u.] | Minimum excitation output | | +$E_{fd}^{\max}$ | [p.u.] | Maximum excitation output | | +$E_{fd}^{\min}$ | [p.u.] | Minimum excitation output | | PowerWorld/PSS/E SEXS_PTI data often gives $T_A/T_B$ as a ratio. GridKit stores $T_A$ and $T_B$ separately, so convert ratio-format data with @@ -76,9 +76,9 @@ so that $\dot E_{fd}$ can be written in piecewise form compactly. \dot E_{fd} &= \begin{cases} f - & \text{if } (E_{fd,\min} < E_{fd} < E_{fd,\max}) & \lor \\ - & \quad (E_{fd} \leq E_{fd,\min} \land f > 0) & \lor \\ - & \quad (E_{fd} \geq E_{fd,\max} \land f < 0) \\ + & \text{if } (E_{fd}^{\min} < E_{fd} < E_{fd}^{\max}) & \lor \\ + & \quad (E_{fd} \leq E_{fd}^{\min} \land f > 0) & \lor \\ + & \quad (E_{fd} \geq E_{fd}^{\max} \land f < 0) \\ 0 & \text{else} \end{cases} \end{aligned} diff --git a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp index ceebbac4a..f09fd3b5f 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/SEXS-PTI/SexsPtiImpl.hpp @@ -199,7 +199,7 @@ namespace GridKit ScalarT func = (-efd + (K_ / Tb_) * (-vr + Ta_ * vtr)) / Te_; ScalarT func_normalized = func / static_cast(100.0); // TODO arbitrary conditioning; see IEEET1 - ScalarT efd_ind = Math::indicator(Efdmin_, Efdmax_, efd, func_normalized); + ScalarT efd_ind = Math::indicator(efd, func_normalized, Efdmin_, Efdmax_); f[0] = -vr_dot + (-vr + Ta_ * vtr) / Tb_ - vtr; f[1] = -efd_dot + efd_ind * func; diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/README.md b/GridKit/Model/PhasorDynamics/Governor/Tgov1/README.md index d7b778bb1..c51dc72b1 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/README.md +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/README.md @@ -18,8 +18,8 @@ $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_{vmax}$ | [p.u.] | Max Valve Position | 1 | -$P_{vmin}$ | [p.u.] | Min Valve Position | 0 | +$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 | ## Model Variables @@ -31,12 +31,12 @@ $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_v$ | [p.u.] | Valve Position (State 2 in Fig. 1)| #### Algebraic Symbol | Units | Description | Note ----------------|--------|-----------------------------------|------- -$P_{m}$ | [p.u.] | Mechnical Power to Generator | Read by a Machine Model +$P_m$ | [p.u.] | Mechnical Power to Generator | Read by a Machine Model ### External Variables @@ -62,13 +62,13 @@ f = \dfrac{1}{T_1}\left[-P_v + \dfrac{1}{R}(P_{ref} - \omega)\right] so that $\dot P_v$ can be written in piecewise form compactly. ```math \begin{aligned} - \dot{P}_{tx} &= P_v - \dfrac{1}{T_3}(P_{tx}+T_2P_v) \\ - \dot{P}_{v} &= + \dot P_{tx} &= P_v - \dfrac{1}{T_3}(P_{tx}+T_2P_v) \\ + \dot P_v &= \begin{cases} f - & \text{if } (P_{vmin} < P_v < P_{vmax}) & \lor \\ - & \quad (P_v \leq P_{vmin} \land f>0) & \lor \\ - & \quad(P_v \geq P_{vmax} \land f<0) \\ + & \text{if } (P_v^{\min} < P_v < P_v^{\max}) & \lor \\ + & \quad (P_v \leq P_v^{\min} \land f>0) & \lor \\ + & \quad(P_v \geq P_v^{\max} \land f<0) \\ 0 & \text{else} \end{cases} @@ -79,26 +79,26 @@ so that $\dot P_v$ can be written in piecewise form compactly. The algebraic equation dictating the mechnical power output. ```math \begin{aligned} - P_{m} &= \dfrac{1}{T_3}(P_{tx}+T_2P_v) - D_t \omega \\ + P_m &= \dfrac{1}{T_3}(P_{tx}+T_2P_v) - D_t \omega \\ \end{aligned} ``` In simulation the piecewise form above is replaced with a smooth approximation where $\phi$ is GridKit's smooth anti-windup indicator. See [CommonMath: Anti-Windup Indicator](../../../../CommonMath.md#anti-windup-indicator) for its definition, behavior, and design rationale. ## Initialization -At steady state we assume that $P_v$ is at or within its limits. This implies the initial conditions are a function of $P_{m}$ which is equal to the electric torque. +At steady state we assume that $P_v$ is at or within its limits. This implies the initial conditions are a function of $P_m$ which is equal to the electric torque. ```math \begin{aligned} - P_{tx} &= (T_3-T_2) P_{m}\\ - P_{v} &= P_{m}\\ - \dot{P}_{tx} &=0\\ - \dot{P}_{v} &=0\\ + P_{tx} &= (T_3-T_2) P_m\\ + P_v &= P_m\\ + \dot P_{tx} &=0\\ + \dot P_v &=0\\ \end{aligned} ``` And if the reference power is a constant parameter, we can determine the value by solving the steady state equations. ```math \begin{aligned} - P_{ref} &= R P_{m}\\ + P_{ref} &= R P_m\\ \end{aligned} ``` diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp index 2d7c7d87a..836e9c3ff 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp @@ -257,12 +257,11 @@ namespace GridKit ScalarT omega = ws[0]; // The 'pre-limit' derivative of Pv - ScalarT func = (-pv + (pref_ - omega) / R_) / T1_; - ScalarT valv_ind = Math::indicator(Pvmin_, Pvmax_, pv, func); + ScalarT func = (-pv + (pref_ - omega) / R_) / T1_; // Internal Differential Equations f[0] = -ptx_dot + pv - (ptx + T2_ * pv) / T3_; - f[1] = -pv_dot + valv_ind * func; + f[1] = -pv_dot + Math::antiwindup(pv, func, Pvmin_, Pvmax_); // Internal Algebraic Equations f[2] = -pmech + (ptx + T2_ * pv) / T3_ - (Dt_ * omega); diff --git a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md index 06dccb66c..72d732a0d 100644 --- a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md +++ b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md @@ -147,7 +147,7 @@ are specified: `Tgov1 ` | the TGOV1 governor model | `pmech`, `speed` | `R`, `T1`, `T2`, `T3`, `Pvmax`, `Pvmin`, `Dt` | `none` `Ieeet1` | the IEEET1 exciter model | `bus`, `speed`, `efd`, `vs`\* | `Tr`, `Ka`, `Ta`, `Ke`, `Te`, `Kf`, `Tf`, `Vrmin`, `Vrmax`, `E1`, `E2`, `Se1`, `Se2`, `Ispdlim` | `efd`, `ksat` `SexsPti` | the SEXS-PTI simplified exciter model | `bus`, `efd`, `vs`\* | `Ta`, `Tb`, `Te`, `K`, `Efdmax`, `Efdmin` | `efd` - `Ieeest` | the IEEEST stabilizer model | `input`, `output`, `cutout`\* | `A1`, `A2`, `A3`, `A4`, `A5`, `A6`, `T1`, `T2`, `T3`, `T4`, `T5`, `T6`, `Ks`, `Lsmin`, `Lsmax`, `Vcl`, `Vcu`, `Tdelay` | `vs` + `Ieeest` | the IEEEST stabilizer model | `input`, `output` | `A1`, `A2`, `A3`, `A4`, `A5`, `A6`, `T1`, `T2`, `T3`, `T4`, `T5`, `T6`, `Ks`, `Lsmin`, `Lsmax`, `Vcl`, `Vcu`, `Tdelay` | `vss` `BusFault` | simple impedance-based fault at a bus | `bus`, `status`\* | `state0`, `R`, `X` | `state`, `ir`, `ii` Ports marked with \* are optional and, if missing, will be assumed to be diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp index b495a49d3..10a137d2f 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/Ieeest.hpp @@ -46,16 +46,14 @@ namespace GridKit V5, ///< Lead-lag 1 output V6, ///< Lead-lag 2 output V7, ///< Unlimited stabilizer signal - VSS, ///< Limited stabilizer signal - VS, ///< Stabilizer output + VSS, ///< Limited stabilizer signal (model output) MAXIMUM, }; /// External variables of a `Ieeest` enum class IeeestExternalVariables : size_t { - U, ///< Stabilizer input signal - VCT, ///< Cutout signal + U, ///< Stabilizer input signal MAXIMUM, }; @@ -134,7 +132,7 @@ namespace GridKit RealT Lsmin_{-0.1}; RealT Lsmax_{0.1}; RealT Vcl_{0}; - RealT Vcu_{1.5}; + RealT Vcu_{0}; RealT Tdelay_{0}; RealT a0_{1}; @@ -154,15 +152,10 @@ namespace GridKit RealT safe_inv_a2_{0}; RealT use_T2_block_{1}; RealT bypass_T2_block_{0}; - RealT safe_inv_T2_{1}; RealT use_T4_block_{1}; RealT bypass_T4_block_{0}; - RealT safe_inv_T4_{1}; RealT use_T6_block_{1}; RealT bypass_T6_block_{0}; - RealT safe_inv_T6_{1}; - RealT use_cutout_{1}; - RealT bypass_cutout_{0}; ComponentSignals signals_; diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestData.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestData.hpp index 00d728a58..49b241832 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestData.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestData.hpp @@ -34,8 +34,8 @@ namespace GridKit Ks, ///< Stabilizer gain Lsmin, ///< Minimum stabilizer output limit Lsmax, ///< Maximum stabilizer output limit - Vcl, ///< Lower input cutout threshold - Vcu, ///< Upper input cutout threshold + Vcl, ///< Lower input cutout threshold (not modeled) + Vcu, ///< Upper input cutout threshold (not modeled) Tdelay, ///< Input time delay (not modeled) }; @@ -45,7 +45,6 @@ namespace GridKit enum class IeeestPorts { input, ///< Unique ID of the stabilizer input signal - cutout, ///< Unique ID of the cutout signal output, ///< Unique ID of the stabilizer output signal }; @@ -54,7 +53,7 @@ namespace GridKit */ enum class IeeestMonitorableVariables { - vs, ///< Stabilizer output + vss, ///< Stabilizer output (limited signal) }; /** diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp index b0c815d90..e428c601a 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/IeeestImpl.hpp @@ -26,7 +26,7 @@ namespace GridKit template Ieeest::Ieeest() { - size_ = 13; + size_ = 12; } template @@ -35,7 +35,7 @@ namespace GridKit { initializeParameters(data); initializeMonitor(); - size_ = 13; + size_ = 12; } template @@ -138,19 +138,12 @@ namespace GridKit use_T2_block_ = static_cast(T2_ != 0.0); bypass_T2_block_ = 1.0 - use_T2_block_; - safe_inv_T2_ = use_T2_block_ / (T2_ + bypass_T2_block_); use_T4_block_ = static_cast(T4_ != 0.0); bypass_T4_block_ = 1.0 - use_T4_block_; - safe_inv_T4_ = use_T4_block_ / (T4_ + bypass_T4_block_); use_T6_block_ = static_cast(T6_ != 0.0); bypass_T6_block_ = 1.0 - use_T6_block_; - safe_inv_T6_ = use_T6_block_ / (T6_ + bypass_T6_block_); - - // Vcl = Vcu = 0 means "cutout disabled" (passthrough: vs = vss). - use_cutout_ = static_cast(Vcl_ != 0.0 || Vcu_ != 0.0); - bypass_cutout_ = 1.0 - use_cutout_; } template @@ -171,12 +164,10 @@ namespace GridKit variable_indices_.resize(size); residual_indices_.resize(size); - ws_.resize(2); - ws_indices_.resize(2); + ws_.resize(1); + ws_indices_.resize(1); ws_[0] = 0.0; ws_indices_[0] = INVALID_INDEX; - ws_[1] = 0.0; - ws_indices_[1] = INVALID_INDEX; for (IdxT j = 0; j < size_; ++j) { @@ -184,10 +175,10 @@ namespace GridKit this->setResidualIndex(j, j); } - if (signals_.template isAssigned()) + if (signals_.template isAssigned()) { - signals_.template getSignalNode()->set( - &y_[12], &(this->getVariableIndex(12))); + signals_.template getSignalNode()->set( + &y_[11], &(this->getVariableIndex(11))); } return 0; @@ -212,15 +203,6 @@ namespace GridKit ret += 1; } - if (signals_.template isAttached()) - { - if (!signals_.template isLinked()) - { - Log::error() << "Ieeest: cutout signal VCT attached with no linked source\n"; - ret += 1; - } - } - if (a4_ == 0 && a3_ == 0 && a2_ == 0 && a1_ != 0) { Log::error() << "Ieeest: a2, a3, and a4 are all zero - no valid notch filter\n"; @@ -249,15 +231,14 @@ namespace GridKit tag_[1] = true; tag_[2] = true; tag_[3] = true; - tag_[4] = true; - tag_[5] = true; - tag_[6] = true; + tag_[4] = (T2_ != 0.0); + tag_[5] = (T4_ != 0.0); + tag_[6] = (T6_ != 0.0); tag_[7] = false; tag_[8] = false; tag_[9] = false; tag_[10] = false; tag_[11] = false; - tag_[12] = false; return 0; } @@ -282,7 +263,6 @@ namespace GridKit ScalarT v6 = y[9]; ScalarT v7 = y[10]; ScalarT vss = y[11]; - ScalarT vs = y[12]; ScalarT x1_dot = yp[0]; ScalarT x2_dot = yp[1]; @@ -292,31 +272,22 @@ namespace GridKit ScalarT x6_dot = yp[5]; ScalarT x7_dot = yp[6]; - ScalarT u = ws[0]; - ScalarT vct = ws[1]; + ScalarT u = ws[0]; 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[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] = -x5_dot + use_T2_block_ * (v4 - x5) * safe_inv_T2_; - f[5] = -x6_dot + use_T4_block_ * (v5 - x6) * safe_inv_T4_; - f[6] = -x7_dot + use_T6_block_ * (v6 - x7) * safe_inv_T6_; - + 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[8] = -v5 + bypass_T2_block_ * v4 + use_T2_block_ * (x5 + T1_ * (v4 - x5) * safe_inv_T2_); - f[9] = -v6 + bypass_T4_block_ * v5 + use_T4_block_ * (x6 + T3_ * (v5 - x6) * safe_inv_T4_); - f[10] = -v7 + bypass_T6_block_ * Ks_ * v6 + use_T6_block_ * Ks_ * T5_ * (v6 - x7) * safe_inv_T6_; - - f[11] = -vss + v7 * Math::indicator_zero(Lsmin_, Lsmax_, v7) - + Lsmin_ * Math::sigmoid(Lsmin_ - v7) - + Lsmax_ * Math::sigmoid(v7 - Lsmax_); - // Cutout: Vcl=Vcu=0 means disabled (passthrough). Otherwise, sigmoid window. - f[12] = -vs + bypass_cutout_ * vss - + use_cutout_ * vss * Math::sigmoid(vct - Vcl_) * Math::sigmoid(Vcu_ - vct); + 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); + f[11] = -vss + Math::clamp(v7, Lsmin_, Lsmax_); return 0; } @@ -330,13 +301,6 @@ namespace GridKit ws_indices_[0] = signals_.template readExternalVariableIndex(); } - ws_[1] = (Vcl_ + Vcu_) / TWO; - if (signals_.template isAttached()) - { - ws_[1] = signals_.template readExternalVariable(); - ws_indices_[1] = signals_.template readExternalVariableIndex(); - } - evaluateInternalResidual(y_.data(), yp_.data(), wb_.data(), ws_.data(), f_.data()); return 0; @@ -352,8 +316,8 @@ namespace GridKit void Ieeest::initializeMonitor() { using Variable = typename model_data_type::MonitorableVariables; - monitor_->set(Variable::vs, [this] - { return y_[12]; }); + monitor_->set(Variable::vss, [this] + { return y_[11]; }); } } // namespace Stabilizer diff --git a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md index 7a5e7b57c..d01c45be3 100644 --- a/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md +++ b/GridKit/Model/PhasorDynamics/Stabilizer/IEEEST/README.md @@ -1,10 +1,7 @@ -# IEEEST +# **IEEE Stabilizer Model (IEEEST)** -The **IEEEST** model is a standard IEEE power system stabilizer used in transient stability simulations. -It consists of a 4th-order notch filter, two lead–lag blocks, a washout block, and an output limiter with input cutout logic. - -Notes: -- The **cutout logic uses** $V_{ct}$ (as labeled in the block diagram), not $u_d$. +Standard IEEE power system stabilizer: 4th-order notch filter, two lead–lag +blocks, washout, and output limiter. ## Block Diagram @@ -14,31 +11,32 @@ Notes: Figure 1: Stabilizer IEEEST model. Figure courtesy of [PowerWorld](https://www.powerworld.com/WebHelp/) - ## Model Parameters -Symbol | Units | Description | Typical Value | Note ------- | ----- | ----------- | -------------- | ---- -$A_1$ | [s] | Notch filter denominator coefficient | 1.013 -$A_2$ | [s²] | Notch filter denominator coefficient | 0.013 -$A_3$ | [s] | Notch filter denominator coefficient | 0.0 -$A_4$ | [s²] | Notch filter denominator coefficient | 0.0 -$A_5$ | [s] | Notch filter numerator coefficient | 1.013 -$A_6$ | [s²] | Notch filter numerator coefficient | 0.113 -$T_1$ | [s] | Lead–lag 1 numerator time constant | 0.0 -$T_2$ | [s] | Lead–lag 1 denominator time constant | 0.02 -$T_3$ | [s] | Lead–lag 2 numerator time constant | 0.0 -$T_4$ | [s] | Lead–lag 2 denominator time constant | 0.0 -$T_5$ | [s] | Washout numerator time constant | 1.65 -$T_6$ | [s] | Washout denominator time constant | 1.65 -$K_s$ | [p.u.] | Stabilizer gain | 3.0 -$L_{s\min}$ | [p.u.] | Minimum stabilizer output limit | 0.1 -$L_{s\max}$ | [p.u.] | Maximum stabilizer output limit | -0.1 -$V_{cl}$ | [p.u.] | Lower input cutout threshold | 0.0 -$V_{cu}$ | [p.u.] | Upper input cutout threshold | 0.0 -$T_{delay}$ | [s] | Input time delay | 0.0 - -### Model Derived Parameters +Symbol | Units | Description | Typical Value +------------|--------|--------------------------------------|-------------- +$A_1$ | [s] | Notch denominator coefficient | 1.013 +$A_2$ | [s²] | Notch denominator coefficient | 0.013 +$A_3$ | [s] | Notch denominator coefficient | 0.0 +$A_4$ | [s²] | Notch denominator coefficient | 0.0 +$A_5$ | [s] | Notch numerator coefficient | 1.013 +$A_6$ | [s²] | Notch numerator coefficient | 0.113 +$T_1$ | [s] | Lead–lag 1 numerator time constant | 0.0 +$T_2$ | [s] | Lead–lag 1 denominator time constant | 0.02 +$T_3$ | [s] | Lead–lag 2 numerator time constant | 0.0 +$T_4$ | [s] | Lead–lag 2 denominator time constant | 0.0 +$T_5$ | [s] | Washout numerator time constant | 1.65 +$T_6$ | [s] | Washout denominator time constant | 1.65 +$K_s$ | [p.u.] | Stabilizer gain | 3.0 +$L_s^{\min}$ | [p.u.] | Minimum stabilizer output limit | -0.1 +$L_s^{\max}$ | [p.u.] | Maximum stabilizer output limit | 0.1 + +The IEEE 421.5 IEEEST also defines a cutout window ($V_{cl}$, $V_{cu}$) and an +input delay ($T_{delay}$). These parameters are accepted for input-format +compatibility but are not modeled here. + +### Derived Parameters + ```math \begin{aligned} a_0 &= 1 \\ @@ -55,71 +53,64 @@ a_4 &= A_2 A_4 #### Differential -Symbol | Units | Description | Note ------- | ----- | ----------- | ---- -$x_1$ | [-] | Notch filter state | -$x_2$ | [-] | Notch filter state | -$x_3$ | [-] | Notch filter state | -$x_4$ | [-] | Notch filter state | -$x_5$ | [-] | Lead–lag 1 state | -$x_6$ | [-] | Lead–lag 2 state | -$x_7$ | [-] | Washout state | +Symbol | Units | Description +----------------------|--------|------------ +$x_1, x_2, x_3, x_4$ | [-] | Notch filter states +$x_5$ | [-] | Lead–lag 1 state +$x_6$ | [-] | Lead–lag 2 state +$x_7$ | [-] | Washout state #### Algebraic -Symbol | Units | Description | Note ------- | ----- | ----------- | ---- -$v_4$ | [p.u.] | Notch filter output | -$v_5$ | [p.u.] | Lead–lag 1 output | -$v_6$ | [p.u.] | Lead–lag 2 output | -$v_7$ | [p.u.] | Unlimited stabilizer signal | -$V_{ss}$ | [p.u.] | Limited stabilizer signal | -$V_s$ | [p.u.] | Stabilizer output | +Symbol | Units | Description +-----------|--------|------------ +$v_4$ | [p.u.] | Notch filter output +$v_5$ | [p.u.] | Lead–lag 1 output +$v_6$ | [p.u.] | Lead–lag 2 output +$v_7$ | [p.u.] | Unlimited stabilizer signal +$V_{ss}$ | [p.u.] | Limited stabilizer signal (model output) ### External Variables -#### Differential - -None. - #### Algebraic -Symbol | Units | Description | Note ------- | ----- | ----------- | ---- -$u$ | [p.u.] | Stabilizer input signal | -$V_{ct}$ | [p.u.] | Cutout signal (compared to $V_{cl},V_{cu}$) | from the block diagram +Symbol | Units | Description +-------|--------|------------ +$u$ | [p.u.] | Stabilizer input signal ## Model Equations ### Differential Equations + ```math \begin{aligned} -\dot{x}_1 &= x_2 \\ -\dot{x}_2 &= x_3 \\ -\dot{x}_3 &= x_4 \\ -\dot{x}_4 &= -\dfrac{a_0}{a_4}x_1 - -\dfrac{a_1}{a_4}x_2 - -\dfrac{a_2}{a_4}x_3 - -\dfrac{a_3}{a_4}x_4 - +\dfrac{1}{a_4}u\\ -\dot{x}_5 &= \dfrac{1}{T_2}(v_4 - x_5) \\ -\dot{x}_6 &= \dfrac{1}{T_4}(v_5 - x_6)\\ -\dot{x}_7 &= \dfrac{1}{T_6}(v_6 - x_7) +0 &= -\dot{x}_1 + x_2 \\ +0 &= -\dot{x}_2 + x_3 \\ +0 &= -\dot{x}_3 + x_4 \\ +0 &= -\dot{x}_4 - \dfrac{a_0}{a_4}x_1 - \dfrac{a_1}{a_4}x_2 - \dfrac{a_2}{a_4}x_3 - \dfrac{a_3}{a_4}x_4 + \dfrac{1}{a_4}u \\ +0 &= -T_2 \dot{x}_5 - x_5 + v_4 \\ +0 &= -T_4 \dot{x}_6 - x_6 + v_5 \\ +0 &= -T_6 \dot{x}_7 - x_7 + v_6 \end{aligned} ``` ### Algebraic Equations + ```math \begin{aligned} 0 &= -v_4 + x_1 + A_5 x_2 + A_6 x_3 \\ -0 &= -v_5 + x_5 + \dfrac{T_1}{T_2}(v_4 - x_5) \\ -0 &= -v_6 + x_6 + \dfrac{T_3}{T_4}(v_5 - x_6) \\ -0 &= -v_7 + K_s \dfrac{T_5}{T_6}(v_6 - x_7) \\ -0 &= -V_{ss} + \min\!\big(\max(v_7, L_{s\min}), L_{s\max}\big) \\ -0 &= -V_s + -\begin{cases} -V_{ss}, & V_{cl} < V_{ct} < V_{cu} \\ -0, & \text{otherwise} -\end{cases} +0 &= -T_2(v_5 - x_5) + T_1(v_4 - x_5) \\ +0 &= -T_4(v_6 - x_6) + T_3(v_5 - x_6) \\ +0 &= -T_6 v_7 + K_s T_5(v_6 - x_7) \\ +0 &= -V_{ss} + \text{clamp}(v_7, L_s^{\min}, L_s^{\max}) \end{aligned} ``` + +The output limiter uses GridKit's smooth +[Clamp](../../../../CommonMath.md#derived-functions). + +## Initialization + +All states and their derivatives initialize to zero. The stabilizer comes +online at rest and produces signal only in response to deviations in the input +$u$. diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index 53f2d990e..43cb02cf6 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -319,16 +319,10 @@ namespace GridKit stabilizer->getSignals().template attachSignalNode(getSignal(input)); } - if (stabdata.ports.contains(IeeestPorts::cutout)) - { - IdxT cutout = stabdata.ports.at(IeeestPorts::cutout); - stabilizer->getSignals().template attachSignalNode(getSignal(cutout)); - } - if (stabdata.ports.contains(IeeestPorts::output)) { IdxT output = stabdata.ports.at(IeeestPorts::output); - stabilizer->getSignals().template assignSignalNode(getSignal(output)); + stabilizer->getSignals().template assignSignalNode(getSignal(output)); } addComponent(stabilizer); diff --git a/docs/Figures/PhasorDynamics_REECA_Diagram.png b/docs/Figures/PhasorDynamics_REECA_Diagram.png new file mode 100644 index 000000000..509adaa72 Binary files /dev/null and b/docs/Figures/PhasorDynamics_REECA_Diagram.png differ diff --git a/docs/Figures/PhasorDynamics_REGCA_Diagram.png b/docs/Figures/PhasorDynamics_REGCA_Diagram.png new file mode 100644 index 000000000..35cb516dd Binary files /dev/null and b/docs/Figures/PhasorDynamics_REGCA_Diagram.png differ diff --git a/docs/Figures/PhasorDynamics_REGCB_Diagram.png b/docs/Figures/PhasorDynamics_REGCB_Diagram.png new file mode 100644 index 000000000..6ed6747d3 Binary files /dev/null and b/docs/Figures/PhasorDynamics_REGCB_Diagram.png differ diff --git a/examples/PhasorDynamics/Large/WECC/README.md b/examples/PhasorDynamics/Large/WECC/README.md index c84cfdedb..373ed3302 100644 --- a/examples/PhasorDynamics/Large/WECC/README.md +++ b/examples/PhasorDynamics/Large/WECC/README.md @@ -44,7 +44,7 @@ The following event types are provided for this case. Only one exciter model is outstanding: - REECB1 - GAST_PTI -- REGC_A +- REGCA - REPCA1 The following examples needs to be constructed with this case. diff --git a/examples/PhasorDynamics/Medium/Hawaii/README.md b/examples/PhasorDynamics/Medium/Hawaii/README.md index 7ffca21bf..9923641ca 100644 --- a/examples/PhasorDynamics/Medium/Hawaii/README.md +++ b/examples/PhasorDynamics/Medium/Hawaii/README.md @@ -42,7 +42,7 @@ The following event types are provided for this case. The following models are not implemented in GridKit and are represented using surrogate models. - IEEEST (In GridKit, not yet added to this case) -- REGC_A +- REGCA - IEEEG1 - GGOV1 - ESST4B diff --git a/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp b/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp index 676d099df..41f6d2ca1 100644 --- a/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp +++ b/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp @@ -1,6 +1,12 @@ #pragma once +#include +#include +#include + +#include #include +#include #include namespace GridKit @@ -10,32 +16,310 @@ namespace GridKit template class SmoothnessIndicatorTests { + private: + using RealT = typename GridKit::ScalarTraits::RealT; + + static constexpr RealT kSmoothTolerance = 1.0e-2; + static constexpr RealT kNearOne = 1.0 - kSmoothTolerance; + static constexpr RealT kRoundoffTolerance = 100.0 * std::numeric_limits::epsilon(); + + static ScalarT scalar(const RealT value) + { + return static_cast(value); + } + + static bool within(const ScalarT value, + const ScalarT expected, + const ScalarT tolerance) + { + return std::abs(value - expected) < tolerance; + } + public: SmoothnessIndicatorTests() = default; ~SmoothnessIndicatorTests() = default; + TestOutcome clamp() + { + TestStatus success = true; + + const ScalarT lower = scalar(-0.25); + const ScalarT upper = scalar(0.75); + const ScalarT tol = scalar(kSmoothTolerance); + + success *= (Math::clamp(scalar(-1.0), lower, upper) < lower + tol); + success *= (Math::clamp(scalar(0.4), lower, upper) > lower); + success *= (Math::clamp(scalar(0.4), lower, upper) < upper); + success *= (Math::clamp(scalar(1.5), lower, upper) > upper - tol); + + return success.report(__func__); + } + + TestOutcome deadband() + { + TestStatus success = true; + + const ScalarT lower = scalar(-0.05); + const ScalarT upper = scalar(0.10); + const ScalarT tol = scalar(kSmoothTolerance); + + const ScalarT below_input = scalar(-1.0); + const ScalarT inside_input = scalar(0.02); + const ScalarT above_input = scalar(1.0); + const ScalarT expected_below = below_input - lower; + const ScalarT expected_above = above_input - upper; + + success *= within(Math::deadband(below_input, lower, upper), expected_below, tol); + success *= (std::abs(Math::deadband(inside_input, lower, upper)) < tol * tol); + success *= within(Math::deadband(above_input, lower, upper), expected_above, tol); + + const ScalarT lower_breakpoint = Math::deadband(lower, lower, upper); + const ScalarT upper_breakpoint = Math::deadband(upper, lower, upper); + + success *= (lower_breakpoint < scalar(0.0)); + success *= (upper_breakpoint > scalar(0.0)); + success *= (std::abs(lower_breakpoint) < tol); + success *= (std::abs(upper_breakpoint) < tol); + + const ScalarT x = scalar(-0.4); + success *= (std::abs(Math::deadband(x, lower, upper) + - (x - Math::clamp(x, lower, upper))) + < scalar(kRoundoffTolerance)); + + const ScalarT far_above_input = scalar(4.0); + const ScalarT far_below_input = scalar(-4.0); + const ScalarT expected_far_above = far_above_input - upper; + const ScalarT expected_far_below = far_below_input - lower; + + success *= std::isfinite(Math::deadband(far_above_input, lower, upper)); + success *= within(Math::deadband(far_above_input, lower, upper), expected_far_above, tol); + success *= std::isfinite(Math::deadband(far_below_input, lower, upper)); + success *= within(Math::deadband(far_below_input, lower, upper), expected_far_below, tol); + + const ScalarT point = scalar(0.25); + const ScalarT above_point = scalar(0.75); + const ScalarT below_point = scalar(-0.25); + success *= (std::abs(Math::deadband(above_point, point, point) - (above_point - point)) + < scalar(kRoundoffTolerance)); + success *= (std::abs(Math::deadband(below_point, point, point) - (below_point - point)) + < scalar(kRoundoffTolerance)); + + return success.report(__func__); + } + + TestOutcome limitIndicators() + { + TestStatus success = true; + + const ScalarT limit_min = scalar(0.0); + const ScalarT limit_max = scalar(3.0); + const ScalarT near_zero = scalar(kSmoothTolerance); + const ScalarT near_one = scalar(kNearOne); + + success *= (Math::above(scalar(1.0), limit_min) > near_one); + success *= (Math::above(scalar(-0.2), limit_min) < near_zero); + success *= (Math::below(scalar(1.0), limit_max) > near_one); + success *= (Math::below(scalar(3.2), limit_max) < near_zero); + + success *= (Math::inside(scalar(1.5), limit_min, limit_max) > near_one); + success *= (Math::inside(scalar(-0.2), limit_min, limit_max) < near_zero); + success *= (Math::inside(scalar(3.2), limit_min, limit_max) < near_zero); + + success *= (Math::outside(scalar(1.5), limit_min, limit_max) < near_zero); + success *= (Math::outside(scalar(-0.2), limit_min, limit_max) > near_one); + success *= (Math::outside(scalar(3.2), limit_min, limit_max) > near_one); + + const ScalarT x = scalar(1.5); + success *= (std::abs(Math::inside(x, limit_min, limit_max) + + Math::outside(x, limit_min, limit_max) + - scalar(1.0)) + < scalar(kRoundoffTolerance)); + + return success.report(__func__); + } + + TestOutcome slew() + { + TestStatus success = true; + + const ScalarT rate = scalar(0.5); + const ScalarT tol = scalar(kSmoothTolerance); + + success *= within(Math::slew(scalar(2.0), rate), rate, tol); + success *= within(Math::slew(scalar(-2.0), rate), -rate, tol); + success *= within(Math::slew(scalar(0.2), rate), scalar(0.2), tol); + success *= within(Math::slew(scalar(-0.2), rate), scalar(-0.2), tol); + + return success.report(__func__); + } + + TestOutcome linseg() + { + TestStatus success = true; + + const ScalarT lower = scalar(0.0); + const ScalarT upper = scalar(2.0); + const ScalarT height = scalar(4.0); + const ScalarT tol = scalar(kSmoothTolerance); + + const ScalarT below_input = scalar(-1.0); + const ScalarT midpoint_input = scalar(1.0); + const ScalarT above_input = scalar(3.0); + const ScalarT expected_mid = height * (midpoint_input - lower) / (upper - lower); + const ScalarT negative_height = -height; + const ScalarT expected_mid_neg = negative_height * (midpoint_input - lower) / (upper - lower); + + success *= (Math::linseg(below_input, lower, upper, height) < tol); + success *= within(Math::linseg(midpoint_input, lower, upper, height), expected_mid, tol); + success *= within(Math::linseg(above_input, lower, upper, height), height, tol); + success *= within(Math::linseg(midpoint_input, lower, upper, negative_height), expected_mid_neg, tol); + + return success.report(__func__); + } + + TestOutcome ramp() + { + TestStatus success = true; + + const ScalarT tol = scalar(kSmoothTolerance); + const ScalarT roundoff = scalar(kRoundoffTolerance); + const ScalarT tau = scalar(1.0 / 240.0); + const ScalarT at_zero = tau * std::log(scalar(2.0)); + const ScalarT far_above = scalar(4.0); + const ScalarT far_below = scalar(-4.0); + const ScalarT lower = scalar(-0.25); + const ScalarT upper = scalar(0.75); + const ScalarT x = scalar(0.4); + const ScalarT smooth_clip = lower + Math::ramp(x - lower) + - Math::ramp(x - upper); + + success *= (Math::ramp(scalar(1.0)) > scalar(kNearOne)); + success *= (Math::ramp(scalar(-1.0)) < tol); + success *= (std::abs(Math::ramp(scalar(0.0)) - at_zero) < roundoff); + success *= (Math::ramp(-tol) > scalar(0.0)); + success *= (Math::ramp(tol) > Math::ramp(scalar(0.0))); + success *= std::isfinite(Math::ramp(far_above)); + success *= (Math::ramp(far_above) > far_above - tol); + success *= std::isfinite(Math::ramp(far_below)); + success *= (Math::ramp(far_below) < roundoff); + + success *= (smooth_clip > lower); + success *= (smooth_clip < upper); + success *= std::isfinite(Math::clamp(far_above, lower, upper)); + success *= (Math::clamp(far_above, lower, upper) < upper + roundoff); + success *= std::isfinite(Math::clamp(far_below, lower, upper)); + success *= (Math::clamp(far_below, lower, upper) > lower - roundoff); + + return success.report(__func__); + } + + TestOutcome minMax() + { + TestStatus success = true; + + const ScalarT high = scalar(2.0); + const ScalarT low = scalar(-1.0); + const ScalarT tol = scalar(kSmoothTolerance); + + success *= within(Math::max(high, low), high, tol); + success *= within(Math::max(low, high), high, tol); + + success *= within(Math::min(high, low), low, tol); + success *= within(Math::min(low, high), low, tol); + + using Variable = GridKit::DependencyTracking::Variable; + + // This mirrors model equations where a differentiable state is + // limited by a plain real parameter. + const RealT parameter_bound = kSmoothTolerance; + const Variable state_below_bound{-1.0, 0}; + const Variable state_above_bound{1.0, 1}; + const auto lower_bounded_state = Math::max(state_below_bound, parameter_bound); + const auto upper_bounded_state = Math::min(state_above_bound, parameter_bound); + + static_assert(std::is_same::type, Variable>::value, + "Mixed-type max should keep the differentiable scalar type."); + static_assert(std::is_same::type, Variable>::value, + "Mixed-type min should keep the differentiable scalar type."); + + success *= (std::abs(lower_bounded_state.getValue() - parameter_bound) < kSmoothTolerance); + success *= (std::abs(upper_bounded_state.getValue() - parameter_bound) < kSmoothTolerance); + + const ScalarT x = scalar(0.4); + const ScalarT y = scalar(-0.7); + success *= (std::abs(Math::min(x, y) + Math::max(x, y) - (x + y)) + < scalar(kRoundoffTolerance)); + + const ScalarT point = scalar(0.25); + const ScalarT bias = std::log(scalar(2.0)) / scalar(240.0); + + success *= (std::abs(Math::max(point, point) - (point + bias)) + < scalar(kRoundoffTolerance)); + success *= (std::abs(Math::min(point, point) - (point - bias)) + < scalar(kRoundoffTolerance)); + + return success.report(__func__); + } + TestOutcome antiWindupIndicator() { TestStatus success = true; - const ScalarT limit_min = 0.0; - const ScalarT limit_max = 3.0; + const ScalarT limit_min = scalar(0.0); + const ScalarT limit_max = scalar(3.0); + const ScalarT inside = scalar(1.5); + const ScalarT above_limit = scalar(3.2); + const ScalarT below_limit = scalar(-0.2); + const ScalarT f_positive = scalar(0.03); + const ScalarT f_negative = -f_positive; + const ScalarT near_zero = scalar(kSmoothTolerance); + const ScalarT near_one = scalar(kNearOne); + + // Inside the limits, both increasing and decreasing dynamics pass. + success *= (Math::indicator(inside, f_positive, limit_min, limit_max) > near_one); + success *= (Math::indicator(inside, f_negative, limit_min, limit_max) > near_one); + + // Above the upper limit, increasing dynamics are blocked. + success *= (Math::indicator(above_limit, f_positive, limit_min, limit_max) < near_zero); + + // Above the upper limit, decreasing dynamics restore the state. + success *= (Math::indicator(above_limit, f_negative, limit_min, limit_max) > near_one); + + // Below the lower limit, decreasing dynamics are blocked. + success *= (Math::indicator(below_limit, f_negative, limit_min, limit_max) < near_zero); + + // Below the lower limit, increasing dynamics restore the state. + success *= (Math::indicator(below_limit, f_positive, limit_min, limit_max) > near_one); - // Inside the limits the indicator passes dynamics through, regardless - // of the sign of f: value is close to 1. - success *= (Math::indicator(limit_min, limit_max, static_cast(1.5), static_cast(0.01)) > static_cast(0.99)); + return success.report(__func__); + } + + TestOutcome antiWindup() + { + TestStatus success = true; - // Above the upper limit with f pushing further out: blocked (≈ 0). - success *= (Math::indicator(limit_min, limit_max, static_cast(3.2), static_cast(0.01)) < static_cast(0.1)); + const ScalarT limit_min = scalar(0.0); + const ScalarT limit_max = scalar(3.0); + const ScalarT inside = scalar(1.5); + const ScalarT above_limit = scalar(3.2); + const ScalarT below_limit = scalar(-0.2); - // Above the upper limit but f pulling back in: passed (≈ 1). - success *= (Math::indicator(limit_min, limit_max, static_cast(3.2), static_cast(-0.01)) > static_cast(0.9)); + const ScalarT f_positive = scalar(0.03); + const ScalarT f_negative = -f_positive; + const ScalarT blocked_tolerance = f_positive * scalar(kSmoothTolerance); - // Below the lower limit with f pushing further out: blocked (≈ 0). - success *= (Math::indicator(limit_min, limit_max, static_cast(-0.2), static_cast(-0.01)) < static_cast(0.1)); + success *= (std::abs(Math::antiwindup(inside, f_positive, limit_min, limit_max) + - Math::indicator(inside, f_positive, limit_min, limit_max) * f_positive) + < scalar(kRoundoffTolerance)); - // Below the lower limit but f pulling back in: passed (≈ 1). - success *= (Math::indicator(limit_min, limit_max, static_cast(-0.2), static_cast(0.01)) > static_cast(0.9)); + success *= (std::abs(Math::antiwindup(above_limit, f_positive, limit_min, limit_max)) < blocked_tolerance); + success *= within(Math::antiwindup(above_limit, f_negative, limit_min, limit_max), + f_negative, + blocked_tolerance); + success *= (std::abs(Math::antiwindup(below_limit, f_negative, limit_min, limit_max)) < blocked_tolerance); + success *= within(Math::antiwindup(below_limit, f_positive, limit_min, limit_max), + f_positive, + blocked_tolerance); return success.report(__func__); } diff --git a/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp b/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp index 3da97577f..233c36ddd 100644 --- a/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp +++ b/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp @@ -6,7 +6,15 @@ int main() GridKit::Testing::SmoothnessIndicatorTests test; + result += test.clamp(); + result += test.deadband(); + result += test.limitIndicators(); + result += test.slew(); + result += test.linseg(); + result += test.ramp(); + result += test.minMax(); result += test.antiWindupIndicator(); + result += test.antiWindup(); return result.summary(); } diff --git a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp index d35d774a6..0c1c77d66 100644 --- a/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp +++ b/tests/UnitTests/PhasorDynamics/StabilizerIeeestTests.hpp @@ -49,24 +49,24 @@ namespace GridKit { TestStatus success = true; - // Create signal nodes for input (u) and output (Vs) + // Create signal nodes for input (u) and output (Vss) PhasorDynamics::SignalNode u_node; - PhasorDynamics::SignalNode vs_node; + PhasorDynamics::SignalNode vss_node; ScalarT u_value{0.0}; - IdxT u_index = 13; // beyond internal variables - ScalarT vs_value{0.0}; - IdxT vs_index = INVALID_INDEX; + IdxT u_index = 12; // beyond internal variables + ScalarT vss_value{0.0}; + IdxT vss_index = INVALID_INDEX; // Link signal nodes to backing storage u_node.set(&u_value, &u_index); - vs_node.set(&vs_value, &vs_index); + vss_node.set(&vss_value, &vss_index); auto data = makeTestData(); PhasorDynamics::Stabilizer::Ieeest stab(data); - // Wire: stabilizer reads u_node as input, writes vs_node as output + // Wire: stabilizer reads u_node as input, writes vss_node as output stab.getSignals().template attachSignalNode(&u_node); - stab.getSignals().template assignSignalNode(&vs_node); + stab.getSignals().template assignSignalNode(&vss_node); stab.allocate(); success *= (stab.verify() == 0); @@ -85,9 +85,9 @@ namespace GridKit } // Verify output signal is linked and reads the correct value - success *= vs_node.linked(); - success *= (vs_node.getVariableIndex() == 12); - success *= isEqual(vs_node.read(), static_cast(0.0), tol); + success *= vss_node.linked(); + success *= (vss_node.getVariableIndex() == 11); + success *= isEqual(vss_node.read(), static_cast(0.0), tol); return success.report(__func__); } @@ -103,20 +103,20 @@ namespace GridKit TestStatus success = true; PhasorDynamics::SignalNode u_node; - PhasorDynamics::SignalNode vs_node; + PhasorDynamics::SignalNode vss_node; ScalarT u_value{0.5}; - IdxT u_index = 13; - ScalarT vs_value{0.0}; - IdxT vs_index = INVALID_INDEX; + IdxT u_index = 12; + ScalarT vss_value{0.0}; + IdxT vss_index = INVALID_INDEX; u_node.set(&u_value, &u_index); - vs_node.set(&vs_value, &vs_index); + vss_node.set(&vss_value, &vss_index); auto data = makeTestData(); PhasorDynamics::Stabilizer::Ieeest stab(data); stab.getSignals().template attachSignalNode(&u_node); - stab.getSignals().template assignSignalNode(&vs_node); + stab.getSignals().template assignSignalNode(&vss_node); stab.allocate(); stab.initialize(); @@ -129,24 +129,23 @@ namespace GridKit 0.28, // f[1]: -x2_dot + x3 0.37, // f[2]: -x3_dot + x4 1.0975, // f[3]: -x4_dot + (-a0*x1 - a1*x2 - a2*x3 - a3*x4 + u) / a4 - 0.25, // f[4]: -x5_dot + (v4 - x5) / T2 - 0.24, // f[5]: -x6_dot + (v5 - x6) / T4 - -0.01, // f[6]: -x7_dot + (v6 - x7) / T6 + 0.25, // f[4]: -T2*x5_dot - x5 + v4 + 0.24, // f[5]: -T4*x6_dot - x6 + v5 + -0.05, // f[6]: -T6*x7_dot - x7 + v6 -0.42, // f[7]: -v4 + x1 + A5*x2 + A6*x3 - -0.25, // f[8]: -v5 + x5 + (T1/T2)*(v4 - x5) - -0.31, // f[9]: -v6 + x6 + (T3/T4)*(v5 - x6) - 1.15, // f[10]: -v7 + Ks*(T5/T6)*(v6 - x7) + -0.25, // f[8]: -T2*(v5 - x5) + T1*(v4 - x5) + -0.31, // f[9]: -T4*(v6 - x6) + T3*(v5 - x6) + 5.75, // f[10]: -T6*v7 + Ks*T5*(v6 - x7) 0.0, // f[11]: limiter (v7=0.05 within [-0.1, 0.1]) - 0.0, // f[12]: cutout (Vct=0.75 within [0.0, 1.5]) }; - // Looser tolerance for f[11],f[12] — smooth sigmoid approximations + // Looser tolerance for f[11] — Math::clamp is a smooth ramp approximation. const auto loose_tol = static_cast(1.0e-4); auto& residual = stab.getResidual(); for (size_t i = 0; i < res_answer.size(); ++i) { - auto test_tol = (i >= 11) ? loose_tol : static_cast(10 * std::numeric_limits::epsilon()); + auto test_tol = (i == 11) ? loose_tol : static_cast(10 * std::numeric_limits::epsilon()); if (!isEqual(residual[i], res_answer[i], test_tol)) { std::cout << "Incorrect result for residual " << i << ": " @@ -157,7 +156,7 @@ namespace GridKit } // Verify output signal reads the stabilizer output - success *= isEqual(vs_node.read(), static_cast(0.05), loose_tol); + success *= isEqual(vss_node.read(), static_cast(0.05), loose_tol); return success.report(__func__); } @@ -196,18 +195,18 @@ namespace GridKit // Set up signal nodes with DependencyTracking scalar type PhasorDynamics::SignalNode u_node; - PhasorDynamics::SignalNode vs_node; + PhasorDynamics::SignalNode vss_node; DepVar u_value{0.5}; - IdxT u_index = 13; - DepVar vs_value{0.0}; - IdxT vs_index = INVALID_INDEX; + IdxT u_index = 12; + DepVar vss_value{0.0}; + IdxT vss_index = INVALID_INDEX; u_node.set(&u_value, &u_index); - vs_node.set(&vs_value, &vs_index); + vss_node.set(&vss_value, &vss_index); PhasorDynamics::Stabilizer::Ieeest stab(ieeestdata); stab.getSignals().template attachSignalNode(&u_node); - stab.getSignals().template assignSignalNode(&vs_node); + stab.getSignals().template assignSignalNode(&vss_node); stab.allocate(); stab.initialize(); @@ -287,18 +286,18 @@ namespace GridKit PhasorDynamics::Stabilizer::IeeestData ieeestdata) { PhasorDynamics::SignalNode u_node; - PhasorDynamics::SignalNode vs_node; + PhasorDynamics::SignalNode vss_node; ScalarT u_value{0.5}; - IdxT u_index = 13; - ScalarT vs_value{0.0}; - IdxT vs_index = INVALID_INDEX; + IdxT u_index = 12; + ScalarT vss_value{0.0}; + IdxT vss_index = INVALID_INDEX; u_node.set(&u_value, &u_index); - vs_node.set(&vs_value, &vs_index); + vss_node.set(&vss_value, &vss_index); PhasorDynamics::Stabilizer::Ieeest stab(ieeestdata); stab.getSignals().template attachSignalNode(&u_node); - stab.getSignals().template assignSignalNode(&vs_node); + stab.getSignals().template assignSignalNode(&vss_node); stab.allocate(); stab.initialize(); @@ -331,7 +330,7 @@ namespace GridKit PhasorDynamics::Stabilizer::IeeestData data; data.device_class = "stabilizer"; data.disambiguation_string = "ieeest_test"; - data.monitored_variables.insert(PhasorDynamics::Stabilizer::IeeestMonitorableVariables::vs); + data.monitored_variables.insert(PhasorDynamics::Stabilizer::IeeestMonitorableVariables::vss); data.parameters[Params::A1] = 0.1; data.parameters[Params::A2] = 0.2; @@ -349,7 +348,7 @@ namespace GridKit data.parameters[Params::Lsmin] = -0.1; data.parameters[Params::Lsmax] = 0.1; data.parameters[Params::Vcl] = 0.0; - data.parameters[Params::Vcu] = 1.5; + data.parameters[Params::Vcu] = 0.0; data.parameters[Params::Tdelay] = 0.0; return data; @@ -372,8 +371,7 @@ namespace GridKit stab.y()[8] = 0.9; // v5 stab.y()[9] = 1.0; // v6 stab.y()[10] = 0.05; // v7 (within limiter range) - stab.y()[11] = 0.05; // Vss - stab.y()[12] = 0.05; // Vs + stab.y()[11] = 0.05; // Vss (model output) stab.yp()[0] = 0.01; // x1_dot stab.yp()[1] = 0.02; // x2_dot @@ -402,7 +400,6 @@ namespace GridKit stab.y()[9].setValue(1.0); stab.y()[10].setValue(0.05); stab.y()[11].setValue(0.05); - stab.y()[12].setValue(0.05); stab.yp()[0].setValue(0.01); stab.yp()[1].setValue(0.02);