diff --git a/GridKit/Model/PhasorDynamics/Controller/REECB/README.md b/GridKit/Model/PhasorDynamics/Controller/REECB/README.md index 0dcf468b8..14617aad2 100644 --- a/GridKit/Model/PhasorDynamics/Controller/REECB/README.md +++ b/GridKit/Model/PhasorDynamics/Controller/REECB/README.md @@ -152,7 +152,6 @@ $e_Q$ | [p.u.] | Reactive-power error $V_Q^\mathrm{PI}$ | [p.u.] | Reactive-power PI output | $e_V^\mathrm{PI}$ | [p.u.] | Voltage-control error | $r_P^\mathrm{ord}$ | [p.u./s] | Limited active-power order rate | Component base -$I_L^{\max}$ | [p.u.] | Current-circle continuation state | Component base $I_L^\mathrm{cap}$ | [p.u.] | Off-axis current capacity | Component base $I_q^{\max}$ | [p.u.] | Reactive-current limit | Component base $I_p^{\max}$ | [p.u.] | Active-current limit | Component base @@ -209,8 +208,7 @@ $P^\mathrm{ref}$ | [p.u.] | Unknown | External active-power reference 0 &= -V_Q^\mathrm{PI}+\text{clamp}(K_\mathrm{qp}e_Q+x_Q^\mathrm{PI};\,V^{\min},V^{\max}) \\ 0 &= -e_V^\mathrm{PI}+s_Q^\mathrm{PI}V_Q^\mathrm{PI}+s_V^\mathrm{ref}Q^\mathrm{ext}-s_QV^\mathrm{meas} \\ 0 &= -r_P^\mathrm{ord}+\text{aslew}\!\left(\dfrac{k_\mathrm{base}P^\mathrm{ref}-P^\mathrm{ord}}{T_\mathrm{pord}};\,R_P^{\min},R_P^{\max}\right) \\ - 0 &= -I_L^{\max}\sqrt{(I_L^{\max})^2+\epsilon_0}+(I^{\max}-I^\mathrm{high})(I^{\max}+I^\mathrm{high}) \\ - 0 &= -I_L^\mathrm{cap}+\dfrac{(I_L^{\max})^2}{\sqrt{(I_L^{\max})^2+\epsilon_0}} \\ + 0 &= -I_L^\mathrm{cap}+\text{sqrtramp}\big((I^{\max}-I^\mathrm{high})(I^{\max}+I^\mathrm{high})\big) \\ 0 &= -I_q^{\max}+s_\mathrm{pq}I_L^\mathrm{cap}+s_\mathrm{pq}^\mathrm{off}I^{\max} \\ 0 &= -I_p^{\max}+s_\mathrm{pq}I^{\max}+s_\mathrm{pq}^\mathrm{off}I_L^\mathrm{cap} \\ 0 &= -I_q^\mathrm{base}+\text{clamp}(K_\mathrm{vp}e_V^\mathrm{PI}+x_V^\mathrm{PI};\,-I_q^{\max},I_q^{\max}) \\ @@ -224,6 +222,7 @@ Here $I^\mathrm{high}=s_\mathrm{pq}k_\mathrm{base}I_p^\mathrm{cmd} +s_\mathrm{pq}^\mathrm{off}k_\mathrm{base}I_q^\mathrm{cmd}$ and $\epsilon_0=100\epsilon_\mathrm{machine}$. +[Appendix C](#appendix-c-sqrtramp) defines `sqrtramp`. CommonMath defines the [`antiwindup`](../../../../CommonMath.md#antiwindup) and [smooth limiter](../../../../CommonMath.md#derived-functions) functions used in @@ -272,23 +271,23 @@ and $I_q=k_\mathrm{base}I_q^\mathrm{cmd}$ be the component-base commands. Let $I_q^\mathrm{need}$ be the smooth-clamp input magnitude required to reproduce $I_q$, including the Volt/VAr recovery margin. Initialization raises -$I^\max$, if needed, to the smallest finite limit satisfying +$I^\max$, if needed, to a finite limit satisfying ```math \begin{aligned} h &= s_\mathrm{pq}I_p+s_\mathrm{pq}^\mathrm{off}|I_q| \\ \ell &= s_\mathrm{pq}I_q^\mathrm{need}+s_\mathrm{pq}^\mathrm{off}I_p \\ I &\ge \text{max}(I^\max,h,\ell,I_q^\mathrm{need}) \\ - x\sqrt{x^2+\epsilon_0} &= (I-h)(I+h) \\ - \dfrac{x^2}{\sqrt{x^2+\epsilon_0}} &\ge \ell. + \text{sqrtramp}\big((I-h)(I+h)\big) &\ge \ell. \end{aligned} ``` -The solution initializes $I^\max\leftarrow I$, $I_L^\max\leftarrow x$, and -$I_L^\mathrm{cap}\leftarrow x^2/\sqrt{x^2+\epsilon_0}$. The operating point is -rejected if no finite solution exists or the reactive-current injection is -incompatible. Q, V, and P limits are expanded as needed; each adjustment logs -a warning. +The solution initializes +$I^\max\leftarrow I$ and +$I_L^\mathrm{cap}\leftarrow\text{sqrtramp}\big((I-h)(I+h)\big)$. The +operating point is rejected if no finite solution exists or the +reactive-current injection is incompatible. Q, V, and P limits are expanded as +needed; each adjustment logs a warning. ```math \begin{aligned} @@ -375,11 +374,11 @@ Output | Units | Description | Note - `initializationAndSignals()` checks initialization, signals, monitors, and power bases. - `initializationDomain()` checks rejected inputs and limit expansion. - `initializationExactness()` checks endpoint and current-circle initialization. -- `residualEquations()` checks the fixed residual answer key. +- `residualEquations()` checks residual equations at a hand-computable state. - `selectorConfigurations()` checks selectors and optional ports. - `voltVarReferenceBase()` checks `qext` units. - `reactiveControl()` checks the reactive-control paths. -- `activeCurrentControl()` checks active-current control and current priority. +- `activeCurrentControl()` checks active-current control, current priority, and the current circle. - `jacobian()` compares the Enzyme and dependency-tracking Jacobians. - `regcaReecbRepca()` checks the public-signal REGCA-REECB-REPCA control loop. - `reecb()` checks construction through the production system-data path. @@ -412,6 +411,26 @@ where $\rho$ is GridKit's smooth reduces to $\text{clamp}(f;\ell,u)$; the smooth form preserves $\text{aslew}(0;\ell,u)=0$. +## Appendix C: `sqrtramp` + +The off-axis capacity approaches $\sqrt{\max(0,s)}$, where +$s=(I^{\max}-I^\mathrm{high})(I^{\max}+I^\mathrm{high})$. Let +$\kappa=\epsilon_0/\mu$, with $\mu$ GridKit's +[smoothing scale](../../../../CommonMath.md#primitives). REECB uses + +```math +\begin{aligned} + h_\kappa(s) &= \tfrac{1}{2}\left(s+\sqrt{s^2+\kappa^2}\right), \\ + \text{sqrtramp}(s) &= + \dfrac{\mu h_\kappa(s)}{\sqrt{\mu^2 h_\kappa(s)+1}}. +\end{aligned} +``` + +Both softenings are differentiable for every finite $s$. The narrow hinge +leaves a negligible positive capacity when the circle is closed or +over-driven, while preserving the current-command dependence through the +transition. Increasing $\mu$ recovers the exact current-circle leg. + [^wecc-reecb-specification]: [WECC REMTF, *Generic Solar Photovoltaic System Dynamic Simulation Model Specification*](https://www.wecc.org/sites/default/files/documents/meeting/2024/WECC-Solar-PV-Dynamic-Model-Specification-September-2012.pdf), September 2012. [^reecb-mva-base]: The [WECC Central Station Photovoltaic Power Plant Model Validation Guideline](https://www.wecc.org/sites/default/files/documents/program/2024/Central%20Station%20Photovoltaic%20Power%20Plant%20Model%20Validation%20Guideline%20June%2017%202015.pdf) diff --git a/GridKit/Model/PhasorDynamics/Controller/REECB/Reecb.hpp b/GridKit/Model/PhasorDynamics/Controller/REECB/Reecb.hpp index fd08c3794..9188a0920 100644 --- a/GridKit/Model/PhasorDynamics/Controller/REECB/Reecb.hpp +++ b/GridKit/Model/PhasorDynamics/Controller/REECB/Reecb.hpp @@ -44,7 +44,6 @@ namespace GridKit VPIQ, ///< \f$V_Q^\mathrm{PI}\f$ Algebraic reactive-power PI output [p.u.] EPIV, ///< \f$e_V^\mathrm{PI}\f$ Algebraic voltage-control error [p.u.] RPORD, ///< \f$r_P^\mathrm{ord}\f$ Algebraic limited active-power order rate [p.u./s] - ILMAX, ///< \f$I_L^\max\f$ Algebraic current-circle continuation state on component base [p.u.] ILCAP, ///< \f$I_L^\mathrm{cap}\f$ Algebraic off-axis current capacity on component base [p.u.] IQMAX, ///< \f$I_q^\max\f$ Algebraic reactive-current limit on component base [p.u.] IPMAX, ///< \f$I_p^\max\f$ Algebraic active-current limit on component base [p.u.] @@ -106,7 +105,7 @@ namespace GridKit using InternalVariablesT = ReecbInternalVariables; using ExternalVariablesT = ReecbExternalVariables; - /// Current-circle regularization and initialization reconstruction tolerance. + /// Initialization reconstruction tolerance. static constexpr RealT INITIALIZATION_TOLERANCE = static_cast(100.0) * std::numeric_limits::epsilon(); @@ -139,12 +138,16 @@ namespace GridKit ScalarT* f); private: + /// Hinge width chosen so closed-circle leakage stays below the + /// initialization tolerance [p.u. current squared]. + static constexpr RealT CURRENT_CIRCLE_KNEE = + INITIALIZATION_TOLERANCE / Math::MU; + struct InitialPoint; struct InitialCurrentLimit { RealT total_limit; - RealT continuation; RealT off_axis_capacity; }; @@ -154,17 +157,18 @@ namespace GridKit /// Smooth anti-windup derivative within a moving symmetric band. [[gnu::always_inline]] static inline ScalarT awband(ScalarT state, ScalarT rate, ScalarT band); - /// Current-circle continuation state for an initial component-base limit. - static RealT circleState(RealT imax, RealT high); + /// Smooth nonnegative root of a current-circle squared radius. + template + [[gnu::always_inline]] static inline ValueT sqrtramp(ValueT x); - /// Off-axis component-base capacity provided by a continuation state. - static RealT capacity(RealT ilmax); + /// Analytic nonnegative seed for inverting sqrtramp(). + static RealT isqrtramp(RealT y); - /// Bisect an initial-limit bracket to its first upper-side point. - template - static RealT bisect(RealT a, RealT b, FuncT below); + /// Overflow-resistant difference of squares. + template + [[gnu::always_inline]] static inline ValueT circleSquare(RealT limit, ValueT high); - /// Solve the smallest feasible initial limit at or above `lower`. + /// Solve a feasible initial limit at or above `lower`. static std::optional solveInitialLimit(RealT lower, RealT high, RealT low); bool buildInitialPoint(InitialPoint& point); diff --git a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp index ffa117d83..6d39b88ec 100644 --- a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp @@ -12,7 +12,6 @@ #include #include #include -#include #include #include @@ -348,7 +347,6 @@ namespace GridKit const auto VPIQ = static_cast(ReecbInternalVariables::VPIQ); const auto EPIV = static_cast(ReecbInternalVariables::EPIV); const auto RPORD = static_cast(ReecbInternalVariables::RPORD); - const auto ILMAX = static_cast(ReecbInternalVariables::ILMAX); const auto ILCAP = static_cast(ReecbInternalVariables::ILCAP); const auto IQMAX = static_cast(ReecbInternalVariables::IQMAX); const auto IPMAX = static_cast(ReecbInternalVariables::IPMAX); @@ -442,7 +440,6 @@ namespace GridKit return false; } const RealT imax = current_limit->total_limit; - const RealT ilmax0 = current_limit->continuation; const RealT ilcap0 = current_limit->off_axis_capacity; RealT iqmax0 = imax; @@ -614,7 +611,7 @@ namespace GridKit const RealT iqcmd_check = Math::clamp(iqraw_check, -iqmax0, iqmax0); const RealT ipcmd_check = Math::clamp(pord0 / vmeas_safe0, ZERO, ipmax0); - if (!std::isfinite(imax) || !std::isfinite(ilmax0) || !std::isfinite(ilcap0) + if (!std::isfinite(imax) || !std::isfinite(ilcap0) || !std::isfinite(iqmax0) || !std::isfinite(ipmax0) || !std::isfinite(ipraw0) || !std::isfinite(iqraw0) || !std::isfinite(pord0) || !std::isfinite(pref0_system) || !std::isfinite(qtarget0) || !std::isfinite(qref0) @@ -654,7 +651,6 @@ namespace GridKit point.variables[VPIQ] = vpiq0; point.variables[EPIV] = epiv0; point.variables[RPORD] = ZERO; - point.variables[ILMAX] = ilmax0; point.variables[ILCAP] = ilcap0; point.variables[IQMAX] = iqmax0; point.variables[IPMAX] = ipmax0; @@ -923,7 +919,6 @@ namespace GridKit const auto VPIQ = static_cast(ReecbInternalVariables::VPIQ); const auto EPIV = static_cast(ReecbInternalVariables::EPIV); const auto RPORD = static_cast(ReecbInternalVariables::RPORD); - const auto ILMAX = static_cast(ReecbInternalVariables::ILMAX); const auto ILCAP = static_cast(ReecbInternalVariables::ILCAP); const auto IQMAX = static_cast(ReecbInternalVariables::IQMAX); const auto IPMAX = static_cast(ReecbInternalVariables::IPMAX); @@ -953,7 +948,6 @@ namespace GridKit const ScalarT vpiq = y[VPIQ]; const ScalarT epiv = y[EPIV]; const ScalarT rpord = y[RPORD]; - const ScalarT ilmax = y[ILMAX]; const ScalarT ilcap = y[ILCAP]; const ScalarT iqmax = y[IQMAX]; const ScalarT ipmax = y[IPMAX]; @@ -984,7 +978,6 @@ namespace GridKit const ScalarT q_pi_state = Kqp_ * eq + xpiq; const ScalarT v_pi_state = Kvp_ * epiv + xpiv; const ScalarT fpord = (pref - pord) / Tpord_; - const ScalarT ilnorm = std::sqrt(ilmax * ilmax + INITIALIZATION_TOLERANCE); // Select before the factored square to avoid 0 * inf on the inactive path. const ScalarT high = pq_on_ * ipcmd + pq_off_ * iqcmd; const ScalarT q_pi_rate = q_pi_on_ * sdip * Math::antiwindup(q_pi_state, Kqi_ * eq, Vmin_, Vmax_); @@ -1012,8 +1005,7 @@ namespace GridKit f[VPIQ] = -vpiq + Math::clamp(q_pi_state, Vmin_, Vmax_); f[EPIV] = -epiv + q_pi_on_ * vpiq + v_ref_on_ * extref - q_on_ * vmeas; f[RPORD] = -rpord + aslew(fpord, dPmin_, dPmax_); - f[ILMAX] = -ilmax * ilnorm + (Imax_ - high) * (Imax_ + high); - f[ILCAP] = -ilcap + (ilmax / ilnorm) * ilmax; + f[ILCAP] = -ilcap + sqrtramp(circleSquare(Imax_, high)); f[IQMAX] = -iqmax + pq_on_ * ilcap + pq_off_ * Imax_; f[IPMAX] = -ipmax + pq_on_ * Imax_ + pq_off_ * ilcap; f[IQBASE] = -iqbase + Math::clamp(v_pi_state, -iqmax, iqmax); @@ -1071,106 +1063,94 @@ namespace GridKit } /** - * @brief Compute the initial current-circle continuation state + * @brief Smooth nonnegative root of a current-circle squared radius * - * Solves the implemented `ILMAX` row for its nonnegative continuation - * state at a total-current limit and priority-axis current. + * Softens both the hinge and root so the current-circle row remains + * differentiable as the circle opens and closes. The scaled norm replaces + * `std::hypot`, which is unavailable for dependency-tracking scalars, and + * avoids squaring a large argument. The conjugate avoids subtractive + * cancellation without branching on the differentiable value. * - * @param[in] imax Total-current limit on the component base. - * @param[in] high Priority-axis current command on the component base. - * @return Initial `ILMAX` state on the component base, or a quiet NaN - * when the circle geometry is invalid or unrepresentable. - * @pre `imax` and `high` are finite and @f$0 \le high \le imax@f$. + * @tparam ValueT Differentiable scalar or plain real. + * @param[in] x Argument. + * @return Softened off-axis current on the component base. */ template - typename Reecb::RealT - Reecb::circleState(RealT imax, RealT high) + template + [[gnu::always_inline]] inline ValueT + Reecb::sqrtramp(ValueT x) { - const RealT rhs = (imax - high) * (imax + high); - if (!std::isfinite(rhs) || rhs < ZERO) - { - return std::numeric_limits::quiet_NaN(); - } - if (rhs == ZERO) - { - return ZERO; - } - - const RealT ratio = INITIALIZATION_TOLERANCE / rhs; - return std::sqrt(rhs) - * std::sqrt(TWO - / (std::hypot(ratio, TWO) + ratio)); + const RealT root_width = ONE / Math::MU; + const RealT knee = CURRENT_CIRCLE_KNEE; + + const ValueT absolute = std::abs(x); + const ValueT normalizer = absolute + knee; + const ValueT scaled_x = absolute / normalizer; + const ValueT scaled_knee = knee / normalizer; + const ValueT magnitude = normalizer + * std::sqrt(scaled_x * scaled_x + + scaled_knee * scaled_knee); + const ValueT half_sum = HALF * magnitude + + HALF * absolute; + const ValueT conjugate = QUARTER * knee * (knee / half_sum); + const ValueT hinged = HALF * x + + HALF * absolute + + conjugate; + // Divide through by MU to avoid overflow for large finite circles. + return hinged / std::sqrt(hinged + root_width * root_width); } /** - * @brief Compute off-axis capacity from a continuation state - * - * This expression must match the `ilcap` calculation during residual - * evaluation so initialization lands on the implemented model. + * @brief Compute an analytic seed for the sqrtramp() inverse * - * @param[in] ilmax Current-circle continuation state on the component base. - * @return Available off-axis current on the component base. + * @param[in] y Requested nonnegative value. + * @return Nonnegative inverse seed. + * @pre `y` is finite and nonnegative. */ template typename Reecb::RealT - Reecb::capacity(RealT ilmax) + Reecb::isqrtramp(RealT y) { - const RealT ilnorm = std::sqrt(ilmax * ilmax + INITIALIZATION_TOLERANCE); - return (ilmax / ilnorm) * ilmax; + if (y <= sqrtramp(ZERO)) + { + return ZERO; + } + + const RealT mu = Math::MU; + const RealT scaled_y = HALF * mu * y; + const RealT hinged = y / mu + * (scaled_y + std::hypot(scaled_y, ONE)); + const RealT square = hinged + - QUARTER * CURRENT_CIRCLE_KNEE + * (CURRENT_CIRCLE_KNEE / hinged); + return std::max(ZERO, square); } /** - * @brief Bisect an initial-limit interval to machine rounding - * - * The upper endpoint is returned because initialization needs the first - * not-below point; solveInitialLimit() separately validates that point as - * finite and feasible. A finite floating-point interval contains finitely - * many representable values, so the loop terminates when no interior - * midpoint remains. + * @brief Compute a difference of squares without overflowing its sum * - * @tparam FuncT Monotone predicate type. - * @param[in] a Lower endpoint, where `below` is true. - * @param[in] b Upper endpoint, where `below` is false. - * @param[in] below Predicate returning true only when finite reconstructed - * capacity lies below the requirement. - * @return The first representable upper-side endpoint. - * @pre `a` and `b` are finite, `a < b`, and `below` is monotone. + * @param[in] limit Total-current limit on the component base. + * @param[in] high Priority-axis current on the component base. + * @return @f$limit^2-high^2@f$. */ template - template - typename Reecb::RealT - Reecb::bisect(RealT a, RealT b, FuncT below) + template + [[gnu::always_inline]] inline ValueT + Reecb::circleSquare(RealT limit, ValueT high) { - while (true) - { - const RealT mid = std::midpoint(a, b); - if (mid <= a || b <= mid) - { - break; - } - - if (below(mid)) - { - a = mid; - } - else - { - b = mid; - } - } - - return b; + return ((HALF * limit - HALF * high) + * (HALF * limit + HALF * high)) + * FOUR; } /** - * @brief Solve the smallest initial total-current limit + * @brief Solve a feasible initial total-current limit * * @param[in] lower Lower bound for the component-base total-current limit. * @param[in] high Priority-axis current command on the component base. * @param[in] low Required off-axis capacity on the component base. - * @return The mutually consistent total limit, continuation state, and - * off-axis capacity, or `std::nullopt` when no finite feasible - * limit is found. + * @return The mutually consistent total limit and off-axis capacity, or + * `std::nullopt` when no finite feasible limit exists. * @pre The arguments are finite and nonnegative, with `lower >= high` and * `lower >= low`. * @warning This function contains conditional branching and may be used @@ -1180,75 +1160,48 @@ namespace GridKit auto Reecb::solveInitialLimit( RealT lower, RealT high, RealT low) -> std::optional { - const RealT ilmax = circleState(lower, high); - const RealT ilcap = capacity(ilmax); - if (!std::isfinite(ilcap)) + RealT imax = lower; + RealT cap = sqrtramp(circleSquare(imax, high)); + if (!std::isfinite(cap)) { return std::nullopt; } - if (!(ilcap < low)) + if (cap >= low) { - return InitialCurrentLimit{lower, ilmax, ilcap}; + return InitialCurrentLimit{imax, cap}; } - const auto below = [high, low](RealT limit) - { - const RealT state = circleState(limit, high); - const RealT cap = capacity(state); - return std::isfinite(cap) && cap < low; - }; - - // Invert ilcap = ilmax^2 / sqrt(ilmax^2 + tolerance), then recover - // imax from imax^2 = high^2 + ilmax * sqrt(ilmax^2 + tolerance). - const RealT delta = std::sqrt(INITIALIZATION_TOLERANCE); - const RealT ilreq = std::sqrt(low) - * std::sqrt(HALF - * (low + std::hypot(low, TWO * delta))); - const RealT seed = std::hypot( - high, std::sqrt(ilreq) * std::sqrt(std::hypot(ilreq, delta))); - - const RealT maximum = std::numeric_limits::max(); - RealT a = lower; - RealT b = seed; - if (!std::isfinite(b)) - { - b = maximum; - } - else - { - b = std::max(b, std::nextafter(a, maximum)); - } - if (!(a < b)) + const RealT square = isqrtramp(low); + if (!std::isfinite(square)) { return std::nullopt; } - while (below(b)) + imax = std::max(lower, std::hypot(high, std::sqrt(square))); + cap = sqrtramp(circleSquare(imax, high)); + + for (int step = 0; + step < std::numeric_limits::digits && cap < low; + ++step) { - a = b; - if (b >= maximum) + const RealT next = std::nextafter(imax, std::numeric_limits::max()); + if (!std::isfinite(next) || next <= imax) { return std::nullopt; } - if (b > maximum / TWO) + imax = next; + cap = sqrtramp(circleSquare(imax, high)); + if (!std::isfinite(cap)) { - b = maximum; - } - else - { - b *= TWO; + return std::nullopt; } } - const RealT result = bisect(a, b, below); - const RealT final_state = circleState(result, high); - const RealT final_cap = capacity(final_state); - if (!std::isfinite(result) || !std::isfinite(final_state) - || !std::isfinite(final_cap) || final_cap < low) + if (!std::isfinite(imax) || !std::isfinite(cap) || cap < low) { return std::nullopt; } - return InitialCurrentLimit{result, final_state, final_cap}; + return InitialCurrentLimit{imax, cap}; } /** diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md b/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md index 269d122cb..7eecfbd7a 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md @@ -3,29 +3,16 @@ REGCA is a first-generation WECC renewable generator/converter model for inverter-coupled resources. -The network equations solve $I_q^\mathrm{extra}$ for whatever nonnegative -current is needed to keep $V_T$ at or below $V_\mathrm{hv}^{\max}$. It is -approximately zero when the limit is inactive: - -```math -0 \le I_q^\mathrm{extra} - \perp \left(V_\mathrm{hv}^{\max} - V_T\right) \ge 0. -``` - -GridKit uses the smooth HVRCM form shown in the algebraic equations below. - ## Notes -- Internal current states and limiter quantities are on component base. -- Signal ports, monitor outputs, branch currents, and branch powers are on system base. -- LVACM uses $V_T$; LVPL uses $V_M$. -- PowerWorld fields `Qmin`, `Khv`, and `Xe` are accepted as optional DYD compatibility fields and are not used by the equations below. +None. ## Block Diagram -![](../../../../../docs/Figures/PhasorDynamics_REGCA_Diagram.png) +![REGCA generator/converter block diagram](../../../../../docs/Figures/PhasorDynamics_REGCA_Diagram.png) -Figure 1: Generator/Converter REGCA model. Figure courtesy of [PowerWorld](https://www.powerworld.com/WebHelp/) +Figure 1: REGCA generator/converter model. Figure courtesy of the +[PowerWorld REGC_A model reference](https://www.powerworld.com/WebHelp/Content/TransientModels_HTML/Machine%20Model%20REGC_A.htm). ## Model Parameters @@ -45,10 +32,10 @@ $V_{L0}$ | [p.u.] | `VL0` | LVPL zero-crossing volt $V_{L1}$ | [p.u.] | `VL1` | LVPL upper breakpoint voltage | 0.9 | Block name: `brkpt` $V_{A0}$ | [p.u.] | `VA0` | LVACM lower breakpoint voltage | 0.4 | Block name: `LVPnt0` $V_{A1}$ | [p.u.] | `VA1` | LVACM upper breakpoint voltage | 0.9 | Block name: `LVPnt1` -$V_\mathrm{hv}^{\max}$ | [p.u.] | `Vhvmax` | Terminal-voltage ceiling for HV reactive management | 1.2 | Block name: `VLim` -$Q^{\min}$ | [p.u.] | `Qmin` | Unused compatibility field | | Optional -$K_\mathrm{hv}$ | [p.u.] | `Khv` | Unused compatibility field | | Optional -$X_\mathrm{e}$ | [p.u.] | `Xe` | Unused compatibility field | | Optional +$V_\mathrm{hv}^{\max}$ | [p.u.] | `Vhvmax` | HV reactive management activation threshold | 1.2 | Block name: `VLim` +$Q^{\min}$ | [p.u.] | `Qmin` | PowerWorld compatibility field | | Optional; accepted but unused +$K_\mathrm{hv}$ | [p.u.] | `Khv` | HV reactive management gain | 0.7 | Optional; defaults to 0.7; block name: `Khv` +$X_\mathrm{e}$ | [p.u.] | `Xe` | PowerWorld compatibility field | | Optional; accepted but unused All listed JSON parameters are required unless marked optional. @@ -73,7 +60,9 @@ every other condition is a configuration error. 0 &\le V_{L0} < V_{L1} \\ 0 - &\le V_{A0} < V_{A1} < V_\mathrm{hv}^{\max} + &\le V_{A0} < V_{A1} < V_\mathrm{hv}^{\max} \\ + 0 \le K_\mathrm{hv} + &< \infty \end{aligned} ``` @@ -84,14 +73,14 @@ every other condition is a configuration error. s_L^\mathrm{off} &= 1 - s_L \\ k_\mathrm{base} - &= \dfrac{S^\mathrm{sys}}{S^\mathrm{base}} + &= \dfrac{S^\mathrm{sys}}{S^\mathrm{base}} \\ + K_L + &= 100 \end{aligned} ``` -Above the upper breakpoint the LVPL characteristic imposes no limit: the -release slope $K_L$ (default $100$ [p.u./p.u.], adjustable through -`setLvplGain`) renders the unlimited region with finite slope, and -the exact characteristic is the limit $K_L \to \infty$. +The fixed slope $K_L$ [p.u./p.u.] approximates the unbounded LVPL +characteristic above $V_{L1}$. ## Model Ports @@ -204,9 +193,7 @@ the sign that enables the corresponding limit. - V_\mathrm{r}(I_q - I_q^\mathrm{extra}) + V_\mathrm{i} I_p\,\text{linseg}(V_T; V_{A0}, V_{A1}, 1) \\ 0 &= -I_q^\mathrm{extra} - + \text{ramp}\!\left( - I_q^\mathrm{extra} - (V_\mathrm{hv}^{\max} - V_T) - \right) \\ + + K_\mathrm{hv}\,\text{ramp}(V_T - V_\mathrm{hv}^{\max}) \\ 0 &= -I_L + \text{linseg}(V_M; V_{L0}, V_{L1}, I_{L1}) + K_L\,\text{ramp}(V_M - V_{L1}) \\ @@ -244,10 +231,8 @@ CommonMath defines the [primitives](../../../../CommonMath.md#primitives) and ### Internal Initialization -REGCA requires $V_{A1} \le V_{T,0} < V_\mathrm{hv}^{\max}$. The lower bound -excludes initialization below the nominal upper LVACM breakpoint. The strict -upper bound is required because the smooth HVRCM constraint has no finite root -at or above the voltage limit. +REGCA requires $V_{A1} \le V_{T,0}$, which excludes initialization below the +nominal upper LVACM breakpoint. With LVPL enabled, REGCA additionally requires $I_{p,0} \le I_{L,0}$. Initialization rejects an operating point above the active-current integrator @@ -272,12 +257,7 @@ to zero: k_\mathrm{base} I_{p,0}^\mathrm{cmd} &= I_{p,0} \\ I_{q,0}^\mathrm{extra} - &\leftarrow \text{nonnegative solution of } - 0 = -I_{q,0}^\mathrm{extra} - + \text{ramp}\!\left( - I_{q,0}^\mathrm{extra} - - (V_\mathrm{hv}^{\max} - V_{T,0}) - \right) \\ + &= K_\mathrm{hv}\,\text{ramp}(V_{T,0} - V_\mathrm{hv}^{\max}) \\ I_{q,0}^\mathrm{cmd} &= \dfrac{Q_0}{V_{T,0}} + \dfrac{I_{q,0}^\mathrm{extra}}{k_\mathrm{base}} \\ @@ -286,12 +266,6 @@ to zero: \end{aligned} ``` -For the closed-form HVRCM solve, let -$x = \mu(V_\mathrm{hv}^{\max} - V_{T,0})$. The implementation switches at -$x = \log 2$ between two algebraically identical forms of -$\log(1 - \exp(-x))$. Their values and derivatives agree at the switch; the -split only avoids cancellation for small $x$. - The remaining algebraic quantities are then initialized as follows: ```math @@ -332,14 +306,14 @@ Output | Units | Description | Note - `validation()` checks construction, monitor creation, parameter validation, bus and signal configuration, and minimum time-constant handling. - `initializationAndSignals()` checks power-flow initialization, base conversion, signal publication, bus injection, and unattached-command latching. - `initializationDomain()` checks rejected and accepted voltage and LVPL initialization boundaries. -- `residualEquations()` checks every model residual against a fixed numerical answer key. +- `residualEquations()` checks every model residual at a hand-computable midpoint state. - `activeCurrentControl()` checks `rrpwr`, enabled and bypassed LVPL behavior, and tracking of a moving LVPL ceiling. - `reactiveCurrentControl()` checks the positive, negative, and unrestricted reactive-current recovery-rate branches. -- `highVoltageManagement()` checks HVRCM initialization, residual values through its transition, and its local derivative. +- `highVoltageManagement()` checks HVRCM initialization, $K_\mathrm{hv}$ loading, threshold behavior, and its local derivative. - `jacobian()` compares the dependency-tracking and Enzyme Jacobians for enabled and bypassed LVPL configurations when Enzyme support is enabled. -Because CommonMath limiters are smooth approximations, tests use smooth-equation -answer keys and reserve $100 \epsilon$ for roundoff checks. +Tests use model identities rather than frozen smoothing decimals and reserve +$100 \epsilon$ for roundoff checks. ## Appendix A: `rrpwr` diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp index 13b73c7f5..0ce153803 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp @@ -102,15 +102,14 @@ namespace GridKit Regca(BusT* bus, const ModelDataT& data); ~Regca(); - int setGridKitComponentID(IdxT component_id) override final; - int allocate() override final; - int verify() const override final; - int initialize() override final; - int tagDifferentiable() override final; - int setAbsoluteTolerance(RealT rel_tol) override final; - void setLvplGain(RealT KL); - int evaluateResidual() override final; - int evaluateJacobian() override final; + int setGridKitComponentID(IdxT component_id) override final; + int allocate() override final; + int verify() const override final; + int initialize() override final; + int tagDifferentiable() override final; + int setAbsoluteTolerance(RealT rel_tol) override final; + int evaluateResidual() override final; + int evaluateJacobian() override final; auto getSignals() -> ComponentSignals(1.0e-3); + static constexpr RealT KL_ = static_cast(100.0); static void logTimeConstantWarning(); BusT* bus_{nullptr}; @@ -225,7 +223,7 @@ namespace GridKit RealT VA0_{0}; RealT VA1_{0}; RealT Vhvmax_{0}; - RealT KL_{100.0}; + RealT Khv_{0.7}; IdxT parameter_error_count_{0}; diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp index 223635077..66138a218 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp @@ -14,8 +14,8 @@ namespace GridKit { namespace Converter { - /// Parameter keys for the REGCA converter model. All keys are - /// required except the optional PowerWorld compatibility fields. + /// Parameter keys for the REGCA converter model. Khv, Qmin, and Xe are + /// optional; all other keys are required. enum class RegcaParameters { p0, ///< Initial active power injection on system base @@ -32,10 +32,10 @@ namespace GridKit VL1, ///< LVPL upper breakpoint voltage VA0, ///< LVACM lower breakpoint voltage VA1, ///< LVACM upper breakpoint voltage - Vhvmax, ///< Terminal-voltage ceiling for HV reactive management + Vhvmax, ///< Terminal-voltage activation threshold for HV reactive management Qmin, ///< Optional PowerWorld compatibility field, accepted and unused - Khv, ///< Optional PowerWorld compatibility field, accepted and unused + Khv, ///< Optional HV reactive management gain, defaults to 0.7 Xe ///< Optional PowerWorld compatibility field, accepted and unused }; diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp index 7b1762b23..1fad805c1 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp @@ -7,6 +7,7 @@ #pragma once #include +#include #include #include @@ -177,55 +178,11 @@ namespace GridKit } /** - * @brief Calculate the initial HVRCM extra current. + * @brief Read parameters from the model data. * - * Solves the smooth HVRCM residual for \f$I_q^\mathrm{extra}\f$ given the - * strictly positive voltage margin \f$V_\mathrm{hv}^{\max} - V_T\f$. - * - * @param[in] dv Strictly positive voltage margin. - * @return Initial HVRCM extra current. - */ - template - scalar_type Regca::initialHvrcmCurrent( - scalar_type dv) const - { - static constexpr RealT log_two = std::numbers::ln2_v; - - const ScalarT x = Math::MU * dv; - - // Both branches evaluate log(1 - exp(-x)) and are algebraically - // identical, so their values and derivatives agree at x = log(2). - // The split only avoids cancellation for small x. - if (x < log_two) - { - return -(log_two - HALF * x - + std::log(std::sinh(HALF * x))) - / Math::MU; - } - - return -std::log1p(-std::exp(-x)) / Math::MU; - } - - /** - * @brief Set the LVPL release slope above the upper breakpoint. - * - * The finite slope is measured in p.u. current per p.u. voltage. The - * unlimited characteristic is recovered as the slope approaches infinity. - * - * @param[in] KL Positive LVPL release slope. - */ - template - void Regca::setLvplGain(RealT KL) - { - KL_ = KL; - } - - /** - * @brief Read the required parameters from the model data. - * - * Missing keys and invalid value types are counted for verify(). Integer - * JSON values are accepted for real parameters, and the LVPL switch accepts - * a boolean or the integer values zero and one. + * Missing required keys and invalid value types are counted for verify(). + * Integer JSON values are accepted for real parameters, and the LVPL switch + * accepts a boolean or the integer values zero and one. * * @param[in] data REGCA parameter data. */ @@ -302,6 +259,10 @@ namespace GridKit load_required_real(Params::VA0, VA0_, "VA0"); load_required_real(Params::VA1, VA1_, "VA1"); load_required_real(Params::Vhvmax, Vhvmax_, "Vhvmax"); + if (data.parameters.contains(Params::Khv)) + { + load_required_real(Params::Khv, Khv_, "Khv"); + } setDerivedParameters(); } @@ -456,7 +417,8 @@ namespace GridKit check(mva_base_ > ZERO, "mva must be positive"); check(Rpmax_ >= ZERO, "Rpmax must be non-negative"); check(IL1_ >= ZERO, "IL1 must be non-negative"); - check(KL_ > ZERO, "LVPL release slope must be positive"); + check(std::isfinite(Khv_) && Khv_ >= ZERO, + "Khv must be finite and non-negative"); check(ZERO <= VL0_ && VL0_ < VL1_, "VL0/VL1 must satisfy 0 <= VL0 < VL1"); check(ZERO <= VA0_ && VA0_ < VA1_ && VA1_ < Vhvmax_, "VA0/VA1/Vhvmax must satisfy 0 <= VA0 < VA1 < Vhvmax"); @@ -491,8 +453,8 @@ namespace GridKit * * @pre allocate() has completed, verify() reports no errors, and the * terminal bus has been initialized. - * @pre \f$V_{A1} \le V_{T,0} < V_\mathrm{hv}^{\max}\f$, and with LVPL - * enabled \f$I_{p,0} \le I_{L,0}\f$. + * @pre \f$V_{A1} \le V_{T,0}\f$, and with LVPL enabled + * \f$I_{p,0} \le I_{L,0}\f$. * @post All internal derivatives are zero. Unattached command ports retain * the resolved setpoints as constant commands. * @return Zero on success; nonzero when the operating point is rejected. @@ -527,12 +489,6 @@ namespace GridKit << "Regca: terminal voltage magnitude must be at least VA1 at initialization\n"; return 1; } - if (vt >= Vhvmax_) - { - Log::error() - << "Regca: terminal voltage magnitude must be below Vhvmax at initialization\n"; - return 1; - } // P0 is a system-base power-flow injection. Resolve the component-base // active current through the LVACM network-interface gain. @@ -548,11 +504,8 @@ namespace GridKit } const ScalarT ipcmd0 = ip0; - // Solve the smooth HVRCM constraint and preserve the requested Q0. The - // Vhvmax check above keeps the voltage margin strictly positive, so the - // solve is always finite. - const ScalarT dv = Vhvmax_ - vt; - const ScalarT iqextra0 = initialHvrcmCurrent(dv); + // Evaluate the HVRCM law and preserve the requested Q0. + const ScalarT iqextra0 = Khv_ * Math::ramp(vt - Vhvmax_); const ScalarT qnet0 = toComponentBase(static_cast(q0_) / vt); const ScalarT iqcmd0 = qnet0 + iqextra0; const ScalarT ir0 = (vi * qnet0 + vr * ip0 * lvacm) / vt; @@ -681,9 +634,9 @@ namespace GridKit + use_rqmin_ * (Math::max(fq, Rqmin_) - fq); const ScalarT fp_limited = rrpwr(ip, fp, Rpmax_); - // The LVPL ceiling IL = linseg(VM) moves with the sensed voltage; its - // rate is the exact chain rule (inside() is the linseg slope mask - // since ramp' = sigmoid), and a pinned Ip tracks the moving ceiling. + // The LVPL ceiling moves with sensed voltage. Its rate is the exact + // chain rule (inside() is the linseg slope mask and ramp' = sigmoid); + // a pinned Ip tracks the moving ceiling. const ScalarT vm_rate = (vt - vm) / TM_; const ScalarT il_rate = (IL1_ / (VL1_ - VL0_) * Math::inside(vm, VL0_, VL1_) + KL_ * Math::sigmoid(vm - VL1_)) @@ -698,7 +651,7 @@ namespace GridKit f[VT] = -vt * vt + vr * vr + vi * vi; f[IR] = -toComponentBase(vt * ir) + vi * qnet + vr * ip * lvacm; f[II] = -toComponentBase(vt * ii) - vr * qnet + vi * ip * lvacm; - f[IQEXTRA] = -iqextra + Math::ramp(iqextra - (Vhvmax_ - vt)); + f[IQEXTRA] = -iqextra + Khv_ * Math::ramp(vt - Vhvmax_); f[IL] = -il + Math::linseg(vm, VL0_, VL1_, IL1_) + KL_ * Math::ramp(vm - VL1_); f[PBR] = -pbr + vr * ir + vi * ii; diff --git a/tests/UnitTests/PhasorDynamics/ControllerReecbTests.hpp b/tests/UnitTests/PhasorDynamics/ControllerReecbTests.hpp index c036227d8..a4c5d34e9 100644 --- a/tests/UnitTests/PhasorDynamics/ControllerReecbTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ControllerReecbTests.hpp @@ -1,6 +1,8 @@ #pragma once +#include #include +#include #include #include #include @@ -75,7 +77,6 @@ namespace GridKit Vars::VPIQ, Vars::EPIV, Vars::RPORD, - Vars::ILMAX, Vars::ILCAP, Vars::IQMAX, Vars::IPMAX, @@ -288,32 +289,37 @@ namespace GridKit success *= fixture.initialize(kInitialIqcmd, kInitialIpcmd); success *= (fixture.evaluate() == 0); - const std::array initial_state{{ - {Vars::VMEAS, 1.0}, - {Vars::PMEAS, 1.5}, - {Vars::XPIQ, 0.0}, - {Vars::XPIV, 0.0}, - {Vars::QV, 1.5}, - {Vars::PORD, 1.5}, - {Vars::VT, 1.0}, - {Vars::VSAFE, 1.0}, - {Vars::SDIP, 1.0}, - {Vars::IQV, 0.0}, - {Vars::QREF, 1.5}, - {Vars::EQ, 0.0}, - {Vars::VPIQ, 0.5}, - {Vars::EPIV, 0.0}, - {Vars::RPORD, 0.0}, - {Vars::ILMAX, 2.0}, - {Vars::ILCAP, 2.0}, - {Vars::IQMAX, 2.0}, - {Vars::IPMAX, 2.5}, - {Vars::IQBASE, 0.0}, - {Vars::IQRAW, 1.5}, - {Vars::IQCMD, 0.75}, - {Vars::IPCMD, 0.75}, - }}; - success *= stateMatches(fixture.reecb, initial_state, "initialization"); + success *= stateMatches(fixture.reecb, + {{Vars::VMEAS, 1.0}, + {Vars::PMEAS, 1.5}, + {Vars::XPIQ, 0.0}, + {Vars::XPIV, 0.0}, + {Vars::QV, 1.5}, + {Vars::PORD, 1.5}, + {Vars::VT, 1.0}, + {Vars::VSAFE, 1.0}, + {Vars::SDIP, 1.0}, + {Vars::IQV, 0.0}, + {Vars::QREF, 1.5}, + {Vars::EQ, 0.0}, + {Vars::VPIQ, 0.5}, + {Vars::EPIV, 0.0}, + {Vars::RPORD, 0.0}, + {Vars::IPMAX, 2.5}, + {Vars::IQBASE, 0.0}, + {Vars::IQRAW, 1.5}, + {Vars::IQCMD, 0.75}, + {Vars::IPCMD, 0.75}}, + "initialization"); + success *= stateMatches(fixture.reecb, + {{Vars::ILCAP, 2.0}}, + "initial current-circle capacity", + kCircleTol); + + const auto* initial_values = fixture.reecb.y().getData(); + success *= scalarMatches(initial_values[index(Vars::IQMAX)], + initial_values[index(Vars::ILCAP)], + "initial reactive-current limit"); success *= scalarPreserved(fixture.iqcmd(), kInitialIqcmd, "preserved iqcmd"); success *= scalarPreserved(fixture.ipcmd(), kInitialIpcmd, "preserved ipcmd"); @@ -370,9 +376,12 @@ namespace GridKit success *= (system_base.evaluate() == 0); success *= stateMatches(system_base.reecb, {{Vars::PMEAS, 0.75}, - {Vars::PORD, 1.5}, - {Vars::ILMAX, 2.0}}, + {Vars::PORD, 1.5}}, "omitted component rating"); + success *= stateMatches(system_base.reecb, + {{Vars::ILCAP, 2.0}}, + "omitted-rating current-circle capacity", + kCircleTol); success *= allResidualsWithinInitTolerance(system_base.reecb); return success.report(__func__); @@ -422,7 +431,7 @@ namespace GridKit Fixture adjusted_imax(expanded_current); success *= adjusted_imax.initialize(0.75, 0.75); success *= (adjusted_imax.evaluate() == 0); - success *= stateMatches(adjusted_imax.reecb, {{Vars::ILMAX, 1.5}}, "adjusted Imax"); + success *= stateMatches(adjusted_imax.reecb, {{Vars::ILCAP, 1.5}}, "adjusted Imax"); success *= allResidualsWithinInitTolerance(adjusted_imax.reecb); auto reactive_pi = data; @@ -590,7 +599,7 @@ namespace GridKit const std::array active_cases{{ {1.0e-6, "near the lower active-current limit"}, {0.75, "interior active-current command"}, - {1.249999, "near the upper active-current limit"}, + {1.25 - 1.0e-6, "near the upper active-current limit"}, }}; // The recovered order limits are widened so the reconstruction, not @@ -659,8 +668,8 @@ namespace GridKit // The reactive command shares the inverse, at both signs. const std::array reactive_commands{{ - static_cast(0.999999), - static_cast(-0.999999), + static_cast(1.0 - 1.0e-6), + static_cast(-1.0 + 1.0e-6), }}; for (const RealT iqcmd : reactive_commands) { @@ -676,7 +685,7 @@ namespace GridKit bool p_priority; RealT iqcmd; RealT ipcmd; - RealT ilmax; + RealT capacity; const char* label; }; @@ -698,7 +707,10 @@ namespace GridKit success *= (boundary.evaluate() == 0); success *= scalarPreserved(boundary.iqcmd(), test_case.iqcmd, test_case.label); success *= scalarPreserved(boundary.ipcmd(), test_case.ipcmd, test_case.label); - success *= stateMatches(boundary.reecb, {{Vars::ILMAX, test_case.ilmax}}, test_case.label); + success *= stateMatches(boundary.reecb, + {{Vars::ILCAP, test_case.capacity}}, + test_case.label, + kCircleTol); success *= allResidualsWithinInitTolerance(boundary.reecb); } @@ -710,7 +722,7 @@ namespace GridKit success *= scalarPreserved(separated.iqcmd(), 5.0e-13, "scale-separated current command"); success *= allResidualsWithinInitTolerance(separated.reecb); - // A low configured Imax requires representable bisection to preserve + // A low configured Imax must still be expanded far enough to preserve // a strict low-priority command. auto capacity_data = exactness_data; capacity_data.parameters[Params::mva] = 100.0; @@ -728,14 +740,15 @@ namespace GridKit success *= (capacity_fixture.evaluate() == 0); success *= scalarPreserved(capacity_fixture.iqcmd(), iqcmd, "low-priority command"); success *= scalarPreserved(capacity_fixture.ipcmd(), ipcmd, "high-priority command"); - const RealT ilmax = static_cast(capacity_fixture.reecb.y().getData()[index(Vars::ILMAX)]); - const RealT ilcap = ilmax * ilmax - / std::sqrt(ilmax * ilmax + ReecbT::INITIALIZATION_TOLERANCE); + const RealT ilcap = static_cast(capacity_fixture.reecb.y().getData()[index(Vars::ILCAP)]); if (ilcap < iqcmd) { std::cout << "REECB low-priority capacity does not include its initial command\n"; success = false; } + success *= scalarMatches(ilcap, + iqcmd, + "minimal low-priority capacity"); success *= allResidualsWithinInitTolerance(capacity_fixture.reecb); } @@ -767,26 +780,30 @@ namespace GridKit success *= exhausted.initialize(0.0, 1.25); success *= (exhausted.evaluate() == 0); success *= scalarPreserved(exhausted.iqcmd(), 0.0, "exhausted reactive-current capacity"); - success *= stateMatches(exhausted.reecb, {{Vars::ILMAX, 0.0}}, "injection does not expand current circle"); - success *= allResidualsWithinInitTolerance(exhausted.reecb); + + success *= stateMatches(exhausted.reecb, + {{Vars::ILCAP, 0.0}}, + "injection does not expand current circle", + kCircleTol); + success *= allResidualsWithinInitTolerance(exhausted.reecb); Log::setVerbosity(previous_verbosity); return success.report(__func__); } - /// Check every residual row against an independent numerical answer key. - /// The expected values are literals, not a second implementation of REECB. + /// Check every residual row at a hand-computable state. TestOutcome residualEquations() { TestStatus success = true; Fixture fixture(makeResidualData(), kStateVr, kStateVi); fixture.attachAllInputs(); - setAnswerKeyInputs(fixture); + setResidualInputs(fixture); success *= fixture.prepare(0.25, 0.35); - setAnswerKeyState(fixture.reecb); + setResidualState(fixture.reecb); success *= (fixture.evaluate() == 0); + const RealT ideal_circle_leg = circleLeg(1.76); const std::array expected_residuals{{ {Vars::VMEAS, 0.99}, {Vars::PMEAS, 0.145}, @@ -803,8 +820,7 @@ namespace GridKit {Vars::VPIQ, 0.0}, {Vars::EPIV, 0.0}, {Vars::RPORD, 0.0}, - {Vars::ILMAX, 0.32}, - {Vars::ILCAP, 0.0}, + {Vars::ILCAP, ideal_circle_leg - 1.20}, {Vars::IQMAX, 0.0}, {Vars::IPMAX, 0.0}, {Vars::IQBASE, 0.0}, @@ -819,14 +835,27 @@ namespace GridKit { if (index(expected_residuals[row].variable) != row) { - std::cout << "REECB residual key position " << row << " names row " - << variableName(expected_residuals[row].variable) << '\n'; + Log::error() << "REECB residual expectation position " << row << " names row " + << variableName(expected_residuals[row].variable) << '\n'; success = false; } } - success *= residualsMatch(fixture.reecb, - expected_residuals, - "independent numerical answer key"); + const auto* residuals = fixture.reecb.getResidual().getData(); + for (const auto& [variable, expected] : expected_residuals) + { + RealT tolerance = kTol; + if (variable == Vars::ILCAP) + { + tolerance = kCircleTol; + } + + success *= variableMatches(residuals[index(variable)], + expected, + "residual", + variable, + "hand-computable state", + tolerance); + } return success.report(__func__); } @@ -882,7 +911,10 @@ namespace GridKit success *= allResidualsWithinInitTolerance(fixture.reecb); success *= scalarPreserved(fixture.iqcmd(), 0.75, "selector iqcmd"); success *= scalarPreserved(fixture.ipcmd(), 0.75, "selector ipcmd"); - success *= stateMatches(fixture.reecb, {{Vars::ILMAX, 2.0}}, "selector ILMAX"); + success *= stateMatches(fixture.reecb, + {{Vars::ILCAP, 2.0}}, + "selector ILCAP", + kCircleTol); // Exactly one reactive path carries the operating point. const auto* y = fixture.reecb.y().getData(); @@ -1210,7 +1242,6 @@ namespace GridKit setState(fixture.reecb, {{Vars::XPIV, test_case.state}, {Vars::EPIV, test_case.reference - 1.0}, - {Vars::ILMAX, 0.5}, {Vars::ILCAP, 0.5}, {Vars::IQMAX, 0.5}}); success *= (fixture.evaluate() == 0); @@ -1249,7 +1280,6 @@ namespace GridKit {Vars::IQV, 0.0}, {Vars::IQRAW, 0.0}, {Vars::IQCMD, 0.0}, - {Vars::ILMAX, 3.0}, {Vars::ILCAP, 3.0}, {Vars::IQMAX, 3.0}}); success *= (fixture.evaluate() == 0); @@ -1300,7 +1330,7 @@ namespace GridKit } /// Check the active-power ramp, its voltage gate and anti-windup, both - /// command limits, the priority circle, and the signed continuation. + /// command limits, and the priority current circle across closure. TestOutcome activeCurrentControl() { TestStatus success = true; @@ -1449,8 +1479,7 @@ namespace GridKit success *= fixture.prepare(0.0, 0.2); setControlState(fixture.reecb); setState(fixture.reecb, - {{Vars::ILMAX, 2.0}, - {Vars::ILCAP, 2.0}, + {{Vars::ILCAP, 2.0}, {Vars::IQMAX, 2.0}, {Vars::IQRAW, test_case.input}, {Vars::IQCMD, 0.0}, @@ -1480,8 +1509,7 @@ namespace GridKit success *= fixture.prepare(0.2, 0.0); setControlState(fixture.reecb); setState(fixture.reecb, - {{Vars::ILMAX, 2.0}, - {Vars::ILCAP, 2.0}, + {{Vars::ILCAP, 2.0}, {Vars::IPMAX, 2.0}, {Vars::IPCMD, 0.0}, {Vars::PORD, test_case.input}}); @@ -1494,19 +1522,20 @@ namespace GridKit { // The priority selector chooses which command consumes the circle. + const RealT capacity_state = 1.20; const std::array, 2> priority_cases{{ - {true, 0.32}, - {false, 0.56}, + {true, 1.76}, + {false, 2.00}, }}; - for (const auto& [p_priority, expected] : priority_cases) + for (const auto& [p_priority, square] : priority_cases) { auto data = makeResidualData(); data.parameters[Params::Pqflag] = p_priority; Fixture fixture(data, kStateVr, kStateVi); fixture.attachAllInputs(); - setAnswerKeyInputs(fixture); + setResidualInputs(fixture); success *= fixture.prepare(0.25, 0.35); - setAnswerKeyState(fixture.reecb); + setResidualState(fixture.reecb); success *= (fixture.evaluate() == 0); const char* label = "Q-priority current circle"; if (p_priority) @@ -1514,8 +1543,9 @@ namespace GridKit label = "P-priority current circle"; } success *= residualsMatch(fixture.reecb, - {{Vars::ILMAX, expected}}, - label); + {{Vars::ILCAP, circleLeg(square) - capacity_state}}, + label, + kCircleTol); } } @@ -1543,67 +1573,88 @@ namespace GridKit ipcmd = limit; } setState(fixture.reecb, - {{Vars::ILMAX, 0.0}, + {{Vars::ILCAP, 0.0}, {Vars::IQCMD, iqcmd}, {Vars::IPCMD, ipcmd}}); success *= (fixture.evaluate() == 0); success *= residualsMatch(fixture.reecb, - {{Vars::ILMAX, 0.0}}, + {{Vars::ILCAP, 0.0}}, "finite selected current circle"); success *= allResidualsFinite(fixture.reecb); } - } - { - // The signed-square continuation keeps a negative capacity iterate - // finite, and its magnitude still bounds the low-priority command. + const RealT open_limit = std::sqrt(maximum) / TWO; + auto open_data = makeData(); + open_data.parameters[Params::mva] = 100.0; + open_data.parameters[Params::Imax] = open_limit; + open_data.parameters[Params::Pqflag] = true; + Fixture open(open_data); + success *= open.prepare(0.0, 0.0); + setControlState(open.reecb); + setState(open.reecb, {{Vars::ILCAP, 0.0}, {Vars::IPCMD, 0.0}}); + success *= (open.evaluate() == 0); + success *= residualsMatch(open.reecb, + {{Vars::ILCAP, open_limit}}, + "finite open current circle", + kCircleTol); + success *= allResidualsFinite(open.reecb); + } + + { + // Probe the current circle above, at, and below closure. auto data = makeData(); data.parameters[Params::Imax] = 1.0; data.parameters[Params::Pqflag] = true; - const std::array continuation_cases{{ - {-0.5, 1.0}, - {0.5, 0.5}, - {0.0, 0.75}, + struct CircleCase + { + RealT ipcmd; + RealT square; + }; + + const std::array circle_cases{{ + {0.25, 0.75}, + {0.5, 0.0}, + {0.75, -1.25}, }}; - for (const auto& test_case : continuation_cases) + for (const auto& test_case : circle_cases) { + const RealT ideal_capacity = circleLeg(test_case.square); + Fixture fixture(data); - success *= fixture.prepare(0.0, 0.25); + success *= fixture.prepare(0.0, test_case.ipcmd); setControlState(fixture.reecb); setState(fixture.reecb, - {{Vars::ILMAX, test_case.input}, - {Vars::ILCAP, 0.0}, + {{Vars::ILCAP, 0.0}, {Vars::IQMAX, 0.0}, {Vars::IQRAW, 0.0}, - {Vars::IPCMD, 0.25}, + {Vars::IPCMD, test_case.ipcmd}, {Vars::IQCMD, 0.0}, {Vars::QV, 1.0}}); success *= (fixture.evaluate() == 0); success *= residualsMatch(fixture.reecb, - {{Vars::ILMAX, test_case.expected}}, - "signed capacity continuation"); + {{Vars::ILCAP, ideal_capacity}}, + "off-axis capacity", + kCircleTol); - RealT expected_capacity = 0.0; - if (test_case.input != ZERO) + const RealT capacity = fixture.reecb.getResidual().getData()[index(Vars::ILCAP)]; + if (!std::isfinite(capacity) || capacity < ZERO) { - expected_capacity = 0.5; + Log::error() << "REECB current circle produced an invalid capacity\n"; + success = false; } - success *= residualsMatch(fixture.reecb, - {{Vars::ILCAP, expected_capacity}}, - "signed off-axis capacity"); - setState(fixture.reecb, {{Vars::ILCAP, expected_capacity}}); + setState(fixture.reecb, {{Vars::ILCAP, capacity}}); success *= (fixture.evaluate() == 0); success *= residualsMatch(fixture.reecb, - {{Vars::IQMAX, expected_capacity}}, + {{Vars::IQMAX, capacity}}, "reactive-current capacity"); setState(fixture.reecb, - {{Vars::IQMAX, expected_capacity}, {Vars::IQRAW, 1.0}}); + {{Vars::IQMAX, capacity}, {Vars::IQRAW, 1.0}}); success *= (fixture.evaluate() == 0); success *= residualsMatch(fixture.reecb, - {{Vars::IQCMD, expected_capacity}}, + {{Vars::IQCMD, capacity}}, "capacity magnitude bound"); success *= allResidualsFinite(fixture.reecb); } @@ -1661,29 +1712,50 @@ namespace GridKit { const char* label; bool p_priority; - RealT ilmax; + RealT imax; + RealT capacity; }; - const std::array current_circle_probes{{ - {"negative signed current-circle capacity", true, -2.0}, - {"zero current-circle capacity", false, 0.0}, + const std::array current_circle_probes{{ + {"open current circle", true, 2.5, 2.0}, + {"closing current circle", true, 0.41, 0.05}, + {"closed current circle", true, 0.4, 0.0}, + {"over-driven current circle", true, 0.2, 0.0}, + {"over-driven current circle on the reactive axis", false, 0.1, 0.0}, }}; for (const auto& probe : current_circle_probes) { auto data = makeJacobianData(); data.parameters[Params::Pqflag] = probe.p_priority; + data.parameters[Params::Imax] = probe.imax; success *= jacobiansMatch( dependencyTrackingJacobian( - data, kNonunitAlpha, success, probe.ilmax), - enzymeJacobian(data, kNonunitAlpha, success, probe.ilmax), + data, kNonunitAlpha, success, probe.capacity), + enzymeJacobian(data, kNonunitAlpha, success, probe.capacity), probe.label); } + auto transition_data = makeJacobianData(); + transition_data.parameters[Params::Pqflag] = true; + transition_data.parameters[Params::Imax] = 0.4; + success *= jacobiansMatch( + dependencyTrackingJacobian( + transition_data, kNonunitAlpha, success, 0.1, 0.2, 0.1), + enzymeJacobian( + transition_data, kNonunitAlpha, success, 0.1, 0.2, 0.1), + "open current-circle transition"); + success *= jacobiansMatch( + dependencyTrackingJacobian( + transition_data, kNonunitAlpha, success, 0.1, 0.2, 0.3), + enzymeJacobian( + transition_data, kNonunitAlpha, success, 0.1, 0.2, 0.3), + "over-driven current-circle transition"); + struct NonlinearProbe { const char* label; RealT mva; - RealT ilmax; + RealT capacity; RealT epiv; }; @@ -1701,8 +1773,8 @@ namespace GridKit data.parameters[Params::mva] = probe.mva; success *= jacobiansMatch( dependencyTrackingJacobian( - data, kNonunitAlpha, success, probe.ilmax, probe.epiv), - enzymeJacobian(data, kNonunitAlpha, success, probe.ilmax, probe.epiv), + data, kNonunitAlpha, success, probe.capacity, probe.epiv), + enzymeJacobian(data, kNonunitAlpha, success, probe.capacity, probe.epiv), probe.label); } @@ -1898,6 +1970,16 @@ namespace GridKit inline static const RealT kUnitSlopeAngle = std::atan(ONE); + /// Exact off-axis leg of the ideal current circle. + static RealT circleLeg(RealT square) + { + return std::sqrt(std::max(square, ZERO)); + } + + /// MU-aware tolerance for ideal current-circle comparisons. + static constexpr RealT kCircleTol = + std::max(ONE / (Math::MU * Math::MU), kTol); + static constexpr size_t kBusVrColumn = index(Vars::MAXIMUM); Data makeMinimalData() const @@ -2048,7 +2130,7 @@ namespace GridKit } template - void setAnswerKeyInputs(Fixture& fixture) const + void setResidualInputs(Fixture& fixture) const { fixture.input(Ext::PE) = static_cast(0.3); fixture.input(Ext::QGEN) = static_cast(-0.1); @@ -2057,11 +2139,10 @@ namespace GridKit fixture.input(Ext::PREF) = static_cast(0.325); } - /// The rich state shared by the residual answer key and the priority - /// circle. Every smooth-transition argument keeps a saturation margin, - /// so each row carries its ideal value. + /// State shared by the residual and priority-circle checks. Each smooth + /// transition retains enough margin to take its ideal value. template - void setAnswerKeyState(PhasorDynamics::Controller::Reecb& reecb) const + void setResidualState(PhasorDynamics::Controller::Reecb& reecb) const { setState(reecb, {{Vars::VMEAS, 0.80}, @@ -2079,7 +2160,6 @@ namespace GridKit {Vars::VPIQ, 1.00}, {Vars::EPIV, 0.20}, {Vars::RPORD, 0.20}, - {Vars::ILMAX, 1.20}, {Vars::ILCAP, 1.20}, {Vars::IQMAX, 1.20}, {Vars::IPMAX, 1.50}, @@ -2118,7 +2198,6 @@ namespace GridKit {Vars::VPIQ, 0.0}, {Vars::EPIV, 0.0}, {Vars::RPORD, 0.0}, - {Vars::ILMAX, 1.4}, {Vars::ILCAP, 1.4}, {Vars::IQMAX, 1.4}, {Vars::IPMAX, 1.5}, @@ -2133,8 +2212,9 @@ namespace GridKit /// probe. The optional states place selected smooth branches explicitly. template void setJacobianState(Fixture& fixture, - RealT ilmax, - RealT epiv = static_cast(0.2)) const + RealT capacity, + RealT epiv = static_cast(0.2), + RealT ipcmd = static_cast(0.2)) const { fixture.input(Ext::PE) = static_cast(0.25); fixture.input(Ext::QGEN) = static_cast(0.5); @@ -2159,14 +2239,13 @@ namespace GridKit {Vars::VPIQ, 1.0}, {Vars::EPIV, epiv}, {Vars::RPORD, 0.05}, - {Vars::ILMAX, ilmax}, - {Vars::ILCAP, std::abs(ilmax)}, - {Vars::IQMAX, std::abs(ilmax)}, - {Vars::IPMAX, std::abs(ilmax)}, + {Vars::ILCAP, capacity}, + {Vars::IQMAX, capacity}, + {Vars::IPMAX, capacity}, {Vars::IQBASE, 0.2}, {Vars::IQRAW, 0.2}, {Vars::IQCMD, 0.1}, - {Vars::IPCMD, 0.2}}); + {Vars::IPCMD, ipcmd}}); fixture.reecb.yp().setDataUpdated(); } @@ -2222,10 +2301,10 @@ namespace GridKit } } - setAnswerKeyInputs(implicit_defaults); - setAnswerKeyInputs(explicit_defaults); - setAnswerKeyState(implicit_defaults.reecb); - setAnswerKeyState(explicit_defaults.reecb); + setResidualInputs(implicit_defaults); + setResidualInputs(explicit_defaults); + setResidualState(implicit_defaults.reecb); + setResidualState(explicit_defaults.reecb); if (implicit_defaults.evaluate() != 0 || explicit_defaults.evaluate() != 0) { success = false; @@ -2383,7 +2462,6 @@ namespace GridKit "VPIQ", "EPIV", "RPORD", - "ILMAX", "ILCAP", "IQMAX", "IPMAX", @@ -2496,7 +2574,8 @@ namespace GridKit bool rowsMatch(const VectorT& vector, const ValuesT& values, const char* what, - const char* context) const + const char* context, + RealT tolerance = kTol) const { bool success = true; const auto* vector_values = vector.getData(); @@ -2506,7 +2585,8 @@ namespace GridKit expected, what, variable, - context)) + context, + tolerance)) { success = false; } @@ -2516,32 +2596,36 @@ namespace GridKit bool residualsMatch(const ReecbT& reecb, std::initializer_list values, - const char* context = "") const + const char* context = "", + RealT tolerance = kTol) const { - return rowsMatch(reecb.getResidual(), values, "residual", context); + return rowsMatch(reecb.getResidual(), values, "residual", context, tolerance); } template bool residualsMatch(const ReecbT& reecb, const std::array& values, - const char* context = "") const + const char* context = "", + RealT tolerance = kTol) const { - return rowsMatch(reecb.getResidual(), values, "residual", context); + return rowsMatch(reecb.getResidual(), values, "residual", context, tolerance); } bool stateMatches(const ReecbT& reecb, std::initializer_list values, - const char* context = "") const + const char* context = "", + RealT tolerance = kTol) const { - return rowsMatch(reecb.y(), values, "state", context); + return rowsMatch(reecb.y(), values, "state", context, tolerance); } template bool stateMatches(const ReecbT& reecb, const std::array& values, - const char* context = "") const + const char* context = "", + RealT tolerance = kTol) const { - return rowsMatch(reecb.y(), values, "state", context); + return rowsMatch(reecb.y(), values, "state", context, tolerance); } bool allResidualsWithinInitTolerance(const ReecbT& reecb) const @@ -2711,15 +2795,16 @@ namespace GridKit dependencyTrackingJacobian(const Data& data, RealT alpha, TestStatus& success, - RealT ilmax = 2.0, - RealT epiv = static_cast(0.2)) const + RealT capacity = 2.0, + RealT epiv = static_cast(0.2), + RealT ipcmd = static_cast(0.2)) const { using DepVar = DependencyTracking::Variable; Fixture fixture(data, kStateVr, kStateVi); fixture.attachAllInputs(); success *= fixture.prepare(0.0, 0.2); - setJacobianState(fixture, ilmax, epiv); + setJacobianState(fixture, capacity, epiv, ipcmd); numberVariables(fixture, alpha); success *= (fixture.evaluate() == 0); @@ -2737,8 +2822,9 @@ namespace GridKit enzymeJacobian(const Data& data, RealT alpha, TestStatus& success, - RealT ilmax = 2.0, - RealT epiv = static_cast(0.2)) const + RealT capacity = 2.0, + RealT epiv = static_cast(0.2), + RealT ipcmd = static_cast(0.2)) const { Fixture fixture(data, kStateVr, kStateVi); fixture.attachAllInputs(); @@ -2749,7 +2835,7 @@ namespace GridKit fixture.bus.setVariableIndex(row, fixture.reecb.size() + row); } - setJacobianState(fixture, ilmax, epiv); + setJacobianState(fixture, capacity, epiv, ipcmd); fixture.reecb.updateTime(0.0, alpha); success *= (fixture.evaluate() == 0); success *= (fixture.reecb.evaluateJacobian() == 0); diff --git a/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp index d6f3537a4..d2524f824 100644 --- a/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include @@ -37,6 +38,8 @@ namespace GridKit static constexpr ScalarT kTol = static_cast(100.0) * std::numeric_limits::epsilon(); + static constexpr RealT kSmoothTol = + std::max(ONE / (Math::MU * Math::MU), kTol); /// Construction, the monitor, and every verify() error class: missing /// and invalid parameters, a null bus, and an unlinked command port. @@ -79,9 +82,15 @@ namespace GridKit success *= invalidParameterCase(bus, Params::VL1, 0.3); success *= invalidParameterCase(bus, Params::VA1, 0.3); - // Vhvmax must lie strictly above VA1: rejected at VA1 and just below. + // Vhvmax must lie strictly above VA1: reject equality and reversed ordering. success *= invalidParameterCase(bus, Params::Vhvmax, kVa1); - success *= invalidParameterCase(bus, Params::Vhvmax, kJustBelowVa1); + success *= invalidParameterCase(bus, Params::Vhvmax, HALF * kVa1); + + // Khv is optional, finite, and non-negative. + success *= invalidParameterCase(bus, Params::Khv, -0.1); + success *= invalidParameterCase(bus, + Params::Khv, + std::numeric_limits::infinity()); // A null bus and an attached command with no linked source count as // configuration errors on the same footing as bad parameters. @@ -134,21 +143,19 @@ namespace GridKit success *= fixture.initialize(); success *= (fixture.evaluate() == 0); - // Fixed answer key for Vr = 0.8, Vi = 0.6, P0 = 0.4, Q0 = -0.1, and - // a 50 MVA component base: IP = 2 P0 / VT and IQ = 2 Q0 / VT on the - // component base, with IR and II returned to the system base. + // P0/Q0 and branch currents remain on system base; current-command + // signals convert the 50 MVA component-base states back to system base. const auto* y = fixture.regca.y().getData(); success *= scalarMatches(y[index(Vars::VM)], 1.0, "VM"); success *= scalarMatches(y[index(Vars::VT)], 1.0, "VT"); - success *= scalarMatches(y[index(Vars::IP)], 0.80000000000025173, "IP"); - success *= scalarMatches(y[index(Vars::IQ)], -0.2, "IQ"); + success *= scalarMatches(y[index(Vars::IQ)] - y[index(Vars::IQEXTRA)], -0.2, "net reactive current"); success *= scalarMatches(y[index(Vars::IR)], 0.26, "IR"); success *= scalarMatches(y[index(Vars::II)], 0.32, "II"); success *= scalarMatches(y[index(Vars::PBR)], 0.4, "PBR"); success *= scalarMatches(y[index(Vars::QBR)], -0.1, "QBR"); - success *= scalarMatches(fixture.ipcmd, 0.40000000000012587, "published ipcmd"); - success *= scalarMatches(fixture.iqcmd, -0.1, "published iqcmd"); + success *= scalarMatches(fixture.ipcmd, 0.5 * y[index(Vars::IP)], "published ipcmd base conversion"); + success *= scalarMatches(fixture.iqcmd, 0.5 * y[index(Vars::IQ)], "published iqcmd base conversion"); success *= scalarMatches(pbranch_node.read(), 0.4, "pbranch signal"); // The accumulated bus injection is the system-base branch current. @@ -165,34 +172,32 @@ namespace GridKit Fixture latched(latch_data); success *= latched.initialize(); + success *= (latched.evaluate() == 0); + success *= allResidualsZero(latched.regca); auto* latched_y = latched.regca.y().getData(); - success *= scalarMatches(latched_y[index(Vars::IP)], - 0.60000000000018872, - "initialized IP"); - success *= scalarMatches(latched_y[index(Vars::IQ)], 0.2, "initialized IQ"); + success *= scalarMatches(latched_y[index(Vars::PBR)], 0.6, "latched PBR"); + success *= scalarMatches(latched_y[index(Vars::QBR)], 0.2, "latched QBR"); - latched_y[index(Vars::IP)] = 0.5; - latched_y[index(Vars::IQ)] = 0.14; + latched_y[index(Vars::IP)] -= 0.1; // arbitrary current displacement + latched_y[index(Vars::IQ)] -= 0.06; // arbitrary current displacement latched.regca.y().setDataUpdated(); success *= (latched.evaluate() == 0); - // f[IP] = (0.6 - 0.5) / Tg and f[IQ] = (0.2 - 0.14) / Tg with - // Tg = 0.2; both rates sit inside every limiter. + // The latched commands restore both displaced states at their ideal + // interior first-order rates. const auto* f = latched.regca.getResidual().getData(); - success *= scalarMatches(f[index(Vars::IP)], - 0.50000000000094358, - "latched active-current rate"); - success *= scalarMatches(f[index(Vars::IQ)], 0.3, "latched reactive-current rate"); + success *= scalarMatches(f[index(Vars::IP)], 0.5, "latched active-current rate", kSmoothTol); + success *= scalarMatches(f[index(Vars::IQ)], 0.3, "latched reactive-current rate", kSmoothTol); return success.report(__func__); } - /// The admissible initialization domain: every rejected operating - /// point, then the accepted boundaries next to them. + /// Rejected and representative admissible initialization points. TestOutcome initializationDomain() { - TestStatus success = true; + TestStatus success = true; + const RealT vm = 0.65; const auto previous_verbosity = Log::verbosity(); // Suppress expected errors from the inadmissible initialization cases below. @@ -203,26 +208,16 @@ namespace GridKit { const char* label; RealT vr; - RealT p0; - RealT il1; }; - // P0 and IL1 default to the makeData() values; each row breaks - // exactly one initialize() guard. - const std::array rejected{{ - {"terminal voltage at Vhvmax", kHvrcmVoltageLimit, 0.0, 1.1}, - {"terminal voltage above Vhvmax", kHvrcmVoltageLimit + 0.1, 0.0, 1.1}, - {"terminal voltage just below VA1", kJustBelowVa1, 0.0, 1.1}, - {"zero terminal voltage", 0.0, 0.0, 1.1}, + const std::array rejected{{ + {"terminal voltage below VA1", HALF * kVa1}, + {"zero terminal voltage", 0.0}, }}; for (const auto& test_case : rejected) { - auto data = makeData(); - data.parameters[Params::p0] = test_case.p0; - data.parameters[Params::IL1] = test_case.il1; - - Fixture fixture(data, test_case.vr); + Fixture fixture(makeData(), test_case.vr); success *= fixture.prepare(); if (fixture.regca.initialize() == 0) { @@ -266,99 +261,93 @@ namespace GridKit success *= scalarMatches(y[index(Vars::QBR)], 0.1, "QBR at the LVACM upper breakpoint"); } - // A ceiling just above the requested current is admissible. + // Inside the LVPL segment, initialization accepts an active current + // strictly below the moving limit. { auto data = makeData(); - data.parameters[Params::p0] = 0.6; - data.parameters[Params::IL1] = 0.61; - - Fixture fixture(data); - success *= fixture.initialize(); - success *= (fixture.evaluate() == 0); - success *= allResidualsZero(fixture.regca); - - const auto* y = fixture.regca.y().getData(); - success *= scalarMatches(y[index(Vars::IP)], - 0.60000000000018872, - "IP below the LVPL ceiling"); - success *= scalarMatches(y[index(Vars::PBR)], 0.6, "PBR below the LVPL ceiling"); - } - - // The active-current state may initialize exactly on the LVPL ceiling. - // VA1 is lowered so the terminal voltage sits inside the LVPL ramp, - // where the release term vanishes. - { - auto data = makeData(); - data.parameters[Params::IL1] = 0.0; + data.parameters[Params::p0] = 0.25; + data.parameters[Params::IL1] = 1.0; data.parameters[Params::VA1] = 0.5; - Fixture fixture(data, 0.6); + Fixture fixture(data, vm); success *= fixture.initialize(); success *= (fixture.evaluate() == 0); success *= allResidualsZero(fixture.regca); - const auto* y = fixture.regca.y().getData(); - success *= scalarMatches(y[index(Vars::IP)], 0.0, "IP at the LVPL ceiling"); - success *= scalarMatches(y[index(Vars::IL)], 0.0, "LVPL ceiling at IP"); + const auto* y = fixture.regca.y().getData(); + if (!(y[index(Vars::IP)] < y[index(Vars::IL)])) + { + std::cout << "IP is not below the active LVPL limit\n"; + success = false; + } + success *= scalarMatches(y[index(Vars::PBR)], 0.25, "PBR below the LVPL limit"); } // Above the upper breakpoint the ceiling releases: an operating point // beyond IL1 initializes at healthy voltage. { - auto data = makeData(); - data.parameters[Params::p0] = 1.3; + const RealT il1 = 1.1; + auto data = makeData(); + data.parameters[Params::p0] = 1.3; + data.parameters[Params::IL1] = il1; Fixture fixture(data); success *= fixture.initialize(); success *= (fixture.evaluate() == 0); success *= allResidualsZero(fixture.regca); - const auto* y = fixture.regca.y().getData(); - success *= scalarMatches(y[index(Vars::IP)], - 1.3000000000004091, - "IP beyond IL1 with the ceiling released"); - success *= scalarMatches(y[index(Vars::PBR)], 1.3, "PBR beyond IL1 with the ceiling released"); + const auto* y = fixture.regca.y().getData(); + if (!(y[index(Vars::IP)] > il1)) + { + std::cout << "IP did not exceed IL1 in the released region\n"; + success = false; + } + success *= scalarMatches(y[index(Vars::PBR)], 1.3, "PBR beyond IL1"); } - // With LVPL bypassed a collapsed ceiling does not constrain the - // initialization. + // Bypassing LVPL admits the same midpoint operating point even when + // its active current exceeds the moving limit. { auto data = makeData(); - data.parameters[Params::p0] = 0.6; + data.parameters[Params::p0] = 0.25; data.parameters[Params::sL] = false; data.parameters[Params::IL1] = 0.2; + data.parameters[Params::VA1] = 0.5; - Fixture fixture(data); + Fixture fixture(data, vm); success *= fixture.initialize(); success *= (fixture.evaluate() == 0); success *= allResidualsZero(fixture.regca); - const auto* y = fixture.regca.y().getData(); - success *= scalarMatches(y[index(Vars::IP)], - 0.60000000000018872, - "IP with LVPL bypassed"); + const auto* y = fixture.regca.y().getData(); + if (!(y[index(Vars::IP)] > y[index(Vars::IL)])) + { + std::cout << "Bypassed IP did not exceed the LVPL limit\n"; + success = false; + } + success *= scalarMatches(y[index(Vars::PBR)], 0.25, "PBR with LVPL bypassed"); } Log::setVerbosity(previous_verbosity); return success.report(__func__); } - /// The complete equation answer key at a hand-computed state. + /// Check every residual row at a hand-computable midpoint state. TestOutcome residualEquations() { TestStatus success = true; Fixture fixture(makeDynamicData(), kStateVr, kStateVi); - fixture.attachIpcmd(kStateIpcmd); + fixture.attachIpcmd(kStateIp); fixture.attachIqcmd(kStateIqcmd); success *= fixture.initialize(); // initialize() published steady-state commands over the attached - // values; restore the commands the answer key assumes. - fixture.ipcmd = kStateIpcmd; + // values; restore the commands used by the residual state. + fixture.ipcmd = kStateIp; fixture.iqcmd = kStateIqcmd; - setAnswerKeyState(fixture.regca); + setResidualState(fixture.regca); success *= (fixture.evaluate() == 0); struct ExpectedResidual @@ -368,19 +357,19 @@ namespace GridKit RealT value; }; - // Each entry is the expected value of the named README equation at - // the answer-key state. + // The LVACM and LVPL states use the midpoint of their breakpoint + // interval, where linseg is exactly one half for every MU. const std::array expected{{ - {Vars::VM, "VM", 0.865}, // -VM' + (VT - VM) / TM - {Vars::IQ, "IQ", 1.52}, // -IQ' + max(fq, Rqmin) - {Vars::IP, "IP", 0.22}, // -IP' + IL' + awmax(IP - IL, fp_limited - IL', 0) - {Vars::VT, "VT", -0.035}, // -VT^2 + Vr^2 + Vi^2 - {Vars::IR, "IR", 0.24999999999974593}, // -VT*IR + Vi*(IQ - IQEXTRA) + Vr*IP*linseg(VT) - {Vars::II, "II", 0.25099999999993317}, // -VT*II - Vr*(IQ - IQEXTRA) + Vi*IP*linseg(VT) - {Vars::IQEXTRA, "IQEXTRA", -0.03}, // smooth HVRCM constraint - {Vars::IL, "IL", 0.35}, // -IL + linseg(VM, VL0, VL1, IL1) - {Vars::PBR, "PBR", 0.0}, // -PBR + Vr*IR + Vi*II - {Vars::QBR, "QBR", 0.0}, // -QBR + Vi*IR - Vr*II + {Vars::VM, "VM", -0.01}, + {Vars::IQ, "IQ", 1.52}, + {Vars::IP, "IP", -0.03}, + {Vars::VT, "VT", 0.5425}, + {Vars::IR, "IR", -0.2875}, + {Vars::II, "II", 0.1265}, + {Vars::IQEXTRA, "IQEXTRA", -0.03}, + {Vars::IL, "IL", 0.35}, + {Vars::PBR, "PBR", 0.0}, + {Vars::QBR, "QBR", 0.0}, }}; const auto& residual = fixture.regca.getResidual(); @@ -389,7 +378,7 @@ namespace GridKit const auto* f = residual.getData(); for (const auto& row : expected) { - success *= scalarMatches(f[index(row.variable)], row.value, row.name); + success *= scalarMatches(f[index(row.variable)], row.value, row.name, kSmoothTol); } return success.report(__func__); @@ -580,19 +569,15 @@ namespace GridKit return success.report(__func__); } - /// High-voltage reactive current management: the initialization root, - /// the residual across the transition, and its local derivative. + /// High-voltage reactive current management: initialization, gain, + /// threshold behavior, and the local derivative. TestOutcome highVoltageManagement() { TestStatus success = true; - // Initialization near Vhvmax solves the smooth constraint for a - // nonzero IQEXTRA. Half the transition margin keeps the root - // distinct from the margin, so a swapped-argument residual cannot - // satisfy this operating point. + // Initialization above Vhvmax activates HVRCM and preserves Q0. { - const RealT margin = 0.5 * hvrcmTransition(); - const RealT terminal_voltage = kHvrcmVoltageLimit - margin; + const RealT terminal_voltage = kHvrcmVoltageLimit + kHvrcmOffset; auto data = makeData(); data.parameters[Params::q0] = 0.1; @@ -602,64 +587,61 @@ namespace GridKit success *= (fixture.evaluate() == 0); success *= allResidualsZero(fixture.regca); - // The zero residual already pins the root to the HVRCM equation. - // What remains to check is that it is a genuine compensation - // above the margin it cancels, and that it leaves Q0 intact. - const auto* y = fixture.regca.y().getData(); - if (!(y[index(Vars::IQEXTRA)] > margin)) + const auto* y = fixture.regca.y().getData(); + const RealT extra_current = y[index(Vars::IQEXTRA)]; + success *= scalarMatches( + extra_current, kHvrcmGain * kHvrcmOffset, "IQEXTRA with the default Khv", kSmoothTol); + success *= scalarMatches( + y[index(Vars::IQ)] - y[index(Vars::IQEXTRA)], 0.1 / terminal_voltage, "IQ preserves Q0 after HVRCM compensation"); + success *= scalarMatches(y[index(Vars::QBR)], 0.1, "QBR"); + } + + // At the activation threshold, the analytical ramp value is ln(2)/MU. + // This also checks that a supplied Khv replaces the default. + { + for (const RealT gain : {0.0, 0.5}) { - std::cout << "IQEXTRA is not above the voltage margin it compensates\n"; - success = false; + auto data = makeData(); + data.parameters[Params::Khv] = gain; + + Fixture fixture(data, kHvrcmVoltageLimit); + success *= fixture.initialize(); + success *= (fixture.evaluate() == 0); + success *= allResidualsZero(fixture.regca); + + const auto* y = fixture.regca.y().getData(); + success *= scalarMatches(y[index(Vars::IQEXTRA)], gain * hvrcmTransition(), "IQEXTRA at a supplied Khv"); } - success *= scalarMatches(y[index(Vars::IQ)] - y[index(Vars::IQEXTRA)], - 0.1 / terminal_voltage, - "IQ preserves Q0 after HVRCM compensation"); - success *= scalarMatches(y[index(Vars::QBR)], 0.1, "QBR"); } - // Numeric answer keys pin both sign and magnitude of the HVRCM row. + // The ramp identity ramp(x) - ramp(-x) = x pins the orientation and + // gain without encoding a MU-specific decimal. { Fixture fixture(makeData()); success *= fixture.initialize(); - struct HvrcmCase + auto* y = fixture.regca.y().getData(); + auto residualAt = [&](RealT voltage, RealT extra_current) { - RealT voltage; - RealT extra_current; - RealT expected_residual; + y[index(Vars::VT)] = voltage; + y[index(Vars::IQEXTRA)] = extra_current; + fixture.regca.y().setDataUpdated(); + success *= (fixture.evaluate() == 0); + return fixture.regca.getResidual().getData()[index(Vars::IQEXTRA)]; }; - const RealT transition = hvrcmTransition(); - - // Away from the transition the smooth tail is below the tolerance - // and the row reduces to -IQEXTRA; at and beyond the ceiling the - // tail is the whole answer. - const std::array cases{{ - {kHvrcmVoltageLimit - 0.2, 0.05, -0.05}, - {kHvrcmVoltageLimit - 0.1, -0.05, 0.05}, - {kHvrcmVoltageLimit - transition, transition, 0.0}, - {kHvrcmVoltageLimit, 0.0, transition}, - {kHvrcmVoltageLimit + 0.1, 0.1, 0.1}, - }}; + const RealT below = residualAt(kHvrcmVoltageLimit - kHvrcmOffset, 0.0); + const RealT at = residualAt(kHvrcmVoltageLimit, 0.0); + const RealT above = residualAt(kHvrcmVoltageLimit + kHvrcmOffset, 0.0); - auto* y = fixture.regca.y().getData(); - for (const auto& test_case : cases) - { - y[index(Vars::VT)] = test_case.voltage; - y[index(Vars::IQEXTRA)] = test_case.extra_current; - fixture.regca.y().setDataUpdated(); + success *= scalarMatches(at, kHvrcmGain * hvrcmTransition(), "HVRCM residual at the threshold"); + success *= scalarMatches(above - below, kHvrcmGain * kHvrcmOffset, "HVRCM symmetric residual difference"); - success *= (fixture.evaluate() == 0); - - const auto* f = fixture.regca.getResidual().getData(); - success *= scalarMatches(f[index(Vars::IQEXTRA)], - test_case.expected_residual, - "HVRCM residual"); - } + const RealT shifted = residualAt(kHvrcmVoltageLimit, 0.1); + success *= scalarMatches(shifted - at, -0.1, "HVRCM extra-current residual shift"); } - // At the transition point both partial derivatives of the HVRCM row - // are exactly one half. + // At the threshold, the voltage sensitivity is half the Khv gain. { using DepVar = DependencyTracking::Variable; @@ -667,8 +649,8 @@ namespace GridKit success *= fixture.initialize(); auto* y = fixture.regca.y().getData(); - y[index(Vars::VT)] = kHvrcmVoltageLimit - hvrcmTransition(); - y[index(Vars::IQEXTRA)] = hvrcmTransition(); + y[index(Vars::VT)] = kHvrcmVoltageLimit; + y[index(Vars::IQEXTRA)] = 0.0; fixture.regca.y().setDataUpdated(); numberVariables(fixture); @@ -678,8 +660,8 @@ namespace GridKit fixture.regca.getResidual().getData()[index(Vars::IQEXTRA)].getDependencies(); const DepVar::DependencyMap expected{{ - {index(Vars::VT), 0.5}, - {index(Vars::IQEXTRA), -0.5}, + {index(Vars::VT), 0.5 * kHvrcmGain}, + {index(Vars::IQEXTRA), -1.0}, }}; success *= isEqual(dependencies, expected, kTol); } @@ -837,23 +819,26 @@ namespace GridKit // Bus and command values shared by the residual and Jacobian fixtures. static constexpr RealT kStateVr = 0.95; static constexpr RealT kStateVi = 0.25; + static constexpr RealT kStateIp = 0.2; static constexpr RealT kStateIpcmd = 0.9; static constexpr RealT kStateIqcmd = 0.1; - static constexpr RealT kVa1 = 0.9; - static constexpr RealT kJustBelowVa1 = kVa1 - 1.0e-6; + static constexpr RealT kVa1 = 0.9; /// Vhvmax in makeData() and in makeDynamicData(). static constexpr RealT kHvrcmVoltageLimit = 1.2; static constexpr RealT kDynamicHvrcmVoltageLimit = 1.3; - /// The smooth ramp at a zero argument. Used as an HVRCM voltage - /// margin, this is the one margin whose constraint root is the margin - /// itself, so the HVRCM row is exactly zero there and its two partial - /// derivatives are exactly one half. - static RealT hvrcmTransition() + static constexpr RealT kHvrcmGain = 0.7; + + /// Offset whose smooth-ramp tail decays by one significand width. + static constexpr RealT kHvrcmOffset = + static_cast(std::numeric_limits::digits) + * std::numbers::ln2_v / Math::MU; + + static constexpr RealT hvrcmTransition() { - return Math::ramp(static_cast(0.0)); + return std::numbers::ln2_v / Math::MU; } Data makeData() @@ -883,7 +868,7 @@ namespace GridKit data.parameters[Params::VL0] = 0.4; data.parameters[Params::VL1] = 0.9; data.parameters[Params::VA0] = 0.4; - data.parameters[Params::VA1] = 0.9; + data.parameters[Params::VA1] = kVa1; data.parameters[Params::Vhvmax] = kHvrcmVoltageLimit; return data; @@ -918,17 +903,17 @@ namespace GridKit } #endif - /// Writes the answer-key state shared by the residual and Jacobian - /// tests. Values only; in the dependency-tracking fixture, number - /// variables after this call. + /// Write the state shared by the residual and Jacobian tests. In a + /// dependency-tracking fixture, number variables after this call. template - void setAnswerKeyState(PhasorDynamics::Converter::Regca& regca) + void setResidualState(PhasorDynamics::Converter::Regca& regca) { - auto* y = regca.y().getData(); - y[index(Vars::VM)] = 0.65; + const RealT v = 0.65; + auto* y = regca.y().getData(); + y[index(Vars::VM)] = v; y[index(Vars::IQ)] = -0.2; - y[index(Vars::IP)] = 0.85; - y[index(Vars::VT)] = 1.0; + y[index(Vars::IP)] = kStateIp; + y[index(Vars::VT)] = v; y[index(Vars::IR)] = 0.5; y[index(Vars::II)] = 0.18; y[index(Vars::IQEXTRA)] = 0.03; @@ -1003,9 +988,12 @@ namespace GridKit return success; } - bool scalarMatches(ScalarT actual, ScalarT expected, const char* label) const + bool scalarMatches(ScalarT actual, + ScalarT expected, + const char* label, + RealT tolerance = kTol) const { - if (isEqual(actual, expected, kTol)) + if (isEqual(actual, expected, tolerance)) { return true; } @@ -1017,19 +1005,19 @@ namespace GridKit } #ifdef GRIDKIT_ENABLE_ENZYME - /// Move the answer-key state to the HVRCM and active-current limiter + /// Move the residual state to the HVRCM and active-current limiter /// transition points, where the Jacobian has the richest structure. template void setJacobianState( PhasorDynamics::Converter::Regca& regca, RealT current) { - setAnswerKeyState(regca); + setResidualState(regca); auto* y = regca.y().getData(); y[index(Vars::IP)] = current; y[index(Vars::IL)] = 0.0; - y[index(Vars::VT)] = kDynamicHvrcmVoltageLimit - hvrcmTransition(); - y[index(Vars::IQEXTRA)] = hvrcmTransition(); + y[index(Vars::VT)] = kDynamicHvrcmVoltageLimit; + y[index(Vars::IQEXTRA)] = kHvrcmGain * hvrcmTransition(); regca.y().setDataUpdated(); } diff --git a/tests/UnitTests/PhasorDynamics/GovernorGastPtiTests.hpp b/tests/UnitTests/PhasorDynamics/GovernorGastPtiTests.hpp index ff691d28d..960bdf41b 100644 --- a/tests/UnitTests/PhasorDynamics/GovernorGastPtiTests.hpp +++ b/tests/UnitTests/PhasorDynamics/GovernorGastPtiTests.hpp @@ -19,6 +19,7 @@ #include #include #include +#include namespace GridKit { @@ -846,12 +847,7 @@ namespace GridKit const auto data = makeResidualData(); - Fixture enzyme(data); - enzyme.attachAllInputs(); - success *= enzyme.initialize(initial_pmech); - - const auto compare = [&](Fixture& enzyme_fixture, - const Data& case_data, + const auto compare = [&](const Data& case_data, RealT pmech, const char* context, std::initializer_list overrides) @@ -859,7 +855,7 @@ namespace GridKit const auto dependency_jacobian = dependencyTrackingJacobian(case_data, pmech, success, overrides); const auto enzyme_jacobian = - enzymeJacobian(enzyme_fixture, success, overrides, context); + enzymeJacobian(case_data, pmech, success, overrides, context); success *= jacobianMatches(enzyme_jacobian, dependency_jacobian, @@ -867,28 +863,23 @@ namespace GridKit kTol); }; - compare(enzyme, - data, + compare(data, initial_pmech, "load-limited Enzyme versus dependency tracking", {}); - compare(enzyme, - data, + compare(data, initial_pmech, "temperature-limited Enzyme versus dependency tracking", {{Internal::VLOAD, 1.5}, {Internal::VTEMP, 0.3}}); - compare(enzyme, - data, + compare(data, initial_pmech, "equal-selector Enzyme versus dependency tracking", {{Internal::VLOAD, 0.9}, {Internal::VTEMP, 0.9}}); - compare(enzyme, - data, + compare(data, initial_pmech, "blocked-response Enzyme versus dependency tracking", {{Internal::XVALVE, 1.6}, {Internal::VLV, 1.85}}); - compare(enzyme, - data, + compare(data, initial_pmech, "restoring-response Enzyme versus dependency tracking", {{Internal::XVALVE, 1.6}, {Internal::VLV, 1.35}}); @@ -903,8 +894,7 @@ namespace GridKit const RealT adjusted_boundary = static_cast( adjusted.gastpti.y().getData()[index(Internal::XVALVE)]); - compare(adjusted, - data, + compare(data, over_rated_pmech, "adjusted-boundary Enzyme versus dependency tracking", {{Internal::XVALVE, adjusted_boundary}, @@ -922,8 +912,7 @@ namespace GridKit Fixture collapsed(collapsed_data); collapsed.attachAllInputs(); success *= collapsed.initialize(initial_pmech); - compare(collapsed, - collapsed_data, + compare(collapsed_data, initial_pmech, "collapsed Enzyme versus dependency tracking", {}); @@ -932,8 +921,7 @@ namespace GridKit collapsed.input(index(External::PREF)) = ZERO; collapsed.seedPmech(over_rated_pmech); success *= (collapsed.gastpti.initialize() == 0); - compare(collapsed, - collapsed_data, + compare(collapsed_data, over_rated_pmech, "reinitialized Enzyme versus dependency tracking", {}); @@ -1710,11 +1698,15 @@ namespace GridKit #ifdef GRIDKIT_ENABLE_ENZYME std::vector enzymeJacobian( - Fixture& fixture, + const Data& data, + RealT pmech, TestStatus& success, std::initializer_list overrides, const char* context) const { + Fixture fixture(data); + fixture.attachAllInputs(); + success *= fixture.initialize(pmech); setAnswerKeyInputs(fixture); setAnswerKeyState(fixture.gastpti); setState(fixture.gastpti, overrides); @@ -1726,32 +1718,13 @@ namespace GridKit return {}; } - auto* coo = fixture.gastpti.getCooJacobian(); - if (coo == nullptr || fixture.gastpti.nnz() != coo->getNnz()) + if (fixture.gastpti.constructCsr() != 0) { - std::cout << "GASTPTI COO structure changed for " << context << '\n'; + std::cout << "GASTPTI CSR construction failed for " << context << '\n'; success = false; return {}; } - - std::vector rows( - static_cast(fixture.gastpti.size())); - const auto* row_indices = coo->getRowData(); - const auto* column_indices = coo->getColData(); - const auto* values = coo->getValues(); - - for (IdxT entry = 0; entry < coo->getNnz(); ++entry) - { - const auto row = static_cast(row_indices[entry]); - if (row >= rows.size()) - { - std::cout << "GASTPTI COO row is invalid for " << context << '\n'; - success = false; - return {}; - } - rows[row][static_cast(column_indices[entry])] += values[entry]; - } - return rows; + return MapFromCsr(fixture.gastpti.getCsrJacobian()); } #endif