From 9f3952461badb6856d7fe139e1b7fe6fbfce2215 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Fri, 26 Jun 2026 08:17:40 -0500 Subject: [PATCH 1/2] Add REGCA Implementation --- CHANGELOG.md | 1 + GridKit/CommonMath.hpp | 30 +- GridKit/Model/PhasorDynamics/CMakeLists.txt | 1 + .../Model/PhasorDynamics/ComponentLibrary.hpp | 1 + .../PhasorDynamics/Converter/CMakeLists.txt | 6 + .../PhasorDynamics/Converter/REECA/README.md | 8 +- .../Converter/REGCA/CMakeLists.txt | 54 + .../PhasorDynamics/Converter/REGCA/README.md | 417 ++++--- .../PhasorDynamics/Converter/REGCA/Regca.cpp | 27 + .../PhasorDynamics/Converter/REGCA/Regca.hpp | 200 ++++ .../Converter/REGCA/RegcaData.hpp | 73 ++ .../REGCA/RegcaDependencyTracking.cpp | 27 + .../Converter/REGCA/RegcaEnzyme.cpp | 124 ++ .../Converter/REGCA/RegcaImpl.hpp | 567 +++++++++ GridKit/Model/PhasorDynamics/INPUT_FORMAT.md | 5 +- .../Model/PhasorDynamics/SystemModelData.hpp | 17 + .../SystemModelDataJSONParser.hpp | 6 + .../Model/PhasorDynamics/SystemModelImpl.hpp | 57 + .../Math/SmoothnessIndicatorTests.hpp | 50 +- .../Math/runSmoothnessIndicatorTests.cpp | 3 +- tests/UnitTests/PhasorDynamics/CMakeLists.txt | 10 + .../PhasorDynamics/ConverterRegcaTests.hpp | 1026 +++++++++++++++++ .../SystemSingleComponentTests.hpp | 51 + .../PhasorDynamics/runConverterRegcaTests.cpp | 28 + .../runSystemSingleComponentTests.cpp | 1 + 25 files changed, 2557 insertions(+), 233 deletions(-) create mode 100644 GridKit/Model/PhasorDynamics/Converter/CMakeLists.txt create mode 100644 GridKit/Model/PhasorDynamics/Converter/REGCA/CMakeLists.txt create mode 100644 GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.cpp create mode 100644 GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp create mode 100644 GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp create mode 100644 GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaDependencyTracking.cpp create mode 100644 GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp create mode 100644 GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp create mode 100644 tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp create mode 100644 tests/UnitTests/PhasorDynamics/runConverterRegcaTests.cpp diff --git a/CHANGELOG.md b/CHANGELOG.md index ba0f1dde6..07391f1e5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -67,6 +67,7 @@ - Added `ConstantSignalSource` component as first use case for `BusToSignalAdapter`. - Added IDA option to suppress algebraic variables in local error tests. - Removed `COO_Matrix` class. +- Added `REGCA` converter model implementation for PhasorDynamics. ## v0.1 diff --git a/GridKit/CommonMath.hpp b/GridKit/CommonMath.hpp index e4a7ec87f..8f9f8c5cf 100644 --- a/GridKit/CommonMath.hpp +++ b/GridKit/CommonMath.hpp @@ -162,7 +162,31 @@ namespace GridKit } /** - * @brief Smooth two-sided deadband function + * @brief Smooth Type 1 no-offset two-sided deadband function + * + * Smooth approximation to a deadband that returns zero inside the band and + * passes the input through unchanged outside the band. + * + * @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 no-offset deadbanded value + */ + template + __attribute__((always_inline)) inline ScalarT deadband1( + const ScalarT x, + const RealT lower, + const RealT upper) + { + assert(lower <= upper); + return x * (sigmoid(lower - x) + sigmoid(x - upper)); + } + + /** + * @brief Smooth Type 2 offset two-sided deadband function * * Smooth approximation to x - min(max(x, lower), upper), composed from the * smooth ramp function. @@ -173,10 +197,10 @@ namespace GridKit * @param[in] x - Input signal * @param[in] lower - Lower breakpoint * @param[in] upper - Upper breakpoint - * @return Smooth deadbanded value + * @return Smooth offset deadbanded value */ template - __attribute__((always_inline)) inline ScalarT deadband( + __attribute__((always_inline)) inline ScalarT deadband2( const ScalarT x, const RealT lower, const RealT upper) diff --git a/GridKit/Model/PhasorDynamics/CMakeLists.txt b/GridKit/Model/PhasorDynamics/CMakeLists.txt index baf76ab90..91bb8d04b 100644 --- a/GridKit/Model/PhasorDynamics/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/CMakeLists.txt @@ -33,6 +33,7 @@ add_subdirectory(Branch) add_subdirectory(Bus) add_subdirectory(BusFault) add_subdirectory(BusToSignalAdapter) +add_subdirectory(Converter) add_subdirectory(Exciter) add_subdirectory(Governor) add_subdirectory(Load) diff --git a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp index 522ff5fd3..6ac2375f3 100644 --- a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp +++ b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include diff --git a/GridKit/Model/PhasorDynamics/Converter/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Converter/CMakeLists.txt new file mode 100644 index 000000000..cafb2cc36 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/CMakeLists.txt @@ -0,0 +1,6 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +add_subdirectory(REGCA) diff --git a/GridKit/Model/PhasorDynamics/Converter/REECA/README.md b/GridKit/Model/PhasorDynamics/Converter/REECA/README.md index 68c8ab96f..d2a7c743f 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REECA/README.md +++ b/GridKit/Model/PhasorDynamics/Converter/REECA/README.md @@ -159,7 +159,7 @@ Symbol | Units | Description | $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` +$V_{\mathrm{err}}$ | [p.u.] | Deadbanded voltage error | Defined by CommonMath `deadband2` $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}}$ @@ -260,7 +260,7 @@ The algebraic targets use CommonMath helper notation where applicable: 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 &= -V_\mathrm{err} + \text{deadband2}(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}) @@ -288,7 +288,7 @@ The algebraic targets use CommonMath helper notation where applicable: 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). +CommonMath defines the helper targets and smooth approximations for [min, max, clamp, deadband2, and outside](../../../../CommonMath.md#derived-functions). ## Initialization @@ -322,7 +322,7 @@ Then evaluate the upstream algebraic chain: \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}}) \\ + V_{\mathrm{err},0} &= \text{deadband2}(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} \\ diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Converter/REGCA/CMakeLists.txt new file mode 100644 index 000000000..3ffe1fef7 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/CMakeLists.txt @@ -0,0 +1,54 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +set(_install_headers Regca.hpp RegcaData.hpp) + +if(GRIDKIT_ENABLE_ENZYME) + gridkit_add_library( + phasor_dynamics_converter_regca + SOURCES RegcaEnzyme.cpp + HEADERS ${_install_headers} + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal + PRIVATE + ClangEnzymeFlags + COMPILE_OPTIONS + PRIVATE + -mllvm + -enzyme-auto-sparsity=1 + -fno-math-errno) +else() + gridkit_add_library( + phasor_dynamics_converter_regca + SOURCES Regca.cpp + HEADERS ${_install_headers} + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal) +endif() + +gridkit_add_library( + phasor_dynamics_converter_regca_dependency_tracking + SOURCES RegcaDependencyTracking.cpp + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal_dependency_tracking) + +target_link_libraries( + phasor_dynamics_components + INTERFACE GridKit::phasor_dynamics_converter_regca) +target_link_libraries( + phasor_dynamics_components_dependency_tracking + INTERFACE GridKit::phasor_dynamics_converter_regca_dependency_tracking) diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md b/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md index 45b63f3a7..951f17dbc 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md @@ -4,10 +4,13 @@ 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}}$. +## 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$. +- HVRCM is represented by algebraic current $I_q^\mathrm{extra}$. +- PowerWorld fields `Qmin`, `Khv`, and `Xe` are not parameters of this GridKit implementation. ## Block Diagram @@ -19,55 +22,66 @@ Figure 1: Generator/Converter REGCA model. Figure courtesy of [PowerWorld](https ## 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` +Symbol | Units | JSON | Description | Typical Value | Note +---------------------------------|----------|----------|-------------------------------------------------------|---------------|------ +$P_0$ | [p.u.] | `P0` | Initial active power injection | 1.0 | System base; required initialization source +$Q_0$ | [p.u.] | `Q0` | Initial reactive power injection | 0.0 | System base; required initialization source +$S^\mathrm{base}$ | [MVA] | `mva` | REGCA component power base | 100.0 | +$T_{\mathrm{g}}$ | [sec] | `Tg` | Converter current-control lag time constant | 0.02 | Block name: `Tg` +$T_M$ | [sec] | `TM` | Terminal voltage sensor time constant | 0.02 | Block name: `Tfltr` +$R_q^{\max}$ | [p.u./s] | `Rqmax` | Reactive-current recovery positive rate limit | 999.0 | Block name: `Iqrmax`; GridKit requires $R_q^{\max} > 0$ +$R_q^{\min}$ | [p.u./s] | `Rqmin` | Reactive-current recovery negative rate limit | -999.0 | Block name: `Iqrmin`; GridKit requires $R_q^{\min} < 0$ +$R_p^{\max}$ | [p.u./s] | `Rpmax` | Active-current magnitude recovery rate limit | 999.0 | Block name: `rrpwr` +$s_L$ | [binary] | `sL` | LVPL switch | 1 | Block name: `LPVLSW` +$I_{L1}$ | [p.u.] | `IL1` | LVPL upper-current ceiling | 1.1 | Block name: `LVPL1` +$V_{L0}$ | [p.u.] | `VL0` | LVPL zero-crossing voltage | 0.4 | Block name: `zerox` +$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` ### Parameter Validation -Implementations should reject or report invalid parameter sets: +Invalid REGCA parameter sets are rejected by the following checks. Let $\epsilon_T=10^{-3}$. ```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 + T &\leftarrow \max\!\left(T, \epsilon_T\right) + \quad T\in\{T_{\mathrm{g}},T_M\} \\ + S^\mathrm{base} + &> 0 \\ + R_p^{\max} + &> 0 \\ + R_q^{\min} + &< 0 < R_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} +\begin{aligned} + s_L^\mathrm{off} + &= 1 - s_L \\ + k_{\mathrm{base}} + &= \dfrac{S^\mathrm{sys}}{S^\mathrm{base}} \\ + M_p + &= 100 R_p^{\max} +\end{aligned} ``` -$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. +$M_p$ is a finite surrogate for inactive active-current rate limits, +not a physical REGCA parameter. ## Model Variables @@ -78,20 +92,22 @@ moderate enough to keep the smooth clamp well conditioned. 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 +$I_q$ | [p.u.] | Reactive-current state | State 1 in Fig. 1 before the `-1` block; component base +$I_p$ | [p.u.] | Active-current state | State 2 in Fig. 1; component 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$ +Symbol | Units | Description | Note +---------------------------|----------|-------------------------------------------------------------|------ +$V_T$ | [p.u.] | Terminal voltage magnitude | +$I_{\mathrm{r}}$ | [p.u.] | Branch-current real component | System base +$I_{\mathrm{i}}$ | [p.u.] | Branch-current imaginary component | System base +$I_q^\mathrm{extra}$ | [p.u.] | Extra inductive current from high-voltage reactive current management | Component base +$I_L$ | [p.u.] | LVPL upper-limit current curve | Component base; function of $V_M$ +$\ell_p$ | [p.u./s] | Smooth active-current lower rate bound | Component base; equivalent to diagram `Rdown` +$u_p$ | [p.u./s] | Smooth active-current upper rate bound | Component base; effective `Rup` +$P^\mathrm{br}$ | [p.u.] | Branch active power | System base +$Q^\mathrm{br}$ | [p.u.] | Branch reactive power | System base ### External Variables @@ -100,210 +116,191 @@ 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 +Symbol | Units | Type | Description | Note +--------------------------------|--------|---------|------------------------------------------------------------------|------ +$V_{\mathrm{r}}$ | [p.u.] | Known | Terminal voltage, real component | Bus input +$V_{\mathrm{i}}$ | [p.u.] | Known | Terminal voltage, imaginary component | Bus input +$I_p^\mathrm{cmd}$ | [p.u.] | Unknown | Active-current command in the terminal-voltage reference frame | Optional signal port `ipcmd`; system base +$I_q^\mathrm{cmd}$ | [p.u.] | Unknown | Reactive-current command in the terminal-voltage reference frame | Optional signal port `iqcmd`; system base -## Model Equations +REGCA initializes the current-command ports from the required `P0`/`Q0` +power-flow injection. If no controller is connected, the resolved initialization +commands are held constant during residual evaluation. -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} -``` +## Model Equations ### 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 +The state equations use CommonMath helper notation. The $I_q$ +limiter branch is selected by the initialized reactive-current command. ```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}}) + 0 &= -\dot V_M + \dfrac{1}{T_M}\left(V_T - V_M\right) \\ + 0 &= -\dot I_q + + \dfrac{1}{T_{\mathrm{g}}}\left(k_{\mathrm{base}}I_q^\mathrm{cmd} - I_q\right) + + \dfrac{1}{T_{\mathrm{g}}} + \begin{cases} + -\rho\left(k_{\mathrm{base}}I_q^\mathrm{cmd} - I_q - T_{\mathrm{g}}R_q^{\max}\right) + & I_{q,0}^\mathrm{cmd} > 0 \\ + \rho\left(T_{\mathrm{g}}R_q^{\min} - \left(k_{\mathrm{base}}I_q^\mathrm{cmd} - I_q\right)\right) + & I_{q,0}^\mathrm{cmd} \le 0 + \end{cases} \\ + 0 &= -\dot I_p + + \ell_p + + \dfrac{1}{T_{\mathrm{g}}} + \rho\left(k_{\mathrm{base}}I_p^\mathrm{cmd} - I_p - T_{\mathrm{g}}\ell_p\right) + - \dfrac{1}{T_{\mathrm{g}}} + \rho\left(k_{\mathrm{base}}I_p^\mathrm{cmd} - I_p - T_{\mathrm{g}}u_p\right) \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}}$. +Here $\rho$ is GridKit's [ramp](../../../../CommonMath.md#primitives). ### 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 &= -V_T I_{\mathrm{r}} + + \dfrac{1}{k_{\mathrm{base}}} + \left[ + V_{\mathrm{i}}\left(I_q - I_q^\mathrm{extra}\right) + + V_{\mathrm{r}}I_p\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1) + \right] \\ + 0 &= -V_T I_{\mathrm{i}} + + \dfrac{1}{k_{\mathrm{base}}} + \left[ + -V_{\mathrm{r}}\left(I_q - I_q^\mathrm{extra}\right) + + V_{\mathrm{i}}I_p\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1) + \right] \\ + 0 &= -I_q^\mathrm{extra} + + \rho(V_T - V_{\mathrm{hv}}^{\max}) \\ 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} + 0 &= -\ell_p + - R_p^{\max} + - \left(M_p - R_p^{\max}\right)\sigma(I_p) \\ + 0 &= -u_p + + M_p\left(1-\sigma(I_p)\right) + + R_p^{\max}\sigma(I_p) + \left(s_L^\mathrm{off} + s_L\sigma(I_L - I_p)\right) \\ + 0 &= -P^\mathrm{br} + + V_{\mathrm{r}} I_{\mathrm{r}} + V_{\mathrm{i}} I_{\mathrm{i}} \\ + 0 &= -Q^\mathrm{br} + + V_{\mathrm{i}} I_{\mathrm{r}} - V_{\mathrm{r}} I_{\mathrm{i}} \end{aligned} ``` -The $V_T$ residual is kept in squared form for smoothness at the origin. +Here $\text{linseg}$ is GridKit's +[linear segment](../../../../CommonMath.md#derived-functions), $\rho$ is the +[ramp](../../../../CommonMath.md#primitives), and $\sigma$ is the +[step function](../../../../CommonMath.md#primitives). ## Network Interface -The bus receives system-base current injections converted from converter-base -REGCA currents: +The bus receives the REGCA branch-current variables directly: ```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}}} + I_{\mathrm{r}}^\mathrm{inj} &:= I_{\mathrm{r}} \\ + I_{\mathrm{i}}^\mathrm{inj} &:= I_{\mathrm{i}} \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: +### External Priors ```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} + V_{\mathrm{r},0}, V_{\mathrm{i},0} + &\leftarrow \text{terminal-bus voltage} \\ + P_0, Q_0 + &\leftarrow \text{power-flow injection on system base} \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. +### Internal Initialization -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. +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: -All internal derivatives initialize to zero. +```math +\begin{aligned} + V_{T,0} + &= \sqrt{V_{\mathrm{r},0}^2 + V_{\mathrm{i},0}^2} \\ + V_{M,0} + &= V_{T,0} \\ + A_0^\mathrm{LVACM} + &= \text{linseg}\left(V_{T,0};\ V_{A0},\ V_{A1},\ 1\right) \\ + I_{p,0}^\mathrm{cmd} + &= + \begin{cases} + \dfrac{P_0}{V_{T,0}A_0^\mathrm{LVACM}} & P_0 \ne 0 \\ + 0 & P_0 = 0 + \end{cases} \\ + I_{q,0}^\mathrm{cmd} + &= \dfrac{Q_0}{V_{T,0}} \\ + I_{p,0} + &= k_{\mathrm{base}}I_{p,0}^\mathrm{cmd} \\ + I_{q,0} + &= k_{\mathrm{base}}I_{q,0}^\mathrm{cmd} \\ + I_{q,0}^\mathrm{extra} + &= 0 \\ + I_{L,0} + &= \text{linseg}\left(V_{T,0};\ V_{L0},\ V_{L1},\ I_{L1}\right) \\ + \ell_{p,0} + &= -R_p^{\max} + - \left(M_p - R_p^{\max}\right)\sigma(I_{p,0}) \\ + u_{p,0} + &= M_p\left(1-\sigma(I_{p,0})\right) + + R_p^{\max}\sigma(I_{p,0}) + \left(s_L^\mathrm{off} + s_L\sigma(I_{L,0} - I_{p,0})\right) \\ + I_{\mathrm{r},0} + &= \dfrac{1}{k_{\mathrm{base}}V_{T,0}} + \left[ + V_{\mathrm{i},0}\left(I_{q,0} - I_{q,0}^\mathrm{extra}\right) + + V_{\mathrm{r},0}I_{p,0}A_0^\mathrm{LVACM} + \right] \\ + I_{\mathrm{i},0} + &= \dfrac{1}{k_{\mathrm{base}}V_{T,0}} + \left[ + -V_{\mathrm{r},0}\left(I_{q,0} - I_{q,0}^\mathrm{extra}\right) + + V_{\mathrm{i},0}I_{p,0}A_0^\mathrm{LVACM} + \right] \\ + P_0^\mathrm{br} + &= V_{\mathrm{r},0}I_{\mathrm{r},0} + V_{\mathrm{i},0}I_{\mathrm{i},0} \\ + Q_0^\mathrm{br} + &= V_{\mathrm{i},0}I_{\mathrm{r},0} - V_{\mathrm{r},0}I_{\mathrm{i},0} +\end{aligned} +``` -## Model Outputs +Here $\text{linseg}$ is GridKit's +[linear segment](../../../../CommonMath.md#derived-functions), and $\sigma$ is +the [step function](../../../../CommonMath.md#primitives). + +### External Solved -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}} + I_p^\mathrm{cmd} + &\leftarrow I_{p,0}^\mathrm{cmd} \\ + I_q^\mathrm{cmd} + &\leftarrow I_{q,0}^\mathrm{cmd} \end{aligned} ``` -Power outputs are positive leaving the converter and entering the bus. + +REGCA writes the resolved current commands to attached `ipcmd` and `iqcmd` +ports. If no controller is attached, those values are used as constant command +inputs. + +## Monitorable Outputs + +Output | Units | Description | Note +-------|--------|-----------------------------|------ +`ir` | [p.u.] | Real current injection | System base; exported through `ibranchr` when assigned +`ii` | [p.u.] | Imaginary current injection | System base; exported through `ibranchi` when assigned +`p` | [p.u.] | Active-power output | System base; exported through `pbranch` when assigned +`q` | [p.u.] | Reactive-power output | System base; exported through `qbranch` when assigned diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.cpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.cpp new file mode 100644 index 000000000..0f6a2aca6 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.cpp @@ -0,0 +1,27 @@ +/** + * @file Regca.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Non-Enzyme instantiation for the REGCA converter model. + */ + +#include "RegcaImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + template + int Regca::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Regca..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Regca; + template class Regca; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp new file mode 100644 index 000000000..dee03de77 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp @@ -0,0 +1,200 @@ +/** + * @file Regca.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Declaration of the REGCA phasor-dynamics converter model. + */ + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + template + class BusBase; + + template + class SignalNode; + } // namespace PhasorDynamics +} // namespace GridKit + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + /// Internal variables of a `Regca` + enum class RegcaInternalVariables : size_t + { + VM, ///< Filtered terminal voltage + IQ, ///< Reactive-current state + IP, ///< Active-current state + VT, ///< Terminal voltage magnitude + IR, ///< Real injected current on system base + II, ///< Imaginary injected current on system base + IQEXTRA, ///< HVRCM extra reactive current + IL, ///< LVPL upper-limit current curve + LP, ///< Active-current lower rate bound + UP, ///< Active-current upper rate bound + PBR, ///< Active-power output on system base + QBR, ///< Reactive-power output on system base + MAXIMUM, + }; + + /// External variables of a `Regca` + enum class RegcaExternalVariables : size_t + { + IPCMD, ///< Active-current command signal on system base + IQCMD, ///< Reactive-current command signal on system base + MAXIMUM, + }; + + template + class Regca : public Component + { + using Component::gridkit_component_id_; + using Component::alpha_; + using Component::abs_tol_; + using Component::f_; + using Component::h_; + using Component::J_cols_buffer_; + using Component::J_rows_buffer_; + using Component::J_vals_buffer_; + using Component::nnz_; + using Component::residual_indices_; + using Component::size_; + using Component::tag_; + using Component::time_; + using Component::va_system_base_; + using Component::variable_indices_; + using Component::wb_; + using Component::y_; + using Component::yp_; + + public: + using ScalarT = scalar_type; + using IdxT = index_type; + using RealT = typename Component::RealT; + using BusT = BusBase; + using SignalT = SignalNode; + using ModelDataT = RegcaData; + using MonitorT = Model::VariableMonitor; + + Regca(BusT* bus); + Regca(BusT* bus, const ModelDataT& data); + ~Regca(); + + int setGridKitComponentID(IdxT) 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& + { + return signals_; + } + + const Model::VariableMonitorBase* getMonitor() const override; + + __attribute__((always_inline)) inline int evaluateInternalResidual( + const ScalarT*, const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); + __attribute__((always_inline)) inline int evaluateBusResidual( + const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); + + private: + void initializeParameters(const ModelDataT& data); + void initializeMonitor(); + void setDerivedParameters(); + + ScalarT toComponentBase(ScalarT value) const + { + return value * va_system_base_ / va_converter_base_; + } + + ScalarT toSystemBase(ScalarT value) const + { + return value / toComponentBase(static_cast(ONE)); + } + + ScalarT lpTarget(ScalarT ip) const; + ScalarT upTarget(ScalarT ip, ScalarT il) const; + + ScalarT& Vr() + { + return bus_->Vr(); + } + + ScalarT& Vi() + { + return bus_->Vi(); + } + + ScalarT& Ir() + { + return bus_->Ir(); + } + + ScalarT& Ii() + { + return bus_->Ii(); + } + + static constexpr RealT TIME_CONSTANT_MINIMUM = static_cast(1.0e-3); + + BusT* bus_{nullptr}; + + RealT P0_{0}; + RealT Q0_{0}; + RealT mva_base_{0}; + RealT Tg_{0}; + RealT TM_{0}; + RealT Rqmax_{0}; + RealT Rqmin_{0}; + RealT Rpmax_{0}; + bool sL_{false}; + RealT IL1_{0}; + RealT VL0_{0}; + RealT VL1_{0}; + RealT VA0_{0}; + RealT VA1_{0}; + RealT Vhvmax_{0}; + IdxT bus_id_{0}; + + IdxT parameter_error_count_{0}; + RealT Mp_{0}; + RealT va_converter_base_{0}; + RealT use_lvpl_{0}; + RealT bypass_lvpl_{1}; + RealT iq_use_upper_{0}; + RealT iq_use_lower_{1}; + + ScalarT ipcmd_set_{0}; + ScalarT iqcmd_set_{0}; + + ComponentSignals signals_; + std::unique_ptr monitor_; + + std::vector ws_; + std::vector ws_indices_; + }; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp new file mode 100644 index 000000000..231343604 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp @@ -0,0 +1,73 @@ +/** + * @file RegcaData.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Modeling data for the REGCA converter model. + */ + +#pragma once + +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + /// Parameter keys for the REGCA converter model. + enum class RegcaParameters + { + P0, ///< Initial active power on system base + Q0, ///< Initial reactive power on system base + mva, ///< MVA base of the REGCA model + Tg, ///< Converter current-control lag time constant + TM, ///< Terminal voltage sensor time constant + Rqmax, ///< Reactive-current recovery positive rate limit + Rqmin, ///< Reactive-current recovery negative rate limit + Rpmax, ///< Active-current magnitude recovery rate limit + sL, ///< LVPL switch + IL1, ///< LVPL upper-current ceiling + VL0, ///< LVPL zero-crossing voltage + VL1, ///< LVPL upper breakpoint voltage + VA0, ///< LVACM lower breakpoint voltage + VA1, ///< LVACM upper breakpoint voltage + Vhvmax ///< Terminal-voltage ceiling for HV reactive management + }; + + /// Ports for the REGCA converter model. + enum class RegcaPorts + { + bus, ///< Terminal bus ID + ipcmd, ///< Optional active-current command signal ID + iqcmd, ///< Optional reactive-current command signal ID + ibranchr, ///< Optional real current measurement output signal ID + ibranchi, ///< Optional imaginary current measurement output signal ID + pbranch, ///< Optional active-power measurement output signal ID + qbranch ///< Optional reactive-power measurement output signal ID + }; + + /// Variables available through the monitor interface. + enum class RegcaMonitorableVariables + { + ir, ///< Real current injection on system base + ii, ///< Imaginary current injection on system base + p, ///< Active power injection on system base + q ///< Reactive power injection on system base + }; + + template + struct RegcaData : public ComponentData + { + RegcaData() = default; + + using Parameters = RegcaParameters; + using Ports = RegcaPorts; + using MonitorableVariables = RegcaMonitorableVariables; + }; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaDependencyTracking.cpp new file mode 100644 index 000000000..bf470cc26 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaDependencyTracking.cpp @@ -0,0 +1,27 @@ +/** + * @file RegcaDependencyTracking.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Dependency-tracking instantiations for the REGCA converter model. + */ + +#include "RegcaImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + template + int Regca::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Regca..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Regca; + template class Regca; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp new file mode 100644 index 000000000..0fc2d1b41 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp @@ -0,0 +1,124 @@ +/** + * @file RegcaEnzyme.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Enzyme sparse Jacobian for the REGCA converter model. + */ + +#include + +#include "RegcaImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + template + int Regca::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Regca..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + + if (J_rows_buffer_ == nullptr) + { + // Reserve space for the dense blocks. Enzyme keeps only structural + // nonzeros for each differentiated block. + auto size = static_cast(size_); + auto bus_size = static_cast(bus_->size()); + auto signal_size = static_cast(ws_.size()); + auto buffer_size = 2 * size * size + 2 * size * bus_size + size * signal_size; + J_rows_buffer_ = new IdxT[buffer_size]; + J_cols_buffer_ = new IdxT[buffer_size]; + J_vals_buffer_ = new RealT[buffer_size]; + } + + using RegcaT = GridKit::PhasorDynamics::Converter::Regca; + using Fn = GridKit::Enzyme::Sparse::MemberFunctions; + + nnz_ = 0; + + GridKit::Enzyme::Sparse::DfDy::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + GridKit::Enzyme::Sparse::DfDyp::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + GridKit::Enzyme::Sparse::DfDwb::eval(this, + f_.size(), + static_cast(bus_->size()), + (this->getResidualIndices()).data(), + (bus_->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + GridKit::Enzyme::Sparse::DfDws::eval(this, + f_.size(), + ws_.size(), + (this->getResidualIndices()).data(), + ws_indices_.data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + GridKit::Enzyme::Sparse::DhDy::eval(this, + static_cast(bus_->size()), + y_.size(), + (bus_->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + this->constructCoo(); + + return 0; + } + + template class Regca; + template class Regca; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp new file mode 100644 index 000000000..ab3fff561 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp @@ -0,0 +1,567 @@ +/** + * @file RegcaImpl.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Definition of the REGCA phasor-dynamics converter model. + */ + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + using Log = ::GridKit::Utilities::Logger; + + template + Regca::Regca(BusT* bus) + : bus_(bus) + { + size_ = static_cast(RegcaInternalVariables::MAXIMUM); + } + + template + Regca::Regca(BusT* bus, const ModelDataT& data) + : bus_(bus), + monitor_(std::make_unique(data)) + { + initializeParameters(data); + initializeMonitor(); + size_ = static_cast(RegcaInternalVariables::MAXIMUM); + } + + template + Regca::~Regca() + { + } + + template + void Regca::setDerivedParameters() + { + Tg_ = std::max(Tg_, TIME_CONSTANT_MINIMUM); + TM_ = std::max(TM_, TIME_CONSTANT_MINIMUM); + Mp_ = static_cast(100.0) * Rpmax_; + use_lvpl_ = ZERO; + bypass_lvpl_ = ONE; + if (sL_) + { + use_lvpl_ = ONE; + bypass_lvpl_ = ZERO; + } + iq_use_upper_ = ZERO; + iq_use_lower_ = ONE; + va_converter_base_ = mva_base_ * static_cast(1.0e6); + } + + template + scalar_type Regca::lpTarget(scalar_type ip) const + { + return -Rpmax_ - (Mp_ - Rpmax_) * Math::sigmoid(ip); + } + + template + scalar_type Regca::upTarget( + scalar_type ip, + scalar_type il) const + { + const ScalarT sigma_ip = Math::sigmoid(ip); + return Mp_ * (ONE - sigma_ip) + + Rpmax_ * sigma_ip * (bypass_lvpl_ + use_lvpl_ * Math::sigmoid(il - ip)); + } + + template + void Regca::initializeParameters(const ModelDataT& data) + { + using Params = typename ModelDataT::Parameters; + using Ports = typename ModelDataT::Ports; + + parameter_error_count_ = 0; + + auto load_required_real = [&](auto key, RealT& target, const char* name) + { + if (!data.parameters.contains(key)) + { + Log::error() << "Regca: missing required parameter '" << name << "'\n"; + ++parameter_error_count_; + return; + } + + const auto& value = data.parameters.at(key); + if (const auto* real_value = std::get_if(&value)) + { + target = *real_value; + } + else if (const auto* index_value = std::get_if(&value)) + { + target = static_cast(*index_value); + } + else + { + Log::error() << "Regca: parameter '" << name << "' must be numeric\n"; + ++parameter_error_count_; + } + }; + + auto load_required_switch = [&](auto key, bool& target, const char* name) + { + if (!data.parameters.contains(key)) + { + Log::error() << "Regca: missing required parameter '" << name << "'\n"; + ++parameter_error_count_; + return; + } + + const auto& value = data.parameters.at(key); + if (const auto* bool_value = std::get_if(&value)) + { + target = *bool_value; + } + else if (const auto* index_value = std::get_if(&value); + index_value && (*index_value == 0 || *index_value == 1)) + { + target = (*index_value == 1); + } + else + { + Log::error() << "Regca: parameter '" << name << "' must be bool or 0/1\n"; + ++parameter_error_count_; + } + }; + + load_required_real(Params::P0, P0_, "P0"); + load_required_real(Params::Q0, Q0_, "Q0"); + load_required_real(Params::mva, mva_base_, "mva"); + load_required_real(Params::Tg, Tg_, "Tg"); + load_required_real(Params::TM, TM_, "TM"); + load_required_real(Params::Rqmax, Rqmax_, "Rqmax"); + load_required_real(Params::Rqmin, Rqmin_, "Rqmin"); + load_required_real(Params::Rpmax, Rpmax_, "Rpmax"); + load_required_switch(Params::sL, sL_, "sL"); + load_required_real(Params::IL1, IL1_, "IL1"); + load_required_real(Params::VL0, VL0_, "VL0"); + load_required_real(Params::VL1, VL1_, "VL1"); + load_required_real(Params::VA0, VA0_, "VA0"); + load_required_real(Params::VA1, VA1_, "VA1"); + load_required_real(Params::Vhvmax, Vhvmax_, "Vhvmax"); + + if (data.ports.contains(Ports::bus)) + { + bus_id_ = data.ports.at(Ports::bus); + } + + setDerivedParameters(); + } + + template + const Model::VariableMonitorBase* Regca::getMonitor() const + { + return monitor_.get(); + } + + template + void Regca::initializeMonitor() + { + using Variable = typename ModelDataT::MonitorableVariables; + auto index = [](RegcaInternalVariables variable) + { + return static_cast(variable); + }; + + monitor_->set(Variable::ir, [this, index] + { return y_[index(RegcaInternalVariables::IR)]; }); + monitor_->set(Variable::ii, [this, index] + { return y_[index(RegcaInternalVariables::II)]; }); + monitor_->set(Variable::p, [this, index] + { return y_[index(RegcaInternalVariables::PBR)]; }); + monitor_->set(Variable::q, [this, index] + { return y_[index(RegcaInternalVariables::QBR)]; }); + } + + template + int Regca::setGridKitComponentID(IdxT component_id) + { + gridkit_component_id_ = component_id; + return 0; + } + + template + int Regca::allocate() + { + size_ = static_cast(RegcaInternalVariables::MAXIMUM); + auto size = static_cast(size_); + + f_.assign(size, ScalarT{0}); + y_.assign(size, ScalarT{0}); + yp_.assign(size, ScalarT{0}); + tag_.assign(size, false); + abs_tol_.assign(size, ScalarT{0}); + variable_indices_.resize(size); + residual_indices_.resize(size); + + wb_.assign(2, ScalarT{0}); + h_.assign(2, ScalarT{0}); + + auto signal_size = static_cast(RegcaExternalVariables::MAXIMUM); + ws_.assign(signal_size, ScalarT{0}); + ws_indices_.assign(signal_size, INVALID_INDEX); + + for (IdxT j = 0; j < size_; ++j) + { + this->setVariableIndex(j, j); + this->setResidualIndex(j, j); + } + + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y_[static_cast(RegcaInternalVariables::IR)], + &(this->getVariableIndex(static_cast(RegcaInternalVariables::IR)))); + } + + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y_[static_cast(RegcaInternalVariables::II)], + &(this->getVariableIndex(static_cast(RegcaInternalVariables::II)))); + } + + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y_[static_cast(RegcaInternalVariables::PBR)], + &(this->getVariableIndex(static_cast(RegcaInternalVariables::PBR)))); + } + + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y_[static_cast(RegcaInternalVariables::QBR)], + &(this->getVariableIndex(static_cast(RegcaInternalVariables::QBR)))); + } + + return 0; + } + + template + int Regca::verify() const + { + int ret = static_cast(parameter_error_count_); + + auto check = [&](bool condition, const char* message) + { + if (!condition) + { + Log::error() << "Regca: " << message << '\n'; + ret += 1; + } + }; + + if (bus_ == nullptr) + { + Log::error() << "Regca: bus pointer is null\n"; + ret += 1; + } + + check(mva_base_ > ZERO, "mva must be positive"); + check(Rpmax_ > ZERO, "Rpmax must be positive"); + check(Rqmin_ < ZERO && ZERO < Rqmax_, "Rqmin < 0 < Rqmax is required"); + check(IL1_ >= ZERO, "IL1 must be non-negative"); + check(ZERO <= VL0_ && VL0_ < VL1_, "VL0/VL1 must satisfy 0 <= VL0 < VL1"); + check(ZERO <= VA0_ && VA0_ < VA1_, "VA0/VA1 must satisfy 0 <= VA0 < VA1"); + check(Vhvmax_ > ZERO, "Vhvmax must be positive"); + + if (signals_.template isAttached()) + { + if (!signals_.template isLinked()) + { + Log::error() << "Regca: ipcmd signal attached with no linked source\n"; + ret += 1; + } + } + + if (signals_.template isAttached()) + { + if (!signals_.template isLinked()) + { + Log::error() << "Regca: iqcmd signal attached with no linked source\n"; + ret += 1; + } + } + + return ret; + } + + template + int Regca::initialize() + { + if (bus_ == nullptr) + { + Log::error() << "Regca: cannot initialize with null bus\n"; + return 1; + } + + if (parameter_error_count_ > 0 || mva_base_ <= ZERO || Rpmax_ <= ZERO + || !(Rqmin_ < ZERO && ZERO < Rqmax_) + || IL1_ < ZERO || !(ZERO <= VL0_ && VL0_ < VL1_) + || !(ZERO <= VA0_ && VA0_ < VA1_) || Vhvmax_ <= ZERO) + { + Log::error() << "Regca: cannot initialize with invalid parameters\n"; + return 1; + } + + const auto VM = static_cast(RegcaInternalVariables::VM); + const auto IQ = static_cast(RegcaInternalVariables::IQ); + const auto IP = static_cast(RegcaInternalVariables::IP); + const auto VT = static_cast(RegcaInternalVariables::VT); + const auto IR = static_cast(RegcaInternalVariables::IR); + const auto II = static_cast(RegcaInternalVariables::II); + const auto IQEXTRA = static_cast(RegcaInternalVariables::IQEXTRA); + const auto IL = static_cast(RegcaInternalVariables::IL); + const auto LP = static_cast(RegcaInternalVariables::LP); + const auto UP = static_cast(RegcaInternalVariables::UP); + const auto PBR = static_cast(RegcaInternalVariables::PBR); + const auto QBR = static_cast(RegcaInternalVariables::QBR); + + const ScalarT vr = Vr(); + const ScalarT vi = Vi(); + const ScalarT vt = std::sqrt(vr * vr + vi * vi); + + if (vt <= ZERO) + { + Log::error() << "Regca: terminal voltage magnitude must be positive at initialization\n"; + return 1; + } + if (vt >= Vhvmax_) + { + Log::error() + << "Regca: terminal voltage magnitude must be below Vhvmax at initialization\n"; + return 1; + } + + // REGCA owns the network terminal and establishes the initial + // converter operating point from the power-flow injection. Controller + // command ports may not have been initialized yet, so initialization + // resolves the commands from P0/Q0 and publishes them on the system + // base to attached ports below. P0/Q0 are given on the system base; + // toComponentBase converts them to the converter base the internal + // states use. + const ScalarT lvacm = Math::linseg(vt, VA0_, VA1_, ONE); + + if (P0_ != ZERO && (vt <= VA0_ || lvacm <= ZERO) ) + { + Log::error() << "Regca: LVACM gain is zero with nonzero initial active power\n"; + return 1; + } + + ScalarT ipcmd0{ZERO}; + if (P0_ != ZERO) + { + ipcmd0 = toComponentBase(static_cast(P0_) / vt) / lvacm; + } + const ScalarT iqcmd0 = toComponentBase(static_cast(Q0_) / vt); + + const ScalarT iqextra0{ZERO}; + const ScalarT qnet0 = iqcmd0 - iqextra0; + const ScalarT ir0 = (vi * qnet0 + vr * ipcmd0 * lvacm) / vt; + const ScalarT ii0 = (-vr * qnet0 + vi * ipcmd0 * lvacm) / vt; + + y_[VM] = vt; + y_[VT] = vt; + y_[IP] = ipcmd0; + y_[IQ] = iqcmd0; + y_[IQEXTRA] = iqextra0; + y_[IL] = Math::linseg(vt, VL0_, VL1_, IL1_); + y_[IR] = toSystemBase(ir0); + y_[II] = toSystemBase(ii0); + y_[LP] = lpTarget(y_[IP]); + y_[UP] = upTarget(y_[IP], y_[IL]); + y_[PBR] = vr * y_[IR] + vi * y_[II]; + y_[QBR] = vi * y_[IR] - vr * y_[II]; + + // Retain the resolved commands as the constant source used during + // residual evaluation when no controller drives the command ports, and + // select the reactive-current rate-limit branch from the command sign. + ipcmd_set_ = toSystemBase(ipcmd0); + iqcmd_set_ = toSystemBase(iqcmd0); + iq_use_upper_ = ZERO; + iq_use_lower_ = ONE; + if (static_cast(iqcmd0) > ZERO) + { + iq_use_upper_ = ONE; + iq_use_lower_ = ZERO; + } + + // Seed attached command nodes with the steady-state values. Controller + // initialization can use these signal values, and unattached ports fall + // back to the constants stored above. + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(ipcmd_set_); + } + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(iqcmd_set_); + } + + std::fill(yp_.begin(), yp_.end(), ZERO); + return 0; + } + + template + int Regca::tagDifferentiable() + { + std::fill(tag_.begin(), tag_.end(), false); + tag_[static_cast(RegcaInternalVariables::VM)] = true; + tag_[static_cast(RegcaInternalVariables::IQ)] = true; + tag_[static_cast(RegcaInternalVariables::IP)] = true; + return 0; + } + + template + int Regca::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + + template + __attribute__((always_inline)) inline int + Regca::evaluateInternalResidual( + const ScalarT* y, + const ScalarT* yp, + const ScalarT* wb, + const ScalarT* ws, + ScalarT* f) + { + const auto VM = static_cast(RegcaInternalVariables::VM); + const auto IQ = static_cast(RegcaInternalVariables::IQ); + const auto IP = static_cast(RegcaInternalVariables::IP); + const auto VT = static_cast(RegcaInternalVariables::VT); + const auto IR = static_cast(RegcaInternalVariables::IR); + const auto II = static_cast(RegcaInternalVariables::II); + const auto IQEXTRA = static_cast(RegcaInternalVariables::IQEXTRA); + const auto IL = static_cast(RegcaInternalVariables::IL); + const auto LP = static_cast(RegcaInternalVariables::LP); + const auto UP = static_cast(RegcaInternalVariables::UP); + const auto PBR = static_cast(RegcaInternalVariables::PBR); + const auto QBR = static_cast(RegcaInternalVariables::QBR); + + const auto IPCMD = static_cast(RegcaExternalVariables::IPCMD); + const auto IQCMD = static_cast(RegcaExternalVariables::IQCMD); + + const ScalarT vm = y[VM]; + const ScalarT iq = y[IQ]; + const ScalarT ip = y[IP]; + const ScalarT vt = y[VT]; + const ScalarT ir = y[IR]; + const ScalarT ii = y[II]; + const ScalarT iqextra = y[IQEXTRA]; + const ScalarT il = y[IL]; + const ScalarT lp = y[LP]; + const ScalarT up = y[UP]; + const ScalarT pbr = y[PBR]; + const ScalarT qbr = y[QBR]; + + const ScalarT vm_dot = yp[VM]; + const ScalarT iq_dot = yp[IQ]; + const ScalarT ip_dot = yp[IP]; + + const ScalarT vr = wb[0]; + const ScalarT vi = wb[1]; + + const ScalarT ipcmd = toComponentBase(ws[IPCMD]); + const ScalarT iqcmd = toComponentBase(ws[IQCMD]); + + const ScalarT iq_error = iqcmd - iq; + const ScalarT ip_error = ipcmd - ip; + + const ScalarT iq_ramp = + -iq_use_upper_ * Math::ramp(iq_error - Tg_ * Rqmax_) + + iq_use_lower_ * Math::ramp(Tg_ * Rqmin_ - iq_error); + + const ScalarT lp_ramp = Math::ramp(ip_error - Tg_ * lp); + const ScalarT up_ramp = Math::ramp(ip_error - Tg_ * up); + const ScalarT lvacm = Math::linseg(vt, VA0_, VA1_, ONE); + const ScalarT qnet = iq - iqextra; + + f[VM] = -vm_dot + (vt - vm) / TM_; + f[IQ] = -iq_dot + iq_error / Tg_ + iq_ramp / Tg_; + f[IP] = -ip_dot + lp + lp_ramp / Tg_ - up_ramp / Tg_; + f[VT] = -vt * vt + vr * vr + vi * vi; + f[IR] = -vt * ir + toSystemBase(vi * qnet + vr * ip * lvacm); + f[II] = -vt * ii + toSystemBase(-vr * qnet + vi * ip * lvacm); + f[IQEXTRA] = -iqextra + Math::ramp(vt - Vhvmax_); + f[IL] = -il + Math::linseg(vm, VL0_, VL1_, IL1_); + f[LP] = -lp + lpTarget(ip); + f[UP] = -up + upTarget(ip, il); + f[PBR] = -pbr + vr * ir + vi * ii; + f[QBR] = -qbr + vi * ir - vr * ii; + + return 0; + } + + template + __attribute__((always_inline)) inline int Regca::evaluateBusResidual( + const ScalarT* y, + [[maybe_unused]] const ScalarT* yp, + [[maybe_unused]] const ScalarT* wb, + ScalarT* h) + { + const auto IR = static_cast(RegcaInternalVariables::IR); + const auto II = static_cast(RegcaInternalVariables::II); + + h[0] = y[IR]; + h[1] = y[II]; + return 0; + } + + template + int Regca::evaluateResidual() + { + const auto IPCMD = static_cast(RegcaExternalVariables::IPCMD); + const auto IQCMD = static_cast(RegcaExternalVariables::IQCMD); + + ws_[IPCMD] = ipcmd_set_; + ws_[IQCMD] = iqcmd_set_; + std::fill(ws_indices_.begin(), ws_indices_.end(), INVALID_INDEX); + + if (signals_.template isAttached()) + { + ws_[IPCMD] = signals_.template readExternalVariable(); + ws_indices_[IPCMD] = + signals_.template readExternalVariableIndex(); + } + + if (signals_.template isAttached()) + { + ws_[IQCMD] = signals_.template readExternalVariable(); + ws_indices_[IQCMD] = + signals_.template readExternalVariableIndex(); + } + + wb_[0] = Vr(); + wb_[1] = Vi(); + + evaluateInternalResidual(y_.data(), yp_.data(), wb_.data(), ws_.data(), f_.data()); + evaluateBusResidual(y_.data(), yp_.data(), wb_.data(), h_.data()); + + Ir() += h_[0]; + Ii() += h_[1]; + + return 0; + } + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md index 203bc7905..f172ae4b4 100644 --- a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md +++ b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md @@ -145,16 +145,17 @@ are specified: `Gensal` | 5th order salient-pole machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Tdop`, `Tdopp`, `Tqopp`, `Xd`, `Xdp`, `Xdpp`, `Xq`, `Xl`, `S10`, `S12`, `mva` | `ir`, `ii`, `p`, `q`, `delta`, `omega`, `speed`, `Eqp`, `psidp`, `psiqpp`, `psidpp`, `vd`, `vq`, `te`, `id`, `iq` `GenClassical` | the classical machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Xdp`, `mva` | `ir`, `ii`, `p`, `q`, `delta`, `omega` `Tgov1 ` | the TGOV1 governor model | `pmech`, `speed` | `Trate`, `R`, `T1`, `T2`, `T3`, `Pvmax`, `Pvmin`, `Dt` | `none` + `Regca` | WECC REGCA renewable generator/converter model | `bus`, `ipcmd`\*, `iqcmd`\*, `ibranchr`\*, `ibranchi`\*, `pbranch`\*, `qbranch`\* | `P0`, `Q0`, `mva`, `Tg`, `TM`, `Rqmax`, `Rqmin`, `Rpmax`, `sL`, `IL1`, `VL0`, `VL1`, `VA0`, `VA1`, `Vhvmax` | `ir`, `ii`, `p`, `q` + `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` | `A1`, `A2`, `A3`, `A4`, `A5`, `A6`, `T1`, `T2`, `T3`, `T4`, `T5`, `T6`, `Ks`, `Lsmin`, `Lsmax`, `Vcl`, `Vcu`, `Tdelay` | `vss` + `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` `BusToSignalAdapter` | signal adapter component for a bus | `bus`, `vr`, `vi`, `ir`, `ii` | | Ports marked with \* are optional and, if missing, will be assumed to be connected to a constant value. This list is subject to change. - For `Branch`, `tap` and `phase` are optional parameters. If omitted, `tap` defaults to `1.0` and `phase` defaults to `0.0` radians. Bus `bus1` is the tap side for off-nominal transformer branches. diff --git a/GridKit/Model/PhasorDynamics/SystemModelData.hpp b/GridKit/Model/PhasorDynamics/SystemModelData.hpp index 182b45f04..7dbdeed11 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelData.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelData.hpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -41,6 +42,7 @@ namespace GridKit using BusDataT = BusData; using BusToSignalAdapterDataT = BusToSignalAdapterData; using BusFaultDataT = BusFaultData; + using RegcaDataT = Converter::RegcaData; using Tgov1DataT = Governor::Tgov1Data; using Ieeet1DataT = Exciter::Ieeet1Data; using SexsPtiDataT = Exciter::SexsPtiData; @@ -106,6 +108,21 @@ namespace GridKit std::vector stabilizer; ///< Stabilizers within the model std::vector constant_source; ///< Constant signal sources within the model std::vector signal; ///< Signal nodes + std::vector bus; ///< Buses within the model + std::vector adapter; ///< bus-to-signal adapters within the model + std::vector branch; ///< Branches within the model + std::vector bus_fault; ///< Bus faults within the model + std::vector regca; ///< REGCA converter instances within the model + std::vector genrou; ///< GENROU instances within the model + std::vector gensal; ///< GENSAL instances within the model + std::vector genclassical; ///< Classical generator instances within the model + std::vector loadz; ///< LoadZ instances within the model + std::vector loadzip; ///< LoadZIP instances within the model + std::vector gov; ///< Governors within the model + std::vector exciter; ///< Exciters within the model + std::vector sexspti; ///< SEXS-PTI exciters within the model + std::vector stabilizer; ///< Stabilizers within the model + std::vector signal; ///< Signal nodes /// Monitor sink specs std::vector monitor_sink; diff --git a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp index 8f1f08987..37049740c 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp @@ -133,6 +133,12 @@ namespace GridKit raw_component.get_to(loadzip); sm.loadzip.push_back(loadzip); } + else if (kind == "Regca") + { + typename SystemModelData::RegcaDataT regca; + raw_component.get_to(regca); + sm.regca.push_back(regca); + } else if (kind == "Tgov1") { typename SystemModelData::Tgov1DataT gov; diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index b020f8f72..cdda9abff 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -42,6 +42,7 @@ namespace GridKit using namespace Governor; using namespace Exciter; using namespace Stabilizer; + using namespace Converter; owns_components_ = true; @@ -105,6 +106,62 @@ namespace GridKit addComponent(adapter); } + // Add REGCA converters + for (const auto& regcadata : data.regca) + { + IdxT bus_index = 0; + if (regcadata.ports.contains(RegcaPorts::bus)) + { + bus_index = regcadata.ports.at(RegcaPorts::bus); + } + + auto* regca = new Regca(getBus(bus_index), regcadata); + + if (regcadata.ports.contains(RegcaPorts::ipcmd)) + { + IdxT ipcmd = regcadata.ports.at(RegcaPorts::ipcmd); + constexpr auto IPCMD = RegcaExternalVariables::IPCMD; + regca->getSignals().template attachSignalNode(getSignal(ipcmd)); + } + + if (regcadata.ports.contains(RegcaPorts::iqcmd)) + { + IdxT iqcmd = regcadata.ports.at(RegcaPorts::iqcmd); + constexpr auto IQCMD = RegcaExternalVariables::IQCMD; + regca->getSignals().template attachSignalNode(getSignal(iqcmd)); + } + + if (regcadata.ports.contains(RegcaPorts::ibranchr)) + { + IdxT ibranchr = regcadata.ports.at(RegcaPorts::ibranchr); + constexpr auto IR = RegcaInternalVariables::IR; + regca->getSignals().template assignSignalNode(getSignal(ibranchr)); + } + + if (regcadata.ports.contains(RegcaPorts::ibranchi)) + { + IdxT ibranchi = regcadata.ports.at(RegcaPorts::ibranchi); + constexpr auto II = RegcaInternalVariables::II; + regca->getSignals().template assignSignalNode(getSignal(ibranchi)); + } + + if (regcadata.ports.contains(RegcaPorts::pbranch)) + { + IdxT pbranch = regcadata.ports.at(RegcaPorts::pbranch); + constexpr auto PBR = RegcaInternalVariables::PBR; + regca->getSignals().template assignSignalNode(getSignal(pbranch)); + } + + if (regcadata.ports.contains(RegcaPorts::qbranch)) + { + IdxT qbranch = regcadata.ports.at(RegcaPorts::qbranch); + constexpr auto QBR = RegcaInternalVariables::QBR; + regca->getSignals().template assignSignalNode(getSignal(qbranch)); + } + + addComponent(regca); + } + // Add branches for (const auto& branchdata : data.branch) { diff --git a/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp b/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp index 1ce2d8b1f..9c85c6388 100644 --- a/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp +++ b/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp @@ -55,7 +55,31 @@ namespace GridKit return success.report(__func__); } - TestOutcome deadband() + TestOutcome deadband1() + { + TestStatus success = true; + + const ScalarT lower = scalar(-0.05); + const ScalarT upper = scalar(0.10); + const ScalarT tol = scalar(kSmoothTolerance); + + // Type 1 (no-offset) deadband returns ~0 inside the band ... + const ScalarT inside_input = scalar(0.02); + success *= (std::abs(Math::deadband1(inside_input, lower, upper)) < tol * tol); + + // ... and passes the input through unchanged (no offset) outside the band. + const ScalarT far_above_input = scalar(4.0); + const ScalarT far_below_input = scalar(-4.0); + + success *= std::isfinite(Math::deadband1(far_above_input, lower, upper)); + success *= within(Math::deadband1(far_above_input, lower, upper), far_above_input, tol); + success *= std::isfinite(Math::deadband1(far_below_input, lower, upper)); + success *= within(Math::deadband1(far_below_input, lower, upper), far_below_input, tol); + + return success.report(__func__); + } + + TestOutcome deadband2() { TestStatus success = true; @@ -69,12 +93,12 @@ namespace GridKit 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); + success *= within(Math::deadband2(below_input, lower, upper), expected_below, tol); + success *= (std::abs(Math::deadband2(inside_input, lower, upper)) < tol * tol); + success *= within(Math::deadband2(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); + const ScalarT lower_breakpoint = Math::deadband2(lower, lower, upper); + const ScalarT upper_breakpoint = Math::deadband2(upper, lower, upper); success *= (lower_breakpoint < scalar(0.0)); success *= (upper_breakpoint > scalar(0.0)); @@ -82,7 +106,7 @@ namespace GridKit success *= (std::abs(upper_breakpoint) < tol); const ScalarT x = scalar(-0.4); - success *= (std::abs(Math::deadband(x, lower, upper) + success *= (std::abs(Math::deadband2(x, lower, upper) - (x - Math::clamp(x, lower, upper))) < scalar(kRoundoffTolerance)); @@ -91,17 +115,17 @@ namespace GridKit 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); + success *= std::isfinite(Math::deadband2(far_above_input, lower, upper)); + success *= within(Math::deadband2(far_above_input, lower, upper), expected_far_above, tol); + success *= std::isfinite(Math::deadband2(far_below_input, lower, upper)); + success *= within(Math::deadband2(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)) + success *= (std::abs(Math::deadband2(above_point, point, point) - (above_point - point)) < scalar(kRoundoffTolerance)); - success *= (std::abs(Math::deadband(below_point, point, point) - (below_point - point)) + success *= (std::abs(Math::deadband2(below_point, point, point) - (below_point - point)) < scalar(kRoundoffTolerance)); return success.report(__func__); diff --git a/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp b/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp index 233c36ddd..aecdc7ae8 100644 --- a/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp +++ b/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp @@ -7,7 +7,8 @@ int main() GridKit::Testing::SmoothnessIndicatorTests test; result += test.clamp(); - result += test.deadband(); + result += test.deadband1(); + result += test.deadband2(); result += test.limitIndicators(); result += test.slew(); result += test.linseg(); diff --git a/tests/UnitTests/PhasorDynamics/CMakeLists.txt b/tests/UnitTests/PhasorDynamics/CMakeLists.txt index 3d2973b62..4eabe21b5 100644 --- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt +++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt @@ -104,6 +104,14 @@ target_link_libraries( GridKit::phasor_dynamics_systemmodel_dependency_tracking GridKit::testing) +add_executable(test_phasor_converter_regca runConverterRegcaTests.cpp) +target_link_libraries( + test_phasor_converter_regca + GridKit::definitions + GridKit::phasor_dynamics_systemmodel + GridKit::phasor_dynamics_systemmodel_dependency_tracking + GridKit::testing) + add_executable(test_phasor_stabilizer_ieeest runStabilizerIeeestTests.cpp) target_link_libraries( test_phasor_stabilizer_ieeest @@ -149,6 +157,7 @@ add_test(NAME PhasorDynamicsGovernorTgov1Test COMMAND test_phasor_governor_tgov1 add_test(NAME PhasorDynamicsExciterIeeet1Test COMMAND test_phasor_exciter_ieeet1) add_test(NAME PhasorDynamicsGensalTest COMMAND test_phasor_gensal) add_test(NAME PhasorDynamicsExciterSexsPtiTest COMMAND test_phasor_exciter_sexspti) +add_test(NAME PhasorDynamicsConverterRegcaTest COMMAND test_phasor_converter_regca) add_test(NAME PhasorDynamicsStabilizerIeeestTest COMMAND test_phasor_stabilizer_ieeest) add_test(NAME PhasorDynamicsGenClassicalTest COMMAND test_phasor_gen_classical) add_test(NAME PhasorDynamicsLoadZTest COMMAND test_phasor_loadz) @@ -171,6 +180,7 @@ install( test_phasor_exciter_ieeet1 test_phasor_gensal test_phasor_exciter_sexspti + test_phasor_converter_regca test_phasor_stabilizer_ieeest test_phasor_gen_classical test_phasor_system diff --git a/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp new file mode 100644 index 000000000..e57915e81 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp @@ -0,0 +1,1026 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + template + class ConverterRegcaTests + { + public: + using ScalarT = scalar_type; + using IdxT = index_type; + using RealT = typename PhasorDynamics::Component::RealT; + + ConverterRegcaTests() = default; + ~ConverterRegcaTests() = default; + + static constexpr ScalarT kTol = static_cast(1.0e-9); + + TestOutcome constructionAndValidation() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + + PhasorDynamics::Converter::Regca minimal(&bus); + success *= (minimal.size() == static_cast(Vars::MAXIMUM)); + success *= (minimal.getMonitor() == nullptr); + success *= (minimal.verify() > 0); + + PhasorDynamics::Converter::Regca configured(&bus, makeData()); + success *= (configured.size() == static_cast(Vars::MAXIMUM)); + success *= (configured.getMonitor() != nullptr); + success *= (configured.verify() == 0); + + return success.report(__func__); + } + + TestOutcome parameterValidation() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + + auto missing = makeData(); + missing.parameters.erase(Params::Tg); + PhasorDynamics::Converter::Regca missing_model(&bus, missing); + success *= (missing_model.verify() > 0); + + auto missing_p0 = makeData(); + missing_p0.parameters.erase(Params::P0); + PhasorDynamics::Converter::Regca missing_p0_model(&bus, missing_p0); + success *= (missing_p0_model.verify() > 0); + + auto missing_q0 = makeData(); + missing_q0.parameters.erase(Params::Q0); + PhasorDynamics::Converter::Regca missing_q0_model(&bus, missing_q0); + success *= (missing_q0_model.verify() > 0); + + auto bad_switch = makeData(); + bad_switch.parameters[Params::sL] = static_cast(2); + PhasorDynamics::Converter::Regca bad_switch_model(&bus, bad_switch); + success *= (bad_switch_model.verify() > 0); + + success *= invalidParameterCase(bus, Params::mva, static_cast(0.0)); + success *= invalidParameterCase(bus, Params::Rpmax, static_cast(0.0)); + success *= invalidParameterCase(bus, Params::Rqmin, static_cast(0.0)); + success *= invalidParameterCase(bus, Params::IL1, static_cast(-0.1)); + success *= invalidParameterCase(bus, Params::VL1, static_cast(0.3)); + success *= invalidParameterCase(bus, Params::VA1, static_cast(0.3)); + success *= invalidParameterCase(bus, Params::Vhvmax, static_cast(0.0)); + + auto zero_time_constants = makeData(); + zero_time_constants.parameters[Params::Tg] = static_cast(0.0); + zero_time_constants.parameters[Params::TM] = static_cast(0.0); + PhasorDynamics::Converter::Regca zero_time_model( + &bus, + zero_time_constants); + success *= (zero_time_model.verify() == 0); + + return success.report(__func__); + } + + TestOutcome initializesFromPowerFlowAndPublishesSignals() + { + TestStatus success = true; + + const ScalarT vr{0.8}; + const ScalarT vi{0.6}; + const ScalarT p0{0.4}; + const ScalarT q0{-0.1}; + const RealT mva{50.0}; + + PhasorDynamics::Bus bus(vr, vi); + bus.allocate(); + bus.initialize(); + + auto data = makeData(); + data.parameters[Params::mva] = mva; + data.parameters[Params::P0] = static_cast(p0); + data.parameters[Params::Q0] = static_cast(q0); + + PhasorDynamics::Converter::Regca regca(&bus, data); + + ScalarT ipcmd_value{-1.0}; + ScalarT iqcmd_value{1.0}; + IdxT ipcmd_index = 20; + IdxT iqcmd_index = 21; + + PhasorDynamics::SignalNode ipcmd_node; + PhasorDynamics::SignalNode iqcmd_node; + PhasorDynamics::SignalNode ir_node; + PhasorDynamics::SignalNode ii_node; + PhasorDynamics::SignalNode pbranch_node; + PhasorDynamics::SignalNode qbranch_node; + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + + regca.getSignals().template attachSignalNode(&ipcmd_node); + regca.getSignals().template attachSignalNode(&iqcmd_node); + regca.getSignals().template assignSignalNode(&ir_node); + regca.getSignals().template assignSignalNode(&ii_node); + regca.getSignals().template assignSignalNode(&pbranch_node); + regca.getSignals().template assignSignalNode(&qbranch_node); + + success *= (regca.allocate() == 0); + success *= (regca.verify() == 0); + success *= (regca.initialize() == 0); + success *= (regca.evaluateResidual() == 0); + + const ScalarT vt = std::sqrt(vr * vr + vi * vi); + const ScalarT lvacm = + Math::linseg(vt, static_cast(0.4), static_cast(0.9), ONE); + const ScalarT ipcmd0 = toComponentBase(p0 / vt, mva) / lvacm; + const ScalarT iqcmd0 = toComponentBase(q0 / vt, mva); + const ScalarT qnet0 = iqcmd0; + const ScalarT ir0 = (vr * ipcmd0 * lvacm + vi * qnet0) / vt; + const ScalarT ii0 = (vi * ipcmd0 * lvacm - vr * qnet0) / vt; + const ScalarT ipcmd_signal = toSystemBase(ipcmd0, mva); + const ScalarT iqcmd_signal = toSystemBase(iqcmd0, mva); + const ScalarT ir_signal = toSystemBase(ir0, mva); + const ScalarT ii_signal = toSystemBase(ii0, mva); + + success *= scalarMatches(regca.y()[index(Vars::VM)], vt, "VM"); + success *= scalarMatches(regca.y()[index(Vars::VT)], vt, "VT"); + success *= scalarMatches(regca.y()[index(Vars::IP)], ipcmd0, "IP"); + success *= scalarMatches(regca.y()[index(Vars::IQ)], iqcmd0, "IQ"); + success *= scalarMatches(regca.y()[index(Vars::IR)], ir_signal, "IR"); + success *= scalarMatches(regca.y()[index(Vars::II)], ii_signal, "II"); + success *= scalarMatches(regca.y()[index(Vars::PBR)], p0, "PBR"); + success *= scalarMatches(regca.y()[index(Vars::QBR)], q0, "QBR"); + success *= scalarMatches(ipcmd_node.read(), ipcmd_signal, "ipcmd signal"); + success *= scalarMatches(iqcmd_node.read(), iqcmd_signal, "iqcmd signal"); + success *= scalarMatches(ir_node.read(), ir_signal, "ir signal"); + success *= scalarMatches(ii_node.read(), ii_signal, "ii signal"); + success *= scalarMatches(pbranch_node.read(), p0, "pbranch signal"); + success *= scalarMatches(qbranch_node.read(), q0, "qbranch signal"); + + success *= allResidualsZero(regca); + return success.report(__func__); + } + + TestOutcome baseSignals() + { + TestStatus success = true; + + const ScalarT p0{0.4}; + const ScalarT q0{-0.1}; + const RealT mva{50.0}; + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + auto data = makeData(); + data.parameters[Params::mva] = mva; + data.parameters[Params::P0] = static_cast(p0); + data.parameters[Params::Q0] = static_cast(q0); + + ScalarT ipcmd_value{99.0}; + ScalarT iqcmd_value{99.0}; + IdxT ipcmd_index = 30; + IdxT iqcmd_index = 31; + + PhasorDynamics::SignalNode ipcmd_node; + PhasorDynamics::SignalNode iqcmd_node; + PhasorDynamics::SignalNode ir_node; + PhasorDynamics::SignalNode ii_node; + PhasorDynamics::SignalNode pbranch_node; + PhasorDynamics::SignalNode qbranch_node; + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + + PhasorDynamics::Converter::Regca regca(&bus, data); + regca.getSignals().template attachSignalNode(&ipcmd_node); + regca.getSignals().template attachSignalNode(&iqcmd_node); + regca.getSignals().template assignSignalNode(&ir_node); + regca.getSignals().template assignSignalNode(&ii_node); + regca.getSignals().template assignSignalNode(&pbranch_node); + regca.getSignals().template assignSignalNode(&qbranch_node); + + success *= (regca.allocate() == 0); + success *= (regca.verify() == 0); + success *= (regca.initialize() == 0); + success *= (regca.evaluateResidual() == 0); + + success *= scalarMatches(regca.y()[index(Vars::IP)], static_cast(0.8), "IP"); + success *= scalarMatches(regca.y()[index(Vars::IQ)], static_cast(-0.2), "IQ"); + success *= scalarMatches(regca.y()[index(Vars::IR)], static_cast(0.4), "IR"); + success *= scalarMatches(regca.y()[index(Vars::II)], static_cast(0.1), "II"); + success *= scalarMatches(regca.y()[index(Vars::PBR)], p0, "PBR"); + success *= scalarMatches(regca.y()[index(Vars::QBR)], q0, "QBR"); + success *= scalarMatches(ipcmd_node.read(), p0, "ipcmd signal"); + success *= scalarMatches(iqcmd_node.read(), q0, "iqcmd signal"); + success *= scalarMatches(ir_node.read(), static_cast(0.4), "ir signal"); + success *= scalarMatches(ii_node.read(), static_cast(0.1), "ii signal"); + success *= scalarMatches(pbranch_node.read(), p0, "pbranch signal"); + success *= scalarMatches(qbranch_node.read(), q0, "qbranch signal"); + success *= allResidualsZero(regca); + + return success.report(__func__); + } + + TestOutcome unconnectedCommandsRemainConstant() + { + TestStatus success = true; + + auto data = makeData(); + data.parameters[Params::P0] = static_cast(0.6); + data.parameters[Params::Q0] = static_cast(0.2); + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + PhasorDynamics::Converter::Regca regca(&bus, data); + + success *= (regca.allocate() == 0); + success *= (regca.verify() == 0); + success *= (regca.initialize() == 0); + success *= (regca.evaluateResidual() == 0); + + success *= isEqual(regca.y()[index(Vars::IP)], static_cast(0.6), kTol); + success *= isEqual(regca.y()[index(Vars::IQ)], static_cast(0.2), kTol); + success *= allResidualsZero(regca); + + return success.report(__func__); + } + + TestOutcome externalCommandsDriveRuntimeResidual() + { + TestStatus success = true; + + auto data = makeData(); + data.parameters[Params::P0] = static_cast(0.5); + data.parameters[Params::Q0] = static_cast(-0.1); + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + ScalarT ipcmd_value{0.0}; + ScalarT iqcmd_value{0.0}; + IdxT ipcmd_index = 22; + IdxT iqcmd_index = 23; + + PhasorDynamics::SignalNode ipcmd_node; + PhasorDynamics::SignalNode iqcmd_node; + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + + PhasorDynamics::Converter::Regca regca(&bus, data); + regca.getSignals().template attachSignalNode(&ipcmd_node); + regca.getSignals().template attachSignalNode(&iqcmd_node); + + success *= (regca.allocate() == 0); + success *= (regca.verify() == 0); + success *= (regca.initialize() == 0); + + ipcmd_value = static_cast(0.55); + iqcmd_value = static_cast(-0.05); + bus.evaluateResidual(); + success *= (regca.evaluateResidual() == 0); + + const ScalarT tg = static_cast(0.02); + const ScalarT ip_error = ipcmd_value - regca.y()[index(Vars::IP)]; + const ScalarT iq_error = iqcmd_value - regca.y()[index(Vars::IQ)]; + + const ScalarT expected_ip_residual = + regca.y()[index(Vars::LP)] + + Math::ramp(ip_error - tg * regca.y()[index(Vars::LP)]) / tg + - Math::ramp(ip_error - tg * regca.y()[index(Vars::UP)]) / tg; + const ScalarT expected_iq_residual = + iq_error / tg + Math::ramp(tg * static_cast(-999.0) - iq_error) / tg; + + success *= isEqual(regca.getResidual()[index(Vars::IP)], expected_ip_residual, kTol); + success *= isEqual(regca.getResidual()[index(Vars::IQ)], expected_iq_residual, kTol); + success *= (regca.getResidual()[index(Vars::IP)] > ZERO); + success *= (regca.getResidual()[index(Vars::IQ)] > ZERO); + + return success.report(__func__); + } + + TestOutcome invalidInitialization() + { + TestStatus success = true; + + auto data = makeData(); + + { + PhasorDynamics::Bus bus(1.3, 0.0); + bus.allocate(); + bus.initialize(); + + PhasorDynamics::Converter::Regca regca(&bus, data); + success *= (regca.allocate() == 0); + success *= (regca.initialize() > 0); + } + + { + PhasorDynamics::Bus bus(0.0, 0.0); + bus.allocate(); + bus.initialize(); + + PhasorDynamics::Converter::Regca regca(&bus, data); + success *= (regca.allocate() == 0); + success *= (regca.initialize() > 0); + } + + { + auto low_voltage_data = data; + low_voltage_data.parameters[Params::P0] = static_cast(0.5); + + PhasorDynamics::Bus bus(0.2, 0.0); + bus.allocate(); + bus.initialize(); + + PhasorDynamics::Converter::Regca regca(&bus, low_voltage_data); + success *= (regca.allocate() == 0); + success *= (regca.initialize() > 0); + } + + { + auto low_voltage_data = data; + low_voltage_data.parameters[Params::P0] = static_cast(0.0); + low_voltage_data.parameters[Params::Q0] = static_cast(0.1); + + PhasorDynamics::Bus bus(0.2, 0.0); + bus.allocate(); + bus.initialize(); + + PhasorDynamics::Converter::Regca regca(&bus, low_voltage_data); + success *= (regca.allocate() == 0); + success *= (regca.initialize() == 0); + } + + return success.report(__func__); + } + + TestOutcome signalVerification() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + PhasorDynamics::Converter::Regca regca(&bus, makeData()); + + PhasorDynamics::SignalNode ipcmd_node; + PhasorDynamics::SignalNode iqcmd_node; + ScalarT ipcmd_value{0.25}; + ScalarT iqcmd_value{-0.10}; + IdxT ipcmd_index = 24; + IdxT iqcmd_index = 25; + + regca.getSignals().template attachSignalNode(&ipcmd_node); + success *= (regca.verify() > 0); + + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + success *= (regca.verify() == 0); + + regca.getSignals().template attachSignalNode(&iqcmd_node); + success *= (regca.verify() > 0); + + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + success *= (regca.verify() == 0); + + return success.report(__func__); + } + + TestOutcome nullBusVerification() + { + TestStatus success = true; + + PhasorDynamics::Converter::Regca regca(nullptr, makeData()); + success *= (regca.verify() > 0); + success *= (regca.initialize() > 0); + + return success.report(__func__); + } + + TestOutcome busInjectionUsesSystemBase() + { + TestStatus success = true; + + const ScalarT vr{0.8}; + const ScalarT vi{0.6}; + const ScalarT p0{0.4}; + const ScalarT q0{-0.1}; + const RealT mva{50.0}; + + PhasorDynamics::Bus bus(vr, vi); + bus.allocate(); + bus.initialize(); + + auto data = makeData(); + data.parameters[Params::mva] = mva; + data.parameters[Params::P0] = static_cast(p0); + data.parameters[Params::Q0] = static_cast(q0); + + PhasorDynamics::Converter::Regca regca(&bus, data); + success *= (regca.allocate() == 0); + success *= (regca.initialize() == 0); + + bus.evaluateResidual(); + success *= (regca.evaluateResidual() == 0); + + success *= isEqual(bus.Ir(), regca.y()[index(Vars::IR)], kTol); + success *= isEqual(bus.Ii(), regca.y()[index(Vars::II)], kTol); + success *= isEqual(bus.Ir(), static_cast(0.26), kTol); + success *= isEqual(bus.Ii(), static_cast(0.32), kTol); + + return success.report(__func__); + } + + TestOutcome residualEquations() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(0.95, 0.25); + bus.allocate(); + bus.initialize(); + + ScalarT ipcmd_value{0.9}; + ScalarT iqcmd_value{0.1}; + IdxT ipcmd_index = 26; + IdxT iqcmd_index = 27; + + PhasorDynamics::SignalNode ipcmd_node; + PhasorDynamics::SignalNode iqcmd_node; + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + + PhasorDynamics::Converter::Regca regca(&bus, makeDynamicData()); + regca.getSignals().template attachSignalNode(&ipcmd_node); + regca.getSignals().template attachSignalNode(&iqcmd_node); + + success *= (regca.allocate() == 0); + success *= (regca.initialize() == 0); + + setResidualState(regca); + bus.evaluateResidual(); + success *= (regca.evaluateResidual() == 0); + + const auto expected = expectedResidualForState(ipcmd_value, iqcmd_value); + success *= vectorMatches(regca.getResidual(), expected, "REGCA residual"); + + return success.report(__func__); + } + + TestOutcome highVoltageReactiveCurrentRoot() + { + TestStatus success = true; + + auto data = makeDynamicData(); + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + PhasorDynamics::Converter::Regca regca(&bus, data); + success *= (regca.allocate() == 0); + success *= (regca.initialize() == 0); + + regca.y()[index(Vars::VT)] = static_cast(1.35); + regca.y()[index(Vars::IQEXTRA)] = + Math::ramp(regca.y()[index(Vars::VT)] - static_cast(1.3)); + + bus.evaluateResidual(); + success *= (regca.evaluateResidual() == 0); + success *= scalarMatches(regca.getResidual()[index(Vars::IQEXTRA)], + static_cast(0.0), + "HVRCM residual root"); + + return success.report(__func__); + } + + TestOutcome limiterBranchCoverage() + { + TestStatus success = true; + + { + auto data = makeDynamicData(); + data.parameters[Params::Q0] = static_cast(0.2); + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + ScalarT iqcmd_value{0.2}; + IdxT iqcmd_index = 28; + PhasorDynamics::SignalNode iqcmd_node; + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + + PhasorDynamics::Converter::Regca regca(&bus, data); + regca.getSignals().template attachSignalNode(&iqcmd_node); + + success *= (regca.allocate() == 0); + success *= (regca.verify() == 0); + success *= (regca.initialize() == 0); + + iqcmd_value = static_cast(0.4); + bus.evaluateResidual(); + success *= (regca.evaluateResidual() == 0); + + const ScalarT tg = static_cast(0.2); + const ScalarT iq_error = iqcmd_value - regca.y()[index(Vars::IQ)]; + const ScalarT expected = + iq_error / tg - Math::ramp(iq_error - tg * static_cast(0.5)) / tg; + + success *= scalarMatches(regca.getResidual()[index(Vars::IQ)], + expected, + "positive IQ branch"); + } + + { + auto data = makeDynamicData(); + data.parameters[Params::sL] = false; + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + PhasorDynamics::Converter::Regca regca(&bus, data); + success *= (regca.allocate() == 0); + success *= (regca.initialize() == 0); + + bus.evaluateResidual(); + success *= (regca.evaluateResidual() == 0); + + const ScalarT ip = regca.y()[index(Vars::IP)]; + const ScalarT sigma_ip = Math::sigmoid(ip); + const RealT rpmax = static_cast(0.7); + const RealT mp = static_cast(100.0) * rpmax; + const ScalarT expected = mp * (ONE - sigma_ip) + rpmax * sigma_ip; + + success *= scalarMatches(regca.y()[index(Vars::UP)], expected, "UP without LVPL"); + success *= scalarMatches(regca.getResidual()[index(Vars::UP)], + static_cast(0.0), + "UP residual without LVPL"); + } + + return success.report(__func__); + } + + TestOutcome jsonParseAndSystemAssembly() + { + TestStatus success = true; + + std::istringstream input(R"json( +{ + "header": { + "format_version": 0, + "format_revision": 1, + "case_name": "REGCA full model", + "case_description": "REGCA parser behavior test", + "case_comments": "", + "freq_base": 60.0, + "va_base": 100000000.0 + }, + "buses": [ + { + "number": 1, + "class": "bus", + "name": "Bus 1", + "init": { "Vr": 1.0, "Vi": 0.0 }, + "v_base": 1.0 + } + ], + "devices": [ + { + "class": "Regca", + "ports": { "bus": 1 }, + "id": "CV1", + "params": { + "P0": 0.0, + "Q0": 0.0, + "mva": 100, + "Tg": 0.02, + "TM": 0.02, + "Rqmax": 999.0, + "Rqmin": -999.0, + "Rpmax": 999.0, + "sL": true, + "IL1": 1.1, + "VL0": 0.4, + "VL1": 0.9, + "VA0": 0.4, + "VA1": 0.9, + "Vhvmax": 1.2 + }, + "mon": ["ir", "ii", "p", "q"] + } + ] +} +)json"); + + auto data = PhasorDynamics::parseSystemModelData(input); + success *= (data.regca.size() == 1); + const auto& regca_data = data.regca[0]; + success *= (regca_data.device_class == "Regca"); + success *= (regca_data.ports.at(PhasorDynamics::Converter::RegcaPorts::bus) + == 1); + success *= (std::get_if(®ca_data.parameters.at(Params::P0)) + != nullptr); + success *= (std::get_if(®ca_data.parameters.at(Params::Q0)) + != nullptr); + success *= (std::get_if(®ca_data.parameters.at(Params::mva)) + != nullptr); + success *= (std::get_if(®ca_data.parameters.at(Params::sL)) + != nullptr); + + PhasorDynamics::SystemModel system(data); + success *= (system.allocate() == 0); + success *= (system.initialize() == 0); + success *= (system.tagDifferentiable() == 0); + success *= (system.evaluateResidual() == 0); + success *= (system.evaluateJacobian() == 0); + success *= (system.size() == 14); + success *= isEqual(system.getResidual()[0], 0.0, static_cast(kTol)); + success *= isEqual(system.getResidual()[1], 0.0, static_cast(kTol)); + for (size_t i = 2; i < system.getResidual().size(); ++i) + { + success *= isEqual(system.getResidual()[i], 0.0, static_cast(kTol)); + } + + return success.report(__func__); + } + +#ifdef GRIDKIT_ENABLE_ENZYME + TestOutcome jacobian() + { + TestStatus success = true; + + auto dependency_tracking_jacobian = dependencyTrackingJacobian(); + auto enzyme_jacobian = enzymeJacobian(); + + success *= (dependency_tracking_jacobian.size() == enzyme_jacobian.size()); + const auto nrows = std::min(dependency_tracking_jacobian.size(), enzyme_jacobian.size()); + for (size_t i = 0; i < nrows; ++i) + { + success *= isEqual(dependency_tracking_jacobian[i], + enzyme_jacobian[i], + static_cast(1.0e-8)); + } + + return success.report(__func__); + } +#endif + + private: + using Params = PhasorDynamics::Converter::RegcaParameters; + using Vars = PhasorDynamics::Converter::RegcaInternalVariables; + using Ext = PhasorDynamics::Converter::RegcaExternalVariables; + + static size_t index(Vars variable) + { + return static_cast(variable); + } + + auto makeData() -> PhasorDynamics::Converter::RegcaData + { + using Ports = PhasorDynamics::Converter::RegcaPorts; + using Mon = PhasorDynamics::Converter::RegcaMonitorableVariables; + + PhasorDynamics::Converter::RegcaData data; + data.device_class = "Regca"; + data.disambiguation_string = "regca_test"; + data.ports[Ports::bus] = 1; + data.monitored_variables.insert(Mon::ir); + data.monitored_variables.insert(Mon::ii); + data.monitored_variables.insert(Mon::p); + data.monitored_variables.insert(Mon::q); + + data.parameters[Params::P0] = static_cast(0.0); + data.parameters[Params::Q0] = static_cast(0.0); + data.parameters[Params::mva] = static_cast(100.0); + data.parameters[Params::Tg] = static_cast(0.02); + data.parameters[Params::TM] = static_cast(0.02); + data.parameters[Params::Rqmax] = static_cast(999.0); + data.parameters[Params::Rqmin] = static_cast(-999.0); + data.parameters[Params::Rpmax] = static_cast(999.0); + data.parameters[Params::sL] = true; + data.parameters[Params::IL1] = static_cast(1.1); + data.parameters[Params::VL0] = static_cast(0.4); + data.parameters[Params::VL1] = static_cast(0.9); + data.parameters[Params::VA0] = static_cast(0.4); + data.parameters[Params::VA1] = static_cast(0.9); + data.parameters[Params::Vhvmax] = static_cast(1.2); + + return data; + } + + auto makeDynamicData() -> PhasorDynamics::Converter::RegcaData + { + auto data = makeData(); + data.parameters[Params::P0] = static_cast(0.6); + data.parameters[Params::Q0] = static_cast(-0.2); + data.parameters[Params::Tg] = static_cast(0.2); + data.parameters[Params::TM] = static_cast(0.4); + data.parameters[Params::Rqmax] = static_cast(0.5); + data.parameters[Params::Rqmin] = static_cast(-0.6); + data.parameters[Params::Rpmax] = static_cast(0.7); + data.parameters[Params::IL1] = static_cast(1.1); + data.parameters[Params::Vhvmax] = static_cast(1.3); + return data; + } + + bool invalidParameterCase(PhasorDynamics::Bus& bus, Params param, RealT value) + { + auto data = makeData(); + data.parameters[param] = value; + PhasorDynamics::Converter::Regca model(&bus, data); + return model.verify() > 0; + } + + ScalarT toComponentBase(ScalarT value, RealT mva) const + { + return value * static_cast(100.0) / mva; + } + + ScalarT toSystemBase(ScalarT value, RealT mva) const + { + return value * mva / static_cast(100.0); + } + + ScalarT activeLower(ScalarT ip, RealT rpmax) const + { + const RealT mp = static_cast(100.0) * rpmax; + return -rpmax - (mp - rpmax) * Math::sigmoid(ip); + } + + ScalarT activeUpper(ScalarT ip, ScalarT il, RealT rpmax, bool use_lvpl) const + { + const RealT mp = static_cast(100.0) * rpmax; + const ScalarT sigma_ip = Math::sigmoid(ip); + ScalarT output = mp * (ONE - sigma_ip) + rpmax * sigma_ip; + if (use_lvpl) + { + output = mp * (ONE - sigma_ip) + rpmax * sigma_ip * Math::sigmoid(il - ip); + } + return output; + } + + bool allResidualsZero(PhasorDynamics::Converter::Regca& regca) const + { + bool success = true; + for (size_t i = 0; i < static_cast(regca.size()); ++i) + { + if (!isEqual(regca.getResidual()[i], static_cast(0.0), kTol)) + { + std::cout << "REGCA residual row " << i << " is " << regca.getResidual()[i] << "\n"; + success = false; + } + if (!isEqual(regca.yp()[i], static_cast(0.0), kTol)) + { + std::cout << "REGCA derivative row " << i << " is " << regca.yp()[i] << "\n"; + success = false; + } + } + return success; + } + + bool vectorMatches(const std::vector& actual, + const std::vector& expected, + const char* label) const + { + bool success = (actual.size() == expected.size()); + const auto n = std::min(actual.size(), expected.size()); + for (size_t i = 0; i < n; ++i) + { + if (!isEqual(actual[i], expected[i], kTol)) + { + std::cout << label << " mismatch at row " << i << ": " + << actual[i] << " != " << expected[i] << "\n"; + success = false; + } + } + return success; + } + + bool scalarMatches(ScalarT actual, ScalarT expected, const char* label) const + { + if (isEqual(actual, expected, kTol)) + { + return true; + } + + std::cout << label << " mismatch: " << actual << " != " << expected << "\n"; + return false; + } + + void setResidualState(PhasorDynamics::Converter::Regca& regca) + { + regca.y()[index(Vars::VM)] = static_cast(0.86); + regca.y()[index(Vars::IQ)] = static_cast(-0.2); + regca.y()[index(Vars::IP)] = static_cast(0.85); + regca.y()[index(Vars::VT)] = static_cast(0.98); + regca.y()[index(Vars::IR)] = static_cast(0.5); + regca.y()[index(Vars::II)] = static_cast(0.18); + regca.y()[index(Vars::IQEXTRA)] = static_cast(0.03); + regca.y()[index(Vars::IL)] = static_cast(0.72); + regca.y()[index(Vars::LP)] = static_cast(-0.4); + regca.y()[index(Vars::UP)] = static_cast(0.3); + regca.y()[index(Vars::PBR)] = static_cast(0.52); + regca.y()[index(Vars::QBR)] = static_cast(-0.046); + + regca.yp()[index(Vars::VM)] = static_cast(0.01); + regca.yp()[index(Vars::IQ)] = static_cast(-0.02); + regca.yp()[index(Vars::IP)] = static_cast(0.03); + } + + std::vector expectedResidualForState(ScalarT ipcmd, ScalarT iqcmd) const + { + constexpr RealT tg = static_cast(0.2); + constexpr RealT tm = static_cast(0.4); + constexpr RealT rqmin = static_cast(-0.6); + constexpr RealT rpmax = static_cast(0.7); + constexpr RealT vl0 = static_cast(0.4); + constexpr RealT vl1 = static_cast(0.9); + constexpr RealT va0 = static_cast(0.4); + constexpr RealT va1 = static_cast(0.9); + constexpr RealT il1 = static_cast(1.1); + constexpr ScalarT vr = static_cast(0.95); + constexpr ScalarT vi = static_cast(0.25); + constexpr ScalarT vm = static_cast(0.86); + constexpr ScalarT iq = static_cast(-0.2); + constexpr ScalarT ip = static_cast(0.85); + constexpr ScalarT vt = static_cast(0.98); + constexpr ScalarT ii = static_cast(0.18); + constexpr ScalarT iqext = static_cast(0.03); + constexpr ScalarT il = static_cast(0.72); + constexpr ScalarT ir = static_cast(0.5); + constexpr ScalarT lp = static_cast(-0.4); + constexpr ScalarT up = static_cast(0.3); + constexpr ScalarT pbr = static_cast(0.52); + constexpr ScalarT qbr = static_cast(-0.046); + constexpr ScalarT vmdot = static_cast(0.01); + constexpr ScalarT iqdot = static_cast(-0.02); + constexpr ScalarT ipdot = static_cast(0.03); + + const ScalarT iq_error = iqcmd - iq; + const ScalarT ip_error = ipcmd - ip; + const ScalarT lvacm = Math::linseg(vt, va0, va1, ONE); + const ScalarT qnet = iq - iqext; + + std::vector expected(static_cast(Vars::MAXIMUM), + static_cast(0.0)); + expected[index(Vars::VM)] = -vmdot + (vt - vm) / tm; + expected[index(Vars::IQ)] = + -iqdot + iq_error / tg + Math::ramp(tg * rqmin - iq_error) / tg; + expected[index(Vars::IP)] = + -ipdot + lp + Math::ramp(ip_error - tg * lp) / tg + - Math::ramp(ip_error - tg * up) / tg; + expected[index(Vars::VT)] = -vt * vt + vr * vr + vi * vi; + expected[index(Vars::IR)] = -vt * ir + vr * ip * lvacm + vi * qnet; + expected[index(Vars::II)] = -vt * ii + vi * ip * lvacm - vr * qnet; + expected[index(Vars::IQEXTRA)] = -iqext + Math::ramp(vt - static_cast(1.3)); + expected[index(Vars::IL)] = -il + Math::linseg(vm, vl0, vl1, il1); + expected[index(Vars::LP)] = -lp + activeLower(ip, rpmax); + expected[index(Vars::UP)] = -up + activeUpper(ip, il, rpmax, true); + expected[index(Vars::PBR)] = -pbr + vr * ir + vi * ii; + expected[index(Vars::QBR)] = -qbr + vi * ir - vr * ii; + return expected; + } + +#ifdef GRIDKIT_ENABLE_ENZYME + void setResidualStateDep( + PhasorDynamics::Converter::Regca& regca, + PhasorDynamics::Bus& bus) + { + bus.y()[0].setValue(0.95); + bus.y()[1].setValue(0.25); + + regca.y()[index(Vars::VM)].setValue(0.86); + regca.y()[index(Vars::IQ)].setValue(-0.2); + regca.y()[index(Vars::IP)].setValue(0.85); + regca.y()[index(Vars::VT)].setValue(0.98); + regca.y()[index(Vars::IR)].setValue(0.5); + regca.y()[index(Vars::II)].setValue(0.18); + regca.y()[index(Vars::IQEXTRA)].setValue(0.03); + regca.y()[index(Vars::IL)].setValue(0.72); + regca.y()[index(Vars::LP)].setValue(-0.4); + regca.y()[index(Vars::UP)].setValue(0.3); + regca.y()[index(Vars::PBR)].setValue(0.52); + regca.y()[index(Vars::QBR)].setValue(-0.046); + + regca.yp()[index(Vars::VM)].setValue(0.01); + regca.yp()[index(Vars::IQ)].setValue(-0.02); + regca.yp()[index(Vars::IP)].setValue(0.03); + } + + std::vector dependencyTrackingJacobian() + { + using DepVar = DependencyTracking::Variable; + + auto data = makeDynamicData(); + data.parameters[Params::mva] = static_cast(50.0); + + PhasorDynamics::Bus bus(DepVar{0.95}, DepVar{0.25}); + PhasorDynamics::Converter::Regca regca(&bus, data); + + PhasorDynamics::SignalNode ipcmd_node; + PhasorDynamics::SignalNode iqcmd_node; + DepVar ipcmd_value{0.9}; + DepVar iqcmd_value{0.1}; + IdxT ipcmd_index = static_cast(regca.size() + bus.size()); + IdxT iqcmd_index = static_cast(regca.size() + bus.size() + 1); + + bus.allocate(); + regca.allocate(); + bus.initialize(); + setResidualStateDep(regca, bus); + + for (IdxT i = 0; i < regca.size(); ++i) + { + regca.y()[static_cast(i)].setVariableNumber(i); + regca.yp()[static_cast(i)].setVariableNumber(i); + } + for (IdxT i = 0; i < bus.size(); ++i) + { + bus.y()[static_cast(i)].setVariableNumber(i + regca.size()); + } + ipcmd_value.setVariableNumber(ipcmd_index); + iqcmd_value.setVariableNumber(iqcmd_index); + + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + regca.getSignals().template attachSignalNode(&ipcmd_node); + regca.getSignals().template attachSignalNode(&iqcmd_node); + + bus.evaluateResidual(); + regca.evaluateResidual(); + + std::vector dependencies( + static_cast(regca.size() + bus.size())); + for (IdxT i = 0; i < regca.size(); ++i) + { + dependencies[static_cast(i)] = + regca.getResidual()[static_cast(i)].getDependencies(); + } + dependencies[static_cast(regca.size())] = bus.Ir().getDependencies(); + dependencies[static_cast(regca.size() + 1)] = bus.Ii().getDependencies(); + + return dependencies; + } + + std::vector enzymeJacobian() + { + auto data = makeDynamicData(); + data.parameters[Params::mva] = static_cast(50.0); + + PhasorDynamics::Bus bus(0.95, 0.25); + PhasorDynamics::Converter::Regca regca(&bus, data); + + PhasorDynamics::SignalNode ipcmd_node; + PhasorDynamics::SignalNode iqcmd_node; + ScalarT ipcmd_value{0.9}; + ScalarT iqcmd_value{0.1}; + IdxT ipcmd_index = static_cast(regca.size() + bus.size()); + IdxT iqcmd_index = static_cast(regca.size() + bus.size() + 1); + + bus.allocate(); + regca.allocate(); + for (IdxT i = 0; i < bus.size(); ++i) + { + bus.setVariableIndex(i, i + regca.size()); + bus.setResidualIndex(i, i + regca.size()); + } + + bus.initialize(); + setResidualState(regca); + regca.updateTime(0.0, 1.0); + + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + regca.getSignals().template attachSignalNode(&ipcmd_node); + regca.getSignals().template attachSignalNode(&iqcmd_node); + + bus.evaluateResidual(); + regca.evaluateResidual(); + regca.evaluateJacobian(); + + regca.constructCsr(); + return MapFromCsr(regca.getCsrJacobian()); + } +#endif + }; + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp b/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp index 2ffcfbc6d..ce5e46c0b 100644 --- a/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp +++ b/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp @@ -165,6 +165,57 @@ namespace GridKit return success.report(__func__); } + TestOutcome regca() + { + TestStatus success = true; + + PhasorDynamics::SystemModel* system = new PhasorDynamics::SystemModel(); + + PhasorDynamics::BusInfinite bus(1.0, 0.0); + system->addBus(&bus); + + PhasorDynamics::Converter::Regca regca(&bus, makeRegcaData()); + system->addComponent(®ca); + + success *= system->allocate() == 0; + success *= system->initialize() == 0; + success *= system->evaluateResidual() == 0; + success *= system->evaluateJacobian() == 0; + success *= system->size() == regca.size(); + + delete system; + system = nullptr; + + return success.report(__func__); + } + + private: + auto makeRegcaData() -> PhasorDynamics::Converter::RegcaData + { + using Params = PhasorDynamics::Converter::RegcaParameters; + + PhasorDynamics::Converter::RegcaData data; + data.device_class = "Regca"; + data.disambiguation_string = "regca_test"; + data.parameters[Params::P0] = static_cast(1.0); + data.parameters[Params::Q0] = static_cast(0.0); + data.parameters[Params::mva] = static_cast(100.0); + data.parameters[Params::Tg] = static_cast(0.02); + data.parameters[Params::TM] = static_cast(0.02); + data.parameters[Params::Rqmax] = static_cast(999.0); + data.parameters[Params::Rqmin] = static_cast(-999.0); + data.parameters[Params::Rpmax] = static_cast(999.0); + data.parameters[Params::sL] = true; + data.parameters[Params::IL1] = static_cast(1.1); + data.parameters[Params::VL0] = static_cast(0.4); + data.parameters[Params::VL1] = static_cast(0.9); + data.parameters[Params::VA0] = static_cast(0.4); + data.parameters[Params::VA1] = static_cast(0.9); + data.parameters[Params::Vhvmax] = static_cast(1.2); + return data; + } + + public: TestOutcome genrou() { TestStatus success = true; diff --git a/tests/UnitTests/PhasorDynamics/runConverterRegcaTests.cpp b/tests/UnitTests/PhasorDynamics/runConverterRegcaTests.cpp new file mode 100644 index 000000000..77d3449a3 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/runConverterRegcaTests.cpp @@ -0,0 +1,28 @@ +#include "ConverterRegcaTests.hpp" + +int main() +{ + GridKit::Testing::TestingResults result; + + GridKit::Testing::ConverterRegcaTests test; + + result += test.constructionAndValidation(); + result += test.parameterValidation(); + result += test.initializesFromPowerFlowAndPublishesSignals(); + result += test.baseSignals(); + result += test.unconnectedCommandsRemainConstant(); + result += test.externalCommandsDriveRuntimeResidual(); + result += test.invalidInitialization(); + result += test.signalVerification(); + result += test.nullBusVerification(); + result += test.busInjectionUsesSystemBase(); + result += test.residualEquations(); + result += test.highVoltageReactiveCurrentRoot(); + result += test.limiterBranchCoverage(); + result += test.jsonParseAndSystemAssembly(); +#ifdef GRIDKIT_ENABLE_ENZYME + result += test.jacobian(); +#endif + + return result.summary(); +} diff --git a/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp b/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp index fc4f07e80..2eb8cecb8 100644 --- a/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp +++ b/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp @@ -14,6 +14,7 @@ int main() result += test.ieeet1(); result += test.load(); result += test.loadZIP(); + result += test.regca(); result += test.genrou(); result += test.genClassical(); result += test.tgov1(); From ab0a8f2baa5d149ae7a88ce972189081252c4892 Mon Sep 17 00:00:00 2001 From: lukelowry Date: Thu, 2 Jul 2026 19:15:59 -0500 Subject: [PATCH 2/2] update REGCA model --- .../PhasorDynamics/Converter/REGCA/README.md | 18 ++++++++-- .../Converter/REGCA/RegcaData.hpp | 33 ++++++++++++++----- .../Converter/REGCA/RegcaImpl.hpp | 6 ++-- .../Model/PhasorDynamics/SystemModelData.hpp | 16 +-------- .../Model/PhasorDynamics/SystemModelImpl.hpp | 28 ++++++++-------- .../PhasorDynamics/ConverterRegcaTests.hpp | 8 +++-- 6 files changed, 63 insertions(+), 46 deletions(-) diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md b/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md index 951f17dbc..6176bdc72 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md @@ -83,6 +83,18 @@ Invalid REGCA parameter sets are rejected by the following checks. Let $\epsilon $M_p$ is a finite surrogate for inactive active-current rate limits, not a physical REGCA parameter. +## Model Ports + +Name | Port | Init | Description +-----------|--------|---------|------ +`bus` | Bus | Known | Terminal bus voltage +`ipcmd` | Input | Unknown | Active-current command input +`iqcmd` | Input | Unknown | Reactive-current command input +`ibranchr` | Output | Known | Branch-current real-component output +`ibranchi` | Output | Known | Branch-current imaginary-component output +`pbranch` | Output | Known | Branch active-power output +`qbranch` | Output | Known | Branch reactive-power output + ## Model Variables ### Internal Variables @@ -210,11 +222,11 @@ The bus receives the REGCA branch-current variables directly: ## Initialization -### External Priors +### Input Initialization ```math \begin{aligned} - V_{\mathrm{r},0}, V_{\mathrm{i},0} + V_{\mathrm{r}}, V_{\mathrm{i}} &\leftarrow \text{terminal-bus voltage} \\ P_0, Q_0 &\leftarrow \text{power-flow injection on system base} @@ -281,7 +293,7 @@ Here $\text{linseg}$ is GridKit's [linear segment](../../../../CommonMath.md#derived-functions), and $\sigma$ is the [step function](../../../../CommonMath.md#primitives). -### External Solved +### Output Initialization ```math \begin{aligned} diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp index 231343604..65a77c786 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp @@ -34,16 +34,29 @@ namespace GridKit Vhvmax ///< Terminal-voltage ceiling for HV reactive management }; - /// Ports for the REGCA converter model. - enum class RegcaPorts + /// Buses for the REGCA converter model. + enum class RegcaBuses : size_t + { + bus, ///< Terminal bus ID + SIZE + }; + + /// Signal inputs for the REGCA converter model. + enum class RegcaSignalInputs : size_t + { + ipcmd, ///< Optional active-current command signal ID + iqcmd, ///< Optional reactive-current command signal ID + SIZE + }; + + /// Signal outputs for the REGCA converter model. + enum class RegcaSignalOutputs : size_t { - bus, ///< Terminal bus ID - ipcmd, ///< Optional active-current command signal ID - iqcmd, ///< Optional reactive-current command signal ID ibranchr, ///< Optional real current measurement output signal ID ibranchi, ///< Optional imaginary current measurement output signal ID pbranch, ///< Optional active-power measurement output signal ID - qbranch ///< Optional reactive-power measurement output signal ID + qbranch, ///< Optional reactive-power measurement output signal ID + SIZE }; /// Variables available through the monitor interface. @@ -59,13 +72,17 @@ namespace GridKit struct RegcaData : public ComponentData { RegcaData() = default; using Parameters = RegcaParameters; - using Ports = RegcaPorts; + using Buses = RegcaBuses; + using SignalInputs = RegcaSignalInputs; + using SignalOutputs = RegcaSignalOutputs; using MonitorableVariables = RegcaMonitorableVariables; }; } // namespace Converter diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp index ab3fff561..f4c3e50c5 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp @@ -85,7 +85,7 @@ namespace GridKit void Regca::initializeParameters(const ModelDataT& data) { using Params = typename ModelDataT::Parameters; - using Ports = typename ModelDataT::Ports; + using Buses = typename ModelDataT::Buses; parameter_error_count_ = 0; @@ -156,9 +156,9 @@ namespace GridKit load_required_real(Params::VA1, VA1_, "VA1"); load_required_real(Params::Vhvmax, Vhvmax_, "Vhvmax"); - if (data.ports.contains(Ports::bus)) + if (data.buses.contains(Buses::bus)) { - bus_id_ = data.ports.at(Ports::bus); + bus_id_ = data.buses.at(Buses::bus); } setDerivedParameters(); diff --git a/GridKit/Model/PhasorDynamics/SystemModelData.hpp b/GridKit/Model/PhasorDynamics/SystemModelData.hpp index 7dbdeed11..dbe151d04 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelData.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelData.hpp @@ -97,6 +97,7 @@ namespace GridKit std::vector adapter; ///< bus-to-signal adapters within the model std::vector branch; ///< Branches within the model std::vector bus_fault; ///< Bus faults within the model + std::vector regca; ///< REGCA converter instances within the model std::vector genrou; ///< GENROU instances within the model std::vector gensal; ///< GENSAL instances within the model std::vector genclassical; ///< Classical generator instances within the model @@ -108,21 +109,6 @@ namespace GridKit std::vector stabilizer; ///< Stabilizers within the model std::vector constant_source; ///< Constant signal sources within the model std::vector signal; ///< Signal nodes - std::vector bus; ///< Buses within the model - std::vector adapter; ///< bus-to-signal adapters within the model - std::vector branch; ///< Branches within the model - std::vector bus_fault; ///< Bus faults within the model - std::vector regca; ///< REGCA converter instances within the model - std::vector genrou; ///< GENROU instances within the model - std::vector gensal; ///< GENSAL instances within the model - std::vector genclassical; ///< Classical generator instances within the model - std::vector loadz; ///< LoadZ instances within the model - std::vector loadzip; ///< LoadZIP instances within the model - std::vector gov; ///< Governors within the model - std::vector exciter; ///< Exciters within the model - std::vector sexspti; ///< SEXS-PTI exciters within the model - std::vector stabilizer; ///< Stabilizers within the model - std::vector signal; ///< Signal nodes /// Monitor sink specs std::vector monitor_sink; diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index cdda9abff..fc38af80a 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -110,51 +110,51 @@ namespace GridKit for (const auto& regcadata : data.regca) { IdxT bus_index = 0; - if (regcadata.ports.contains(RegcaPorts::bus)) + if (regcadata.buses.contains(RegcaBuses::bus)) { - bus_index = regcadata.ports.at(RegcaPorts::bus); + bus_index = regcadata.buses.at(RegcaBuses::bus); } auto* regca = new Regca(getBus(bus_index), regcadata); - if (regcadata.ports.contains(RegcaPorts::ipcmd)) + if (regcadata.signal_inputs.contains(RegcaSignalInputs::ipcmd)) { - IdxT ipcmd = regcadata.ports.at(RegcaPorts::ipcmd); + IdxT ipcmd = regcadata.signal_inputs.at(RegcaSignalInputs::ipcmd); constexpr auto IPCMD = RegcaExternalVariables::IPCMD; regca->getSignals().template attachSignalNode(getSignal(ipcmd)); } - if (regcadata.ports.contains(RegcaPorts::iqcmd)) + if (regcadata.signal_inputs.contains(RegcaSignalInputs::iqcmd)) { - IdxT iqcmd = regcadata.ports.at(RegcaPorts::iqcmd); + IdxT iqcmd = regcadata.signal_inputs.at(RegcaSignalInputs::iqcmd); constexpr auto IQCMD = RegcaExternalVariables::IQCMD; regca->getSignals().template attachSignalNode(getSignal(iqcmd)); } - if (regcadata.ports.contains(RegcaPorts::ibranchr)) + if (regcadata.signal_outputs.contains(RegcaSignalOutputs::ibranchr)) { - IdxT ibranchr = regcadata.ports.at(RegcaPorts::ibranchr); + IdxT ibranchr = regcadata.signal_outputs.at(RegcaSignalOutputs::ibranchr); constexpr auto IR = RegcaInternalVariables::IR; regca->getSignals().template assignSignalNode(getSignal(ibranchr)); } - if (regcadata.ports.contains(RegcaPorts::ibranchi)) + if (regcadata.signal_outputs.contains(RegcaSignalOutputs::ibranchi)) { - IdxT ibranchi = regcadata.ports.at(RegcaPorts::ibranchi); + IdxT ibranchi = regcadata.signal_outputs.at(RegcaSignalOutputs::ibranchi); constexpr auto II = RegcaInternalVariables::II; regca->getSignals().template assignSignalNode(getSignal(ibranchi)); } - if (regcadata.ports.contains(RegcaPorts::pbranch)) + if (regcadata.signal_outputs.contains(RegcaSignalOutputs::pbranch)) { - IdxT pbranch = regcadata.ports.at(RegcaPorts::pbranch); + IdxT pbranch = regcadata.signal_outputs.at(RegcaSignalOutputs::pbranch); constexpr auto PBR = RegcaInternalVariables::PBR; regca->getSignals().template assignSignalNode(getSignal(pbranch)); } - if (regcadata.ports.contains(RegcaPorts::qbranch)) + if (regcadata.signal_outputs.contains(RegcaSignalOutputs::qbranch)) { - IdxT qbranch = regcadata.ports.at(RegcaPorts::qbranch); + IdxT qbranch = regcadata.signal_outputs.at(RegcaSignalOutputs::qbranch); constexpr auto QBR = RegcaInternalVariables::QBR; regca->getSignals().template assignSignalNode(getSignal(qbranch)); } diff --git a/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp index e57915e81..1efe02f7e 100644 --- a/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp @@ -635,8 +635,10 @@ namespace GridKit success *= (data.regca.size() == 1); const auto& regca_data = data.regca[0]; success *= (regca_data.device_class == "Regca"); - success *= (regca_data.ports.at(PhasorDynamics::Converter::RegcaPorts::bus) + success *= (regca_data.buses.at(PhasorDynamics::Converter::RegcaBuses::bus) == 1); + success *= regca_data.signal_inputs.empty(); + success *= regca_data.signal_outputs.empty(); success *= (std::get_if(®ca_data.parameters.at(Params::P0)) != nullptr); success *= (std::get_if(®ca_data.parameters.at(Params::Q0)) @@ -696,13 +698,13 @@ namespace GridKit auto makeData() -> PhasorDynamics::Converter::RegcaData { - using Ports = PhasorDynamics::Converter::RegcaPorts; + using Buses = PhasorDynamics::Converter::RegcaBuses; using Mon = PhasorDynamics::Converter::RegcaMonitorableVariables; PhasorDynamics::Converter::RegcaData data; data.device_class = "Regca"; data.disambiguation_string = "regca_test"; - data.ports[Ports::bus] = 1; + data.buses[Buses::bus] = 1; data.monitored_variables.insert(Mon::ir); data.monitored_variables.insert(Mon::ii); data.monitored_variables.insert(Mon::p);