diff --git a/CHANGELOG.md b/CHANGELOG.md index e552d2cf3..4141999ae 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -78,6 +78,7 @@ - Remove unnecessary data copying while evaluating `PowerElectronics` models, speeding up large simulations by up to 3x - Added `HYGOV` governor model implementation for PhasorDynamics. - Added `REPCA` controller model implementation for PhasorDynamics. +- Added `REECB` electrical-control model implementation for PhasorDynamics. ## v0.1 diff --git a/GridKit/CommonMath.md b/GridKit/CommonMath.md index 74e79289e..33b9f073f 100644 --- a/GridKit/CommonMath.md +++ b/GridKit/CommonMath.md @@ -67,8 +67,8 @@ q(x)=x^2\,\sigma(x) | `linseg` | Saturated linear segment contribution | `REGCA`, `REECA` | | `above` | Above-lower-limit indicator | `REPCA` | | `below` | Below-upper-limit indicator | - | -| `inside` | Interior pulse indicator | - | -| `outside` | Outside-band indicator | `REECA`, `REECB` | +| `inside` | Interior pulse indicator | `REECB` | +| `outside` | Outside-band indicator | `REECA` | | `antiwindup` | Anti-windup limited derivative | `IEEET1`, `SEXS-PTI`, `TGOV1`, `REECA`, `REECB`, `REPCA` | ### `max` diff --git a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp index 21b7210ff..4d2d81aa5 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/Controller/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Controller/CMakeLists.txt index be3287ed4..d3ff1437e 100644 --- a/GridKit/Model/PhasorDynamics/Controller/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/Controller/CMakeLists.txt @@ -3,4 +3,5 @@ # - Luke Lowery # ]] +add_subdirectory(REECB) add_subdirectory(REPCA) diff --git a/GridKit/Model/PhasorDynamics/Controller/README.md b/GridKit/Model/PhasorDynamics/Controller/README.md index ce5f57486..aac18d485 100644 --- a/GridKit/Model/PhasorDynamics/Controller/README.md +++ b/GridKit/Model/PhasorDynamics/Controller/README.md @@ -7,4 +7,5 @@ directly contributing to the network equations. ## Types +- Renewable Energy Electrical Control Model REECB (See [REECB](REECB/README.md)) - Renewable Energy Plant Control Model REPCA (See [REPCA](REPCA/README.md)) diff --git a/GridKit/Model/PhasorDynamics/Controller/REECB/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Controller/REECB/CMakeLists.txt new file mode 100644 index 000000000..50f8862f3 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Controller/REECB/CMakeLists.txt @@ -0,0 +1,59 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +set(_install_headers Reecb.hpp ReecbData.hpp) + +if(GRIDKIT_ENABLE_ENZYME) + gridkit_add_library( + phasor_dynamics_controller_reecb + SOURCES ReecbEnzyme.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) + + if(CMAKE_CXX_COMPILER_ID STREQUAL "Clang" AND CMAKE_CXX_COMPILER_VERSION VERSION_GREATER_EQUAL 19) + # Work around EnzymeAD/Enzyme#3101. + target_compile_options(phasor_dynamics_controller_reecb PRIVATE -fno-builtin-tan) + endif() +else() + gridkit_add_library( + phasor_dynamics_controller_reecb + SOURCES Reecb.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_controller_reecb_dependency_tracking + SOURCES ReecbDependencyTracking.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_controller_reecb) +target_link_libraries( + phasor_dynamics_components_dependency_tracking + INTERFACE GridKit::phasor_dynamics_controller_reecb_dependency_tracking) diff --git a/GridKit/Model/PhasorDynamics/Controller/REECB/README.md b/GridKit/Model/PhasorDynamics/Controller/REECB/README.md new file mode 100644 index 000000000..2de5597f9 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Controller/REECB/README.md @@ -0,0 +1,404 @@ +# **Renewable Energy Electrical Control Model (REECB)** + +REECB is a WECC renewable electrical-control model with power-factor, +reactive-power, voltage, and active-power command paths for an +inverter-coupled resource. + +## Notes + +- Current commands and power signals are on system base. +- Internal control states and the reactive-power, active-power, and current + limits are on REECB component base. +- REECB uses `mva` as its component power base. +- In direct-voltage mode ($s_Q=1$, $s_V=0$) `qext` carries a terminal-voltage + reference instead of a system-base reactive power. +- REECB contributes no bus current injection. +- GridKit does not apply generator-level Governor Response Limits to `Pmin` or + `Pmax`. + +> [!WARNING] +> GridKit does not yet inherit `mva` from the associated REGCA model. Set it +> explicitly to the REGCA component base; omitting it falls back to the system +> base and is correct only when those bases match.[^reecb-mva-base] + +## Block Diagram + +![REECB electrical-control block diagram](../../../../../docs/Figures/PhasorDynamics/REECB/diagram.png) + +Figure 1: REECB electrical-control model. Figure courtesy of the +[PowerWorld REEC_B model reference](https://www.powerworld.com/WebHelp/Content/TransientModels_HTML/Exciter%20REEC_B.htm). + +## Model Parameters + +Symbol | Units | JSON | Description | Typical Value | Note +------------------------------------|-----------|----------|---------------------------------------------------------|---------------|----- +$S^\mathrm{base}$ | [MVA] | `mva` | REECB component power base | 100.0 | System power base when omitted +$s_\mathrm{pf}$ | [boolean] | `PfFlag` | Power-factor control selector | `false` | `true` = power-factor control, `false` = reactive-power control +$s_V$ | [boolean] | `VFlag` | Voltage-reference selector under $s_Q=1$ | `false` | `true` = cascaded Q-PI voltage command, `false` = direct external voltage reference +$s_Q$ | [boolean] | `QFlag` | Reactive-path selector | `false` | `true` = Volt/VAr PI control, `false` = reactive-current lag +$s_\mathrm{pq}$ | [boolean] | `Pqflag` | Converter current-priority selector | `false` | `true` = P priority, `false` = Q priority +$T_\mathrm{rv}$ | [sec] | `Trv` | Voltage-measurement filter time constant | 0.02 | State 1 in Fig. 1 +$T_\mathrm{p}$ | [sec] | `Tp` | Electrical-power measurement filter time constant | 0.0 | State 2 in Fig. 1 +$V^\mathrm{ref}$ | [p.u.] | `Vref0` | Reactive-current-injection voltage reference | $V_T$ | Initialized from terminal voltage when omitted +$V_\mathrm{dip}$ | [p.u.] | `Vdip` | Low-voltage threshold for the voltage-band gate | 0.85 | +$V_\mathrm{up}$ | [p.u.] | `Vup` | High-voltage threshold for the voltage-band gate | 1.15 | +$D_1^\mathrm{db}$ | [p.u.] | `dbd1` | Lower deadband threshold for voltage-error response | 0.0 | +$D_2^\mathrm{db}$ | [p.u.] | `dbd2` | Upper deadband threshold for voltage-error response | 0.0 | +$K_\mathrm{qv}$ | [p.u.] | `kqv` | Reactive-current injection gain | 5.0 | +$I_{q,\mathrm{inj}}^{\min}$ | [p.u.] | `Iql1` | Minimum reactive-current injection | -1.1 | +$I_{q,\mathrm{inj}}^{\max}$ | [p.u.] | `Iqh1` | Maximum reactive-current injection | 1.1 | +$Q^{\max}$ | [p.u.] | `Qmax` | Maximum reactive-power control output | 0.436 | +$Q^{\min}$ | [p.u.] | `Qmin` | Minimum reactive-power control output | -0.436 | +$K_\mathrm{qp}$ | [p.u.] | `Kqp` | Reactive-power controller proportional gain | 0.0 | +$K_\mathrm{qi}$ | [p.u./s] | `Kqi` | Reactive-power controller integral gain | 0.1 | +$V^{\max}$ | [p.u.] | `Vmax` | Maximum voltage-control output | 1.1 | +$V^{\min}$ | [p.u.] | `Vmin` | Minimum voltage-control output | 0.9 | +$K_\mathrm{vp}$ | [p.u.] | `Kvp` | Voltage controller proportional gain | 18.0 | +$K_\mathrm{vi}$ | [p.u./s] | `Kvi` | Voltage controller integral gain | 5.0 | +$T_\mathrm{iq}$ | [sec] | `Tiq` | Reactive-current command lag time constant | 0.02 | State 5 in Fig. 1 +$T_\mathrm{pord}$ | [sec] | `Tpord` | Active-power order filter time constant | 0.02 | State 6 in Fig. 1 +$R_P^{\max}$ | [p.u./s] | `dPmax` | Positive active-power order ramp-rate limit | 99.0 | +$R_P^{\min}$ | [p.u./s] | `dPmin` | Negative active-power order ramp-rate limit | -99.0 | +$P^{\max}$ | [p.u.] | `Pmax` | Maximum active-power order | 1.0 | +$P^{\min}$ | [p.u.] | `Pmin` | Minimum active-power order | 0.0 | +$I^{\max}$ | [p.u.] | `Imax` | Maximum converter current | 1.3 | + +All parameters are optional. An omitted parameter starts from its Typical +Value; the time-constant floor below is then applied. Real-valued parameters +accept real or integer JSON values; selectors require Boolean JSON values. + +### Parameter Validation + +Invalid REECB parameter sets are rejected by the following checks: + +```math +\begin{aligned} + S^\mathrm{base} &> 0,\quad \text{when provided} \\ + T_\mathrm{rv},T_\mathrm{p},T_\mathrm{iq},T_\mathrm{pord} &\ge 0 \\ + V_\mathrm{dip} &< V_\mathrm{up} \\ + D_1^\mathrm{db} &\le 0 \le D_2^\mathrm{db} \\ + K_\mathrm{qv},K_\mathrm{qp},K_\mathrm{qi},K_\mathrm{vp},K_\mathrm{vi} &\ge 0 \\ + I_{q,\mathrm{inj}}^{\min} &\le I_{q,\mathrm{inj}}^{\max} \\ + Q^{\min} &\le Q^{\max} \\ + V^{\min} &\le V^{\max} \\ + R_P^{\min} &< 0 < R_P^{\max} \\ + P^{\min} &\le P^{\max} \\ + I^{\max} &> 0. +\end{aligned} +``` + +Enabling both `PfFlag` and `QFlag` logs an atypical-configuration warning. + +### Model Derived Parameters + +Let $\epsilon_T=10^{-3}\ \mathrm{s}$. A time constant below $\epsilon_T$ is +raised to that floor in place, so every equation below uses the raised value: + +```math +\begin{aligned} + T_x &\leftarrow \text{max}(T_x,\epsilon_T), && x\in\{\mathrm{rv},\mathrm{p},\mathrm{iq},\mathrm{pord}\} \\ + s_\mathrm{pf}^\mathrm{off} &= 1 - s_\mathrm{pf} \\ + s_Q^\mathrm{off} &= 1 - s_Q \\ + s_Q^\mathrm{PI} &= s_Q s_V \\ + s_V^\mathrm{ref} &= s_Q(1-s_V) \\ + s_Q^\mathrm{ref} &= 1 - s_V^\mathrm{ref} \\ + s_\mathrm{pq}^\mathrm{off} &= 1 - s_\mathrm{pq} \\ + k_\mathrm{base} &= \dfrac{S^\mathrm{sys}}{S^\mathrm{base}}. +\end{aligned} +``` + +Multiplying by $k_\mathrm{base}$ converts system base to component base. + +## Model Ports + +Name | Port | Init | Description +---------|--------|---------|------------ +`bus` | Bus | Known | Terminal-bus voltage +`pe` | Input | Known | Active-power feedback +`qgen` | Input | Known | Reactive-power feedback +`qext` | Input | Unknown | Volt/VAr reference +`pfaref` | Input | Unknown | Power-factor angle reference +`pref` | Input | Unknown | Active-power reference +`iqcmd` | Output | Known | Reactive-current command +`ipcmd` | Output | Known | Active-current command + +`bus` is required; signal ports are optional and must be linked when attached. +`Known` ports are seeded before `initialize()` and preserved by it. `Unknown` +inputs are resolved during initialization and written to attached signal +storage, or retained as constant inputs when the port is unattached. + +## Model Variables + +### Internal Variables + +#### Differential + +Symbol | Units | Description | Note +------------------------|--------|-------------------------------------|----- +$V^\mathrm{meas}$ | [p.u.] | Filtered terminal voltage | State 1 in Fig. 1 +$P^\mathrm{meas}$ | [p.u.] | Filtered electrical power | State 2 in Fig. 1; component base +$x_Q^\mathrm{PI}$ | [p.u.] | Reactive-power PI controller state | State 3 in Fig. 1 +$x_V^\mathrm{PI}$ | [p.u.] | Voltage-control PI controller state | State 4 in Fig. 1; component-base current +$Q_V$ | [p.u.] | Reactive-current command lag state | State 5 in Fig. 1; component base +$P^\mathrm{ord}$ | [p.u.] | Filtered active-power order | State 6 in Fig. 1; component base + +#### Algebraic + +Symbol | Units | Description | Note +---------------------|--------|-----------------------------------------------|----- +$V_T$ | [p.u.] | Terminal voltage magnitude | +$I_L^{\max}$ | [p.u.] | Current-circle continuation state | Component base +$I_q^\mathrm{cmd}$ | [p.u.] | Reactive-current command output | System base +$I_p^\mathrm{cmd}$ | [p.u.] | Active-current command output | System base + +### External Variables + +#### Differential + +None. + +#### Algebraic + +Symbol | Units | Init | 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 +$P_e$ | [p.u.] | Known | Electrical active-power feedback | Optional signal port `pe`; system base +$Q^\mathrm{gen}$ | [p.u.] | Known | Reactive-power feedback | Optional signal port `qgen`; system base +$Q^\mathrm{ext}$ | [p.u.] | Unknown | External Volt/VAr reference | Optional signal port `qext`; terminal voltage in direct-voltage mode +$\phi^\mathrm{ref}$ | [rad] | Unknown | Power-factor angle reference | Optional signal port `pfaref` +$P^\mathrm{ref}$ | [p.u.] | Unknown | External active-power reference | Optional signal port `pref`; system base + +## Model Equations + +For readability, define: + +```math +\begin{aligned} + V_\mathrm{safe}^\mathrm{meas} &= \text{max}(V^\mathrm{meas},0.01) \\ + s_\mathrm{dip} &= \text{inside}(V_T;\,V_\mathrm{dip},V_\mathrm{up}) \\ + e_V^\mathrm{db} &= \text{deadband2}(V^\mathrm{ref}-V^\mathrm{meas};\,D_1^\mathrm{db},D_2^\mathrm{db}) \\ + I_q^\mathrm{inj} &= \text{clamp}(K_\mathrm{qv}e_V^\mathrm{db};\,I_{q,\mathrm{inj}}^{\min},I_{q,\mathrm{inj}}^{\max}) \\ + Q^\mathrm{ref} &= s_Q^\mathrm{ref}(s_\mathrm{pf}P^\mathrm{meas}\tan(\phi^\mathrm{ref})+s_\mathrm{pf}^\mathrm{off}k_\mathrm{base}Q^\mathrm{ext}) \\ + e_Q &= \text{clamp}(Q^\mathrm{ref};\,Q^{\min},Q^{\max})-k_\mathrm{base}Q^\mathrm{gen} \\ + V_Q^\mathrm{PI} &= \text{clamp}(K_\mathrm{qp}e_Q+x_Q^\mathrm{PI};\,V^{\min},V^{\max}) \\ + e_V^\mathrm{PI} &= s_Q^\mathrm{PI}V_Q^\mathrm{PI}+s_V^\mathrm{ref}Q^\mathrm{ext}-s_QV^\mathrm{meas} \\ + f_P^\mathrm{ord} &= \dfrac{1}{T_\mathrm{pord}}(k_\mathrm{base}P^\mathrm{ref}-P^\mathrm{ord}) \\ + r_P^\mathrm{ord} &= \text{aslew}(f_P^\mathrm{ord};\,R_P^{\min},R_P^{\max}) \\ + N_L &= \sqrt{(I_L^{\max})^2+\epsilon_0},\qquad I_L^\mathrm{cap}=\dfrac{(I_L^{\max})^2}{N_L} \\ + I_q^{\max} &= s_\mathrm{pq}I_L^\mathrm{cap}+s_\mathrm{pq}^\mathrm{off}I^{\max} \\ + I_p^{\max} &= s_\mathrm{pq}I^{\max}+s_\mathrm{pq}^\mathrm{off}I_L^\mathrm{cap} \\ + I_q^\mathrm{base} &= \text{clamp}(K_\mathrm{vp}e_V^\mathrm{PI}+x_V^\mathrm{PI};\,-I_q^{\max},I_q^{\max}) \\ + I_q^\mathrm{raw} &= s_QI_q^\mathrm{base}+s_Q^\mathrm{off}Q_V+I_q^\mathrm{inj}. +\end{aligned} +``` + +CommonMath defines the [`antiwindup`](../../../../CommonMath.md#antiwindup) and +[smooth limiter](../../../../CommonMath.md#derived-functions) functions used in +these equations. [Appendix B](#appendix-b-aslew) defines `aslew`. + +### Differential Equations + +```math +\begin{aligned} + 0 &= -\dot{V}^\mathrm{meas} + \dfrac{1}{T_\mathrm{rv}}(V_T-V^\mathrm{meas}) \\ + 0 &= -\dot{P}^\mathrm{meas} + \dfrac{1}{T_\mathrm{p}}(k_\mathrm{base}P_e-P^\mathrm{meas}) \\ + 0 &= -\dot{x}_Q^\mathrm{PI} + s_Q^\mathrm{PI}s_\mathrm{dip}\,\text{antiwindup}(K_\mathrm{qp}e_Q+x_Q^\mathrm{PI},K_\mathrm{qi}e_Q;\,V^{\min},V^{\max}) \\ + 0 &= -\dot{x}_V^\mathrm{PI} + s_Qs_\mathrm{dip}\,\text{antiwindup}(K_\mathrm{vp}e_V^\mathrm{PI}+x_V^\mathrm{PI},K_\mathrm{vi}e_V^\mathrm{PI};\,-I_q^{\max},I_q^{\max}) \\ + 0 &= -\dot{Q}_V + \dfrac{1}{T_\mathrm{iq}}s_Q^\mathrm{off}s_\mathrm{dip}\left(\dfrac{Q^\mathrm{ref}}{V_\mathrm{safe}^\mathrm{meas}}-Q_V\right) \\ + 0 &= -\dot{P}^\mathrm{ord} + s_\mathrm{dip}\,\text{antiwindup}(P^\mathrm{ord},r_P^\mathrm{ord};\,P^{\min},P^{\max}). +\end{aligned} +``` + +### Algebraic Equations + +```math +\begin{aligned} + 0 &= -V_T^2+V_\mathrm{r}^2+V_\mathrm{i}^2 \\ + 0 &= -I_L^{\max}N_L+(I^{\max})^2-s_\mathrm{pq}(k_\mathrm{base}I_p^\mathrm{cmd})^2-s_\mathrm{pq}^\mathrm{off}(k_\mathrm{base}I_q^\mathrm{cmd})^2 \\ + 0 &= -k_\mathrm{base}I_q^\mathrm{cmd}+\text{clamp}(I_q^\mathrm{raw};\,-I_q^{\max},I_q^{\max}) \\ + 0 &= -k_\mathrm{base}I_p^\mathrm{cmd}+\text{clamp}\left(\dfrac{P^\mathrm{ord}}{V_\mathrm{safe}^\mathrm{meas}};\,0,I_p^{\max}\right). +\end{aligned} +``` + +Here $\epsilon_0=100\epsilon_\mathrm{machine}$ regularizes the `ILMAX` row at +zero remaining capacity. + +## Initialization + +REECB reconstructs a steady operating point. Arbitrary-state restart is unsupported. + +### Input Initialization + +```math +\begin{aligned} + V_\mathrm{r},V_\mathrm{i} &\leftarrow \text{terminal-bus voltage} \\ + I_q^\mathrm{cmd},I_p^\mathrm{cmd} &\leftarrow \text{owned current-command variables} \\ + P_e &\leftarrow \text{attached active-power feedback},\quad \text{if attached} \\ + Q^\mathrm{gen} &\leftarrow \text{attached reactive-power feedback},\quad \text{if attached}. +\end{aligned} +``` + +### Internal Initialization + +Initialization resolves the steady-state quantities in dependency order; all +internal derivatives start at zero. Let $I_p=k_\mathrm{base}I_p^\mathrm{cmd}$ +and $I_q=k_\mathrm{base}I_q^\mathrm{cmd}$ be the component-base commands. +[Appendix A](#appendix-a-iclamp) defines the initialization-only `iclamp`. + +```math +\begin{aligned} + V_T &\leftarrow \sqrt{V_\mathrm{r}^2+V_\mathrm{i}^2} \\ + V^\mathrm{ref} &\leftarrow V_T,\quad \text{if omitted} \\ + V^\mathrm{meas} &\leftarrow V_T \\ + V_\mathrm{safe}^\mathrm{meas} &\leftarrow \text{max}(V^\mathrm{meas},0.01) \\ + P_e &\leftarrow V_\mathrm{safe}^\mathrm{meas}I_p^\mathrm{cmd},\quad \text{if unattached} \\ + Q^\mathrm{gen} &\leftarrow V_\mathrm{safe}^\mathrm{meas}I_q^\mathrm{cmd},\quad \text{if unattached} \\ + P^\mathrm{meas} &\leftarrow k_\mathrm{base}P_e \\ + e_V^\mathrm{db} &\leftarrow \text{deadband2}(V^\mathrm{ref}-V^\mathrm{meas};\,D_1^\mathrm{db},D_2^\mathrm{db}) \\ + I_q^\mathrm{inj} &\leftarrow \text{clamp}(K_\mathrm{qv}e_V^\mathrm{db};\,I_{q,\mathrm{inj}}^{\min},I_{q,\mathrm{inj}}^{\max}) \\ + I_L^{\max}N_L &\leftarrow (I^{\max})^2-s_\mathrm{pq}I_p^2-s_\mathrm{pq}^\mathrm{off}I_q^2,\qquad I_L^{\max}\ge0 \\ + I_q^{\max} &\leftarrow s_\mathrm{pq}I_L^\mathrm{cap}+s_\mathrm{pq}^\mathrm{off}I^{\max},\qquad I_p^{\max}\leftarrow s_\mathrm{pq}I^{\max}+s_\mathrm{pq}^\mathrm{off}I_L^\mathrm{cap}. +\end{aligned} +``` + +Initialization raises $I^\max$ when needed to reproduce the supplied current +commands; incompatible reactive-current injection is rejected. Q, V, and P +limits are expanded as needed to include their initialized values, and each +adjustment logs a warning. + +```math +\begin{aligned} + I_q^\mathrm{raw} &\leftarrow \text{iclamp}(I_q;\,-I_q^{\max},I_q^{\max}) \\ + I_q^\mathrm{ctrl} &\leftarrow I_q^\mathrm{raw}-I_q^\mathrm{inj} \\ + P^\mathrm{ord} &\leftarrow V_\mathrm{safe}^\mathrm{meas}\text{iclamp}(I_p;\,0,I_p^{\max}) \\ + Q^\mathrm{target} &\leftarrow + \begin{cases} + V_\mathrm{safe}^\mathrm{meas}I_q^\mathrm{ctrl} & s_Q=0 \\ + \text{iclamp}(k_\mathrm{base}Q^\mathrm{gen};\,Q^{\min},Q^{\max}) & s_Qs_V=1 \\ + 0 & \text{otherwise} + \end{cases}. +\end{aligned} +``` + +```math +\begin{aligned} + \phi^\mathrm{ref} &\leftarrow + \begin{cases} + \arctan(Q^\mathrm{target}/P^\mathrm{meas}) & s_\mathrm{pf}=1\ \land\ P^\mathrm{meas}\ne0 \\ + 0 & \text{otherwise} + \end{cases} \\ + Q^\mathrm{ext} &\leftarrow + \begin{cases} + V^\mathrm{meas} & s_V^\mathrm{ref}=1 \\ + 0 & s_V^\mathrm{ref}=0\ \land\ s_\mathrm{pf}=1 \\ + Q^\mathrm{target}/k_\mathrm{base} & s_V^\mathrm{ref}=0\ \land\ s_\mathrm{pf}=0 + \end{cases} \\ + Q^\mathrm{ref} &\leftarrow s_Q^\mathrm{ref}(s_\mathrm{pf}P^\mathrm{meas}\tan(\phi^\mathrm{ref})+s_\mathrm{pf}^\mathrm{off}k_\mathrm{base}Q^\mathrm{ext}). +\end{aligned} +``` + +```math +\begin{aligned} + e_Q &\leftarrow \text{clamp}(Q^\mathrm{ref};\,Q^{\min},Q^{\max})-k_\mathrm{base}Q^\mathrm{gen} \\ + x_Q^\mathrm{PI} &\leftarrow + \begin{cases} + \text{iclamp}(V^\mathrm{meas};\,V^{\min},V^{\max})-K_\mathrm{qp}e_Q & s_Qs_V=1 \\ + 0 & s_Qs_V=0 + \end{cases} \\ + V_Q^\mathrm{PI} &\leftarrow \text{clamp}(K_\mathrm{qp}e_Q+x_Q^\mathrm{PI};\,V^{\min},V^{\max}) \\ + e_V^\mathrm{PI} &\leftarrow s_Q^\mathrm{PI}V_Q^\mathrm{PI}+s_V^\mathrm{ref}Q^\mathrm{ext}-s_QV^\mathrm{meas} \\ + x_V^\mathrm{PI} &\leftarrow + \begin{cases} + -K_\mathrm{vp}e_V^\mathrm{PI} & s_Q=1\ \land\ I_q^{\max}\le\epsilon_0 \\ + \text{iclamp}(I_q^\mathrm{ctrl};\,-I_q^{\max},I_q^{\max})-K_\mathrm{vp}e_V^\mathrm{PI} & s_Q=1\ \land\ I_q^{\max}>\epsilon_0 \\ + 0 & s_Q=0 + \end{cases} \\ + Q_V &\leftarrow + \begin{cases} + 0 & s_Q=1 \\ + Q^\mathrm{ref}/V_\mathrm{safe}^\mathrm{meas} & s_Q=0 + \end{cases}. +\end{aligned} +``` + +Invalid or non-equilibrium operating points are rejected before any state, +limit, latch, or signal is changed. + +### Output Initialization + +```math +\begin{aligned} + \phi^\mathrm{ref} &\leftarrow + \begin{cases} + \arctan(Q^\mathrm{target}/P^\mathrm{meas}) & s_\mathrm{pf}=1\ \land\ P^\mathrm{meas}\ne0 \\ + 0 & \text{otherwise} + \end{cases} \\ + Q^\mathrm{ext} &\leftarrow + \begin{cases} + V^\mathrm{meas} & s_V^\mathrm{ref}=1 \\ + 0 & s_V^\mathrm{ref}=0\ \land\ s_\mathrm{pf}=1 \\ + Q^\mathrm{target}/k_\mathrm{base} & s_V^\mathrm{ref}=0\ \land\ s_\mathrm{pf}=0 + \end{cases} \\ + P^\mathrm{ref} &\leftarrow \dfrac{P^\mathrm{ord}}{k_\mathrm{base}} +\end{aligned} +``` + +REECB writes the resolved references to attached signal inputs; unattached +ports retain them as constant inputs. + +## Monitorable Outputs + +Output | Units | Description | Note +--------|--------|---------------------------------|----- +`iqcmd` | [p.u.] | Reactive-current command output | $I_q^\mathrm{cmd}$ (system base) +`ipcmd` | [p.u.] | Active-current command output | $I_p^\mathrm{cmd}$ (system base) +`vmeas` | [p.u.] | Filtered terminal voltage | $V^\mathrm{meas}$ +`pmeas` | [p.u.] | Filtered electrical power | $P^\mathrm{meas}$ (component base) + +## Testing + +- `validation()` checks configuration and defaults. +- `initializationAndSignals()` checks initialization, signals, monitors, and power bases. +- `initializationDomain()` checks rejected inputs and limit expansion. +- `initializationExactness()` checks endpoint and current-circle initialization. +- `residualEquations()` checks the fixed residual answer key. +- `selectorConfigurations()` checks selectors and optional ports. +- `voltVarReferenceBase()` checks `qext` units. +- `reactiveControl()` checks the reactive-control paths. +- `activeCurrentControl()` checks active-current control and current priority. +- `dependencyTracking()` checks sparse dependencies. +- `jacobian()` compares the Enzyme and dependency-tracking Jacobians. +- `regcaReecb()` checks REGCA-REECB signal wiring. +- `reecb()` checks construction through the production system-data path. + +## Appendix A: `iclamp` + +For $\ell + int Reecb::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Reecb...\n"; + Log::misc() << "Jacobian evaluation is not implemented!\n"; + return 0; + } + + template class Reecb; + template class Reecb; + } // namespace Controller + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Controller/REECB/Reecb.hpp b/GridKit/Model/PhasorDynamics/Controller/REECB/Reecb.hpp new file mode 100644 index 000000000..7e000aac0 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Controller/REECB/Reecb.hpp @@ -0,0 +1,241 @@ +/** + * @file Reecb.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Declaration of the REECB electrical-control model. + */ + +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + template + class BusBase; + + namespace Controller + { + /// Internal variables and residual rows of a `Reecb`. + enum class ReecbInternalVariables : size_t + { + VMEAS, ///< \f$V^\mathrm{meas}\f$ Differential filtered terminal voltage [p.u.] + PMEAS, ///< \f$P^\mathrm{meas}\f$ Differential filtered electrical power on component base [p.u.] + XPIQ, ///< \f$x_Q^\mathrm{PI}\f$ Differential reactive-power PI state [p.u.] + XPIV, ///< \f$x_V^\mathrm{PI}\f$ Differential voltage-control PI state on component base [p.u.] + QV, ///< \f$Q_V\f$ Differential reactive-current command lag state on component base [p.u.] + PORD, ///< \f$P^\mathrm{ord}\f$ Differential filtered active-power order on component base [p.u.] + VT, ///< \f$V_T\f$ Algebraic terminal-voltage magnitude [p.u.] + ILMAX, ///< \f$I_L^\max\f$ Algebraic current-circle continuation state on component base [p.u.] + IQCMD, ///< \f$I_q^\mathrm{cmd}\f$ Algebraic reactive-current command output on system base [p.u.] + IPCMD, ///< \f$I_p^\mathrm{cmd}\f$ Algebraic active-current command output on system base [p.u.] + MAXIMUM ///< Number of REECB internal variables and residual rows + }; + + /// External signal variables read or initialized by a `Reecb`. + enum class ReecbExternalVariables : size_t + { + PE, ///< \f$P_e\f$ Optional Known active-power feedback input on system base [p.u.] + QGEN, ///< \f$Q^\mathrm{gen}\f$ Optional Known reactive-power feedback input on system base [p.u.] + QEXT, ///< \f$Q^\mathrm{ext}\f$ Optional Unknown Volt/VAr reference input: system-base reactive power [p.u.], or the terminal-voltage reference [p.u.] when \f$s_Q=1\f$ and \f$s_V=0\f$ + PFAREF, ///< \f$\phi^\mathrm{ref}\f$ Optional Unknown power-factor angle-reference input [rad] + PREF, ///< \f$P^\mathrm{ref}\f$ Optional Unknown active-power reference input on system base [p.u.] + MAXIMUM ///< Number of REECB external signal variables + }; + + /** + * @brief WECC renewable electrical controller with reactive-power, + * voltage, and active-power command paths (REECB). + * + * @tparam scalar_type Plain real or differentiable scalar type. + * @tparam index_type Integer index type. + */ + template + class Reecb : public Component + { + using Component::gridkit_component_id_; + using Component::alpha_; + using Component::allocated_; + using Component::abs_tol_; + using Component::f_; + 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::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 ModelDataT = ReecbData; + using MonitorT = Model::VariableMonitor; + using InternalVariablesT = ReecbInternalVariables; + using ExternalVariablesT = ReecbExternalVariables; + + /// Current-circle regularization and initialization reconstruction tolerance. + static constexpr RealT INITIALIZATION_TOLERANCE = + static_cast(100.0) * std::numeric_limits::epsilon(); + + Reecb(BusT* bus); + Reecb(BusT* bus, const ModelDataT& data); + ~Reecb(); + + int setGridKitComponentID(IdxT component_id) override final; + int allocate() override final; + int verify() const override final; + int initialize() override final; + int tagDifferentiable() override final; + int setAbsoluteTolerance(RealT rel_tol) override final; + int evaluateResidual() override final; + int evaluateJacobian() override final; + + auto getSignals() + -> ComponentSignals&; + + const Model::VariableMonitorBase* getMonitor() const override; + + [[gnu::always_inline]] inline int evaluateInternalResidual( + const ScalarT* y, + const ScalarT* yp, + const ScalarT* wb, + const ScalarT* ws, + ScalarT* f); + + private: + /// Smooth asymmetric slew-rate limiter. + [[gnu::always_inline]] static inline ScalarT aslew(ScalarT rate, RealT lower, RealT upper); + + /// Smooth anti-windup derivative within a moving symmetric band. + [[gnu::always_inline]] static inline ScalarT awband(ScalarT state, ScalarT rate, ScalarT band); + + /// Current-circle continuation state for an initial component-base limit. + static RealT circleState(RealT imax, RealT high); + + /// Off-axis component-base capacity provided by a continuation state. + static RealT capacity(RealT ilmax); + + /// Bisect an initial-limit bracket to its first upper-side point. + template + static RealT bisect(RealT a, RealT b, FuncT below); + + /// Solve the smallest feasible initial limit at or above `lower`. + static RealT solveInitialLimit(RealT lower, RealT high, RealT low); + + void loadRealParameter(const ModelDataT& data, + ReecbParameters parameter, + RealT& target, + const char* name); + void loadBooleanParameter(const ModelDataT& data, + ReecbParameters parameter, + bool& target, + const char* name); + bool floorTimeConstant(RealT& value, const char* name); + void initializeParameters(const ModelDataT& data); + void initializeMonitor(); + void setDerivedParameters(); + + static RealT logOneMinusExp(RealT x); + bool iclamp(RealT output, RealT lower, RealT upper, RealT& input) const; + RealT componentPowerBase() const; + + template + [[gnu::always_inline]] inline ValueT toComponentBase(ValueT value) const; + + template + ValueT toSystemBase(ValueT value) const; + + ScalarT& Vr(); + ScalarT& Vi(); + + static constexpr RealT TIME_CONSTANT_MINIMUM = static_cast(1.0e-3); + static constexpr RealT VMEAS_MINIMUM = static_cast(0.01); + + BusT* bus_{nullptr}; + + // Input parameters + RealT mva_base_{0}; + bool PfFlag_{false}; + bool VFlag_{false}; + bool QFlag_{false}; + bool Pqflag_{false}; + RealT Trv_{0.02}; + RealT Tp_{0}; + RealT Vref0_{0}; + RealT Vdip_{0.85}; + RealT Vup_{1.15}; + RealT dbd1_{0}; + RealT dbd2_{0}; + RealT kqv_{5.0}; + RealT Iql1_{-1.1}; + RealT Iqh1_{1.1}; + RealT Qmax_{0.436}; + RealT Qmin_{-0.436}; + RealT Kqp_{0}; + RealT Kqi_{0.1}; + RealT Vmax_{1.1}; + RealT Vmin_{0.9}; + RealT Kvp_{18.0}; + RealT Kvi_{5.0}; + RealT Tiq_{0.02}; + RealT Tpord_{0.02}; + RealT dPmax_{99.0}; + RealT dPmin_{-99.0}; + RealT Pmax_{1}; + RealT Pmin_{0}; + RealT Imax_{1.3}; + + bool mva_given_{false}; + bool Vref0_given_{false}; + IdxT parameter_error_count_{0}; + + // Derived parameters + RealT va_component_base_{0}; + RealT pf_on_{0}; + RealT pf_off_{1}; + RealT q_on_{0}; + RealT q_off_{1}; + RealT q_pi_on_{0}; + RealT v_ref_on_{0}; + RealT q_ref_on_{1}; + RealT pq_on_{0}; + RealT pq_off_{1}; + + // Unattached signal setpoints + ScalarT pe_set_{0}; + ScalarT qgen_set_{0}; + ScalarT qext_set_{0}; + ScalarT pfaref_set_{0}; + ScalarT pref_set_{0}; + + ComponentSignals signals_; + std::unique_ptr monitor_; + + // Local copies of signal variables + std::vector ws_; + std::vector ws_indices_; + }; + } // namespace Controller + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbData.hpp b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbData.hpp new file mode 100644 index 000000000..f14613d23 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbData.hpp @@ -0,0 +1,114 @@ +/** + * @file ReecbData.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Modeling data for the REECB electrical-control model. + */ + +#pragma once + +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Controller + { + /// Parameters for REECB. + enum class ReecbParameters + { + mva, ///< \f$S^\mathrm{base}\f$ Component power base [MVA] + PfFlag, ///< \f$s_\mathrm{pf}\f$ Power-factor control selector: true = power-factor control, false = reactive-power control [boolean] + VFlag, ///< \f$s_V\f$ Voltage-reference selector under \f$s_Q=1\f$: true = cascaded Q-PI voltage command, false = direct external voltage reference [boolean] + QFlag, ///< \f$s_Q\f$ Reactive-path selector: true = Volt/VAr PI control, false = reactive-current lag [boolean] + Pqflag, ///< \f$s_\mathrm{pq}\f$ Converter current-priority selector: true = P priority, false = Q priority [boolean] + Trv, ///< \f$T_\mathrm{rv}\f$ Voltage-measurement filter time constant [sec] + Tp, ///< \f$T_\mathrm{p}\f$ Electrical-power measurement filter time constant [sec] + Vref0, ///< \f$V^\mathrm{ref}\f$ Reactive-current-injection voltage reference [p.u.] + Vdip, ///< \f$V_\mathrm{dip}\f$ Low-voltage threshold for the voltage-band gate [p.u.] + Vup, ///< \f$V_\mathrm{up}\f$ High-voltage threshold for the voltage-band gate [p.u.] + dbd1, ///< \f$D_1^\mathrm{db}\f$ Lower voltage-error deadband threshold [p.u.] + dbd2, ///< \f$D_2^\mathrm{db}\f$ Upper voltage-error deadband threshold [p.u.] + kqv, ///< \f$K_\mathrm{qv}\f$ Reactive-current injection gain [p.u.] + Iql1, ///< \f$I_{q,\mathrm{inj}}^\min\f$ Minimum reactive-current injection on component base [p.u.] + Iqh1, ///< \f$I_{q,\mathrm{inj}}^\max\f$ Maximum reactive-current injection on component base [p.u.] + Qmax, ///< \f$Q^\max\f$ Maximum reactive-power control output on component base [p.u.] + Qmin, ///< \f$Q^\min\f$ Minimum reactive-power control output on component base [p.u.] + Kqp, ///< \f$K_\mathrm{qp}\f$ Reactive-power proportional gain [p.u.] + Kqi, ///< \f$K_\mathrm{qi}\f$ Reactive-power integral gain [p.u./s] + Vmax, ///< \f$V^\max\f$ Maximum voltage-control output [p.u.] + Vmin, ///< \f$V^\min\f$ Minimum voltage-control output [p.u.] + Kvp, ///< \f$K_\mathrm{vp}\f$ Voltage-control proportional gain [p.u.] + Kvi, ///< \f$K_\mathrm{vi}\f$ Voltage-control integral gain [p.u./s] + Tiq, ///< \f$T_\mathrm{iq}\f$ Reactive-current command lag time constant [sec] + Tpord, ///< \f$T_\mathrm{pord}\f$ Active-power order filter time constant [sec] + dPmax, ///< \f$R_P^\max\f$ Positive active-power ramp-rate limit on component base [p.u./s] + dPmin, ///< \f$R_P^\min\f$ Negative active-power ramp-rate limit on component base [p.u./s] + Pmax, ///< \f$P^\max\f$ Maximum active-power order limit on component base [p.u.] + Pmin, ///< \f$P^\min\f$ Minimum active-power order limit on component base [p.u.] + Imax ///< \f$I^\max\f$ Maximum converter current on component base [p.u.] + }; + + /// Buses for the REECB electrical-control model. + enum class ReecbBuses : size_t + { + bus, ///< \f$V_\mathrm{r},V_\mathrm{i}\f$ Required Known terminal-bus voltage [p.u.] + SIZE ///< Number of REECB bus ports + }; + + /// Signal inputs for the REECB electrical-control model. + enum class ReecbSignalInputs : size_t + { + pe, ///< \f$P_e\f$ Optional Known active-power feedback input on system base [p.u.] + qgen, ///< \f$Q^\mathrm{gen}\f$ Optional Known reactive-power feedback input on system base [p.u.] + qext, ///< \f$Q^\mathrm{ext}\f$ Optional Unknown Volt/VAr reference input: system-base reactive power [p.u.], or the terminal-voltage reference [p.u.] when \f$s_Q=1\f$ and \f$s_V=0\f$ + pfaref, ///< \f$\phi^\mathrm{ref}\f$ Optional Unknown power-factor angle-reference input [rad] + pref, ///< \f$P^\mathrm{ref}\f$ Optional Unknown active-power reference input on system base [p.u.] + SIZE ///< Number of REECB signal-input ports + }; + + /// Signal outputs for the REECB electrical-control model. + enum class ReecbSignalOutputs : size_t + { + iqcmd, ///< \f$I_q^\mathrm{cmd}\f$ Optional Known reactive-current command output on system base [p.u.] + ipcmd, ///< \f$I_p^\mathrm{cmd}\f$ Optional Known active-current command output on system base [p.u.] + SIZE ///< Number of REECB signal-output ports + }; + + /// Variables available through the monitor interface. + enum class ReecbMonitorableVariables + { + iqcmd, ///< \f$I_q^\mathrm{cmd}\f$ Reactive-current command output on system base [p.u.] + ipcmd, ///< \f$I_p^\mathrm{cmd}\f$ Active-current command output on system base [p.u.] + vmeas, ///< \f$V^\mathrm{meas}\f$ Filtered terminal voltage [p.u.] + pmeas ///< \f$P^\mathrm{meas}\f$ Filtered electrical power on component base [p.u.] + }; + + /** + * @brief Model data for REECB parameters, bus and signal ports, and monitored variables. + * + * @tparam real_type Real parameter value type. + * @tparam index_type Integer index type. + * + * @see Reecb + */ + template + struct ReecbData : public ComponentData + { + ReecbData() = default; + + using Parameters = ReecbParameters; + using Buses = ReecbBuses; + using SignalInputs = ReecbSignalInputs; + using SignalOutputs = ReecbSignalOutputs; + using MonitorableVariables = ReecbMonitorableVariables; + }; + } // namespace Controller + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbDependencyTracking.cpp new file mode 100644 index 000000000..1b2f97f32 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbDependencyTracking.cpp @@ -0,0 +1,31 @@ +/** + * @file ReecbDependencyTracking.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Dependency-tracking instantiations for the REECB electrical-control model. + */ + +#include "ReecbImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Controller + { + /** + * @brief Report that DependencyTracking exposes structure through the + * residual rather than a separately assembled Jacobian. + */ + template + int Reecb::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Reecb...\n"; + Log::misc() << "Jacobian evaluation is not implemented!\n"; + return 0; + } + + template class Reecb; + template class Reecb; + } // namespace Controller + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbEnzyme.cpp b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbEnzyme.cpp new file mode 100644 index 000000000..c71231ee8 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbEnzyme.cpp @@ -0,0 +1,122 @@ +/** + * @file ReecbEnzyme.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Enzyme sparse Jacobian for the REECB electrical-control model. + */ + +#include + +#include "ReecbImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Controller + { + /** + * @brief Assemble the sparse REECB Jacobian with Enzyme. + * + * Differentiates the internal residual with respect to state, derivative, + * terminal-bus, and linked signal variables, then constructs the model + * COO matrix. + * + * @pre allocate() has sized the model and Jacobian index maps. + * @pre evaluateResidual() has refreshed the current bus/signal values and + * signal indices. + * @pre The containing solver has set the current integration coefficient + * and global variable/residual indices. + */ + template + int Reecb::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Reecb...\n"; + Log::misc() << "Jacobian evaluation is experimental!\n"; + + if (J_rows_buffer_ == nullptr) + { + const auto size = static_cast(size_); + const auto bus_size = static_cast(bus_->size()); + const auto signal_size = ws_.size(); + const auto buffer_size = 2 * size * size + 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 ModelT = GridKit::PhasorDynamics::Controller::Reecb; + using Fn = GridKit::Enzyme::Sparse::MemberFunctions; + + nnz_ = 0; + + GridKit::Enzyme::Sparse::DfDy::eval( + this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + this->getResidualIndices().data(), + this->getVariableIndices().data(), + y_.getData(), + yp_.getData(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + GridKit::Enzyme::Sparse::DfDyp::eval( + this, + static_cast(f_.getSize()), + static_cast(y_.getSize()), + this->getResidualIndices().data(), + this->getVariableIndices().data(), + y_.getData(), + yp_.getData(), + wb_.data(), + ws_.data(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + GridKit::Enzyme::Sparse::DfDwb::eval( + this, + static_cast(f_.getSize()), + static_cast(bus_->size()), + this->getResidualIndices().data(), + bus_->getVariableIndices().data(), + y_.getData(), + yp_.getData(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + GridKit::Enzyme::Sparse::DfDws::eval( + this, + static_cast(f_.getSize()), + ws_.size(), + this->getResidualIndices().data(), + ws_indices_.data(), + y_.getData(), + yp_.getData(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + this->constructCoo(); + + return 0; + } + + template class Reecb; + template class Reecb; + } // namespace Controller + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp new file mode 100644 index 000000000..edd7daf35 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Controller/REECB/ReecbImpl.hpp @@ -0,0 +1,1468 @@ +/** + * @file ReecbImpl.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Definition of the REECB electrical-control model. + */ + +#pragma once + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Controller + { + /// Logger used for REECB diagnostics. + using Log = ::GridKit::Utilities::Logger; + + /** + * @brief Construct REECB with its documented parameter defaults + * + * The terminal bus is retained, the model is sized, and no monitor or + * signal connection is created. + * + * @param[in] bus Terminal bus measured by the controller. + */ + template + Reecb::Reecb(BusT* bus) + : bus_(bus) + { + size_ = static_cast(ReecbInternalVariables::MAXIMUM); + setDerivedParameters(); + } + + /** + * @brief Construct REECB from model data + * + * @param[in] bus Terminal bus measured by the controller. + * @param[in] data Model parameters and monitor selections. + */ + template + Reecb::Reecb(BusT* bus, const ModelDataT& data) + : bus_(bus), + monitor_(std::make_unique(data)) + { + initializeParameters(data); + initializeMonitor(); + size_ = static_cast(ReecbInternalVariables::MAXIMUM); + } + + /** + * @brief Destroy the electrical controller and its optional variable monitor. + */ + template + Reecb::~Reecb() + { + } + + /** + * @brief Set the component ID + * + * @param[in] component_id Identifier assigned by the system model. + */ + template + int Reecb::setGridKitComponentID(IdxT component_id) + { + gridkit_component_id_ = component_id; + return 0; + } + + /** + * @brief Allocate model vectors and wire assigned current-command outputs + * + * Sizes the state, residual, bus, and signal-interface buffers, initializes + * identity index maps, and points assigned command nodes at the internal + * system-base states that REECB publishes. Repeated allocation reuses the + * existing model vectors and signal links. + */ + template + int Reecb::allocate() + { + const auto IQCMD = static_cast(ReecbInternalVariables::IQCMD); + const auto IPCMD = static_cast(ReecbInternalVariables::IPCMD); + + if (!allocated_) + { + this->allocateVectors(size_); + } + const auto size = static_cast(size_); + + tag_.assign(size, false); + variable_indices_.resize(size); + residual_indices_.resize(size); + + wb_.assign(2, ScalarT{0}); + + const auto signal_size = static_cast(ReecbExternalVariables::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); + } + + auto* y = y_.getData(); + + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y[IQCMD], + &(this->getVariableIndex(static_cast(IQCMD)))); + } + + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y[IPCMD], + &(this->getVariableIndex(static_cast(IPCMD)))); + } + + allocated_ = true; + return 0; + } + + /** + * @brief Validate the REECB configuration + * + * Checks parameter-loading errors, finiteness and static relationships, + * system/component bases and conversion ratios, the terminal bus, and + * attached optional signals. Operating-point feasibility is checked by + * initialize(). + * + * @return Number of configuration errors; zero when valid. + */ + template + int Reecb::verify() const + { + int ret = static_cast(parameter_error_count_); + + auto check = [&](bool condition, const char* message) + { + if (!condition) + { + Log::error() << "Reecb: " << message << '\n'; + ret += 1; + } + }; + + check(bus_ != nullptr, "terminal bus is required"); + + const RealT component_power_base = componentPowerBase(); + const bool valid_component_base = std::isfinite(component_power_base) && component_power_base > ZERO; + const bool valid_system_base = std::isfinite(va_system_base_) && va_system_base_ > ZERO; + check(valid_component_base, "component power base must be finite and positive"); + check(valid_system_base, "system power base must be finite and positive"); + if (valid_component_base && valid_system_base) + { + const RealT system_to_component = va_system_base_ / component_power_base; + const RealT component_to_system = component_power_base / va_system_base_; + check( + std::isfinite(system_to_component) + && system_to_component > ZERO + && std::isfinite(component_to_system) + && component_to_system > ZERO, + "system/component power-base conversion ratios must be finite and positive"); + } + + check(std::isfinite(Trv_), "Trv must be finite"); + check(std::isfinite(Tp_), "Tp must be finite"); + check(std::isfinite(Vref0_), "Vref0 must be finite"); + + const bool finite_voltage_thresholds = std::isfinite(Vdip_) && std::isfinite(Vup_); + check(finite_voltage_thresholds, "Vdip and Vup must be finite"); + if (finite_voltage_thresholds) + { + check(Vdip_ < Vup_, "Vdip must be less than Vup"); + } + + const bool finite_voltage_deadband = std::isfinite(dbd1_) && std::isfinite(dbd2_); + check(finite_voltage_deadband, "dbd1 and dbd2 must be finite"); + if (finite_voltage_deadband) + { + check(dbd1_ <= ZERO && ZERO <= dbd2_, "dbd1 <= 0 <= dbd2 is required"); + } + + check(std::isfinite(kqv_) && kqv_ >= ZERO, "kqv must be finite and non-negative"); + + const bool finite_injection_limits = std::isfinite(Iql1_) && std::isfinite(Iqh1_); + check(finite_injection_limits, "Iql1 and Iqh1 must be finite"); + if (finite_injection_limits) + { + check(Iql1_ <= Iqh1_, "Iql1 must be less than or equal to Iqh1"); + } + + const bool finite_reactive_limits = std::isfinite(Qmin_) && std::isfinite(Qmax_); + check(finite_reactive_limits, "Qmin and Qmax must be finite"); + if (finite_reactive_limits) + { + check(Qmin_ <= Qmax_, "Qmin must be less than or equal to Qmax"); + } + + check(std::isfinite(Kqp_) && Kqp_ >= ZERO, "Kqp must be finite and non-negative"); + check(std::isfinite(Kqi_) && Kqi_ >= ZERO, "Kqi must be finite and non-negative"); + + const bool finite_voltage_limits = std::isfinite(Vmin_) && std::isfinite(Vmax_); + check(finite_voltage_limits, "Vmin and Vmax must be finite"); + if (finite_voltage_limits) + { + check(Vmin_ <= Vmax_, "Vmin must be less than or equal to Vmax"); + } + + check(std::isfinite(Kvp_) && Kvp_ >= ZERO, "Kvp must be finite and non-negative"); + check(std::isfinite(Kvi_) && Kvi_ >= ZERO, "Kvi must be finite and non-negative"); + check(std::isfinite(Tiq_), "Tiq must be finite"); + check(std::isfinite(Tpord_), "Tpord must be finite"); + + const bool finite_ramp_limits = std::isfinite(dPmin_) && std::isfinite(dPmax_); + check(finite_ramp_limits, "dPmin and dPmax must be finite"); + if (finite_ramp_limits) + { + check(dPmin_ < ZERO && ZERO < dPmax_, "dPmin < 0 < dPmax is required"); + } + + const bool finite_active_limits = std::isfinite(Pmin_) && std::isfinite(Pmax_); + check(finite_active_limits, "Pmin and Pmax must be finite"); + if (finite_active_limits) + { + check(Pmin_ <= Pmax_, "Pmin must be less than or equal to Pmax"); + } + + check(std::isfinite(Imax_) && Imax_ > ZERO, "Imax must be finite and positive"); + + auto check_optional_signal = [&](const char* name) + { + if (signals_.template isAttached() && !signals_.template isLinked()) + { + Log::error() << "Reecb: " << name << " signal attached with no linked source\n"; + ret += 1; + } + }; + + check_optional_signal.template operator()("pe"); + check_optional_signal.template operator()("qgen"); + check_optional_signal.template operator()("qext"); + check_optional_signal.template operator()("pfaref"); + check_optional_signal.template operator()("pref"); + + return ret; + } + + /** + * @brief Initialize REECB from the initial current commands and feedback + * + * Preserves the system-base command states, consumes attached initialized + * power feedback or reconstructs unattached feedback, and constructs the + * remaining states and reference setpoints at a steady operating point. + * + * @pre allocate() has completed. + * @pre verify() reports a valid parameter and port configuration. + * @pre The terminal bus and current-command states are initialized. + * + * @post On failure no state, derivative, parameter, or signal storage is + * modified. + * + * @return 0 on success; nonzero when allocation, configuration, initial-value, limiter inversion, or steady-state checks fail. + */ + template + int Reecb::initialize() + { + const auto VMEAS = static_cast(ReecbInternalVariables::VMEAS); + const auto PMEAS = static_cast(ReecbInternalVariables::PMEAS); + const auto XPIQ = static_cast(ReecbInternalVariables::XPIQ); + const auto XPIV = static_cast(ReecbInternalVariables::XPIV); + const auto QV = static_cast(ReecbInternalVariables::QV); + const auto PORD = static_cast(ReecbInternalVariables::PORD); + const auto VT = static_cast(ReecbInternalVariables::VT); + const auto ILMAX = static_cast(ReecbInternalVariables::ILMAX); + const auto IQCMD = static_cast(ReecbInternalVariables::IQCMD); + const auto IPCMD = static_cast(ReecbInternalVariables::IPCMD); + + if (!allocated_) + { + Log::error() << "Reecb: allocate must complete before initialize\n"; + return 1; + } + + if (verify() > 0) + { + Log::error() << "Reecb: cannot initialize with invalid configuration\n"; + return 1; + } + + auto* y = y_.getData(); + + const RealT ipcmd0_system = static_cast(y[IPCMD]); + const RealT iqcmd0_system = static_cast(y[IQCMD]); + const RealT ipcmd0 = toComponentBase(ipcmd0_system); + const RealT iqcmd0 = toComponentBase(iqcmd0_system); + const RealT vr0 = static_cast(Vr()); + const RealT vi0 = static_cast(Vi()); + const RealT vt0 = std::sqrt(vr0 * vr0 + vi0 * vi0); + const RealT vmeas0 = vt0; + const RealT vmeas_safe0 = Math::max(vmeas0, VMEAS_MINIMUM); + + RealT pe0_system = toSystemBase(ipcmd0 * vmeas_safe0); + RealT qgen0_system = toSystemBase(iqcmd0 * vmeas_safe0); + + if (signals_.template isAttached()) + { + pe0_system = static_cast( + signals_.template readExternalVariable()); + } + if (signals_.template isAttached()) + { + qgen0_system = static_cast( + signals_.template readExternalVariable()); + } + + const RealT pmeas0 = toComponentBase(pe0_system); + const RealT qgen0 = toComponentBase(qgen0_system); + RealT vref0 = vmeas0; + if (Vref0_given_) + { + vref0 = Vref0_; + } + + if (!std::isfinite(vr0) || !std::isfinite(vi0) || !std::isfinite(vt0) || !std::isfinite(vmeas_safe0) + || !std::isfinite(ipcmd0) || !std::isfinite(iqcmd0) || !std::isfinite(pmeas0) + || !std::isfinite(qgen0) || !std::isfinite(vref0)) + { + Log::error() << "Reecb: initial bus, command, and feedback values must be finite\n"; + return 1; + } + if (vt0 <= ZERO) + { + Log::error() << "Reecb: initial terminal-voltage magnitude must be positive\n"; + return 1; + } + + if (ipcmd0 < ZERO) + { + Log::error() << "Reecb: initial active-current command must be non-negative\n"; + return 1; + } + + const RealT verr0 = Math::deadband2(vref0 - vmeas0, dbd1_, dbd2_); + const RealT iqv0 = Math::clamp(kqv_ * verr0, Iql1_, Iqh1_); + const RealT iqabs0 = std::abs(iqcmd0); + RealT iqneed0 = iqabs0; + if (QFlag_ && iqabs0 > ZERO) + { + iqneed0 += std::numbers::ln2_v / Math::MU + INITIALIZATION_TOLERANCE; + } + + const RealT high0 = pq_on_ * ipcmd0 + pq_off_ * iqabs0; + const RealT low0 = pq_on_ * iqneed0 + pq_off_ * ipcmd0; + // Q priority uses Imax directly for reactive current, so include the + // smooth-clamp recovery margin carried by iqneed0. + const RealT imax = solveInitialLimit( + std::max({Imax_, high0, low0, iqneed0}), high0, low0); + const RealT ilmax0 = circleState(imax, high0); + const RealT ilcap0 = capacity(ilmax0); + if (!std::isfinite(imax) || !std::isfinite(ilmax0) + || !std::isfinite(ilcap0) || ilcap0 < low0) + { + Log::error() << "Reecb: adjusted Imax cannot include the initial current commands\n"; + return 1; + } + + const RealT iqmax0 = pq_on_ * ilcap0 + pq_off_ * imax; + const RealT ipmax0 = pq_on_ * imax + pq_off_ * ilcap0; + + RealT ipraw0 = ZERO; + RealT iqraw0 = ZERO; + if (!iclamp(ipcmd0, ZERO, ipmax0, ipraw0) + || !iclamp(iqcmd0, -iqmax0, iqmax0, iqraw0)) + { + Log::error() << "Reecb: initial current commands cannot be reproduced by their limiters\n"; + return 1; + } + + const RealT iqctl0 = iqraw0 - iqv0; + const RealT pord0 = vmeas_safe0 * ipraw0; + RealT qmin = Qmin_; + RealT qmax = Qmax_; + RealT vmin = Vmin_; + RealT vmax = Vmax_; + if (q_pi_on_ != ZERO) + { + qmin = std::min(Qmin_, qgen0); + qmax = std::max(Qmax_, qgen0); + vmin = std::min(Vmin_, vmeas0); + vmax = std::max(Vmax_, vmeas0); + + const RealT infinity = std::numeric_limits::infinity(); + if (qmin == qgen0 && qmin < qmax) + { + qmin = std::nextafter(qmin, -infinity); + } + if (qmax == qgen0 && qmin < qmax) + { + qmax = std::nextafter(qmax, infinity); + } + if (vmin == vmeas0 && vmin < vmax) + { + vmin = std::nextafter(vmin, -infinity); + } + if (vmax == vmeas0 && vmin < vmax) + { + vmax = std::nextafter(vmax, infinity); + } + } + const RealT pmin = std::min(Pmin_, pord0); + const RealT pmax = std::max(Pmax_, pord0); + const RealT pref0_system = toSystemBase(pord0); + + RealT qtarget0 = ZERO; + if (!QFlag_) + { + qtarget0 = iqctl0 * vmeas_safe0; + } + else if (VFlag_ && !iclamp(qgen0, qmin, qmax, qtarget0)) + { + Log::error() << "Reecb: reactive-power limiter has no finite steady input\n"; + return 1; + } + + RealT qref0 = ZERO; + RealT qext0_port = ZERO; + RealT pfaref0 = ZERO; + + if (v_ref_on_ != ZERO) + { + qext0_port = vmeas0; + } + else if (PfFlag_) + { + if (pmeas0 == ZERO && qtarget0 != ZERO) + { + Log::error() << "Reecb: power-factor mode cannot reproduce the reactive target at zero active power\n"; + return 1; + } + if (pmeas0 != ZERO) + { + pfaref0 = std::atan(qtarget0 / pmeas0); + } + qref0 = pmeas0 * std::tan(pfaref0); + if (std::abs(qref0 - qtarget0) > std::abs(qtarget0) * INITIALIZATION_TOLERANCE) + { + Log::error() << "Reecb: power-factor angle cannot reproduce the reactive target\n"; + return 1; + } + } + else + { + qext0_port = toSystemBase(qtarget0); + qref0 = toComponentBase(qext0_port); + } + + const RealT eq0 = Math::clamp(qref0, qmin, qmax) - qgen0; + RealT xpiq0 = ZERO; + if (q_pi_on_ != ZERO) + { + RealT vpiq_input0 = ZERO; + if (!iclamp(vmeas0, vmin, vmax, vpiq_input0)) + { + Log::error() << "Reecb: voltage limiter has no finite steady input\n"; + return 1; + } + xpiq0 = vpiq_input0 - Kqp_ * eq0; + } + + const RealT vpiq0 = Math::clamp(Kqp_ * eq0 + xpiq0, vmin, vmax); + const RealT epiv0 = q_pi_on_ * vpiq0 + v_ref_on_ * qext0_port - q_on_ * vmeas0; + RealT qv0 = ZERO; + RealT xpiv0 = ZERO; + + if (QFlag_) + { + if (iqmax0 <= INITIALIZATION_TOLERANCE) + { + xpiv0 = -Kvp_ * epiv0; + } + else + { + RealT iqctl_input0 = ZERO; + if (!iclamp(iqctl0, -iqmax0, iqmax0, iqctl_input0)) + { + Log::error() << "Reecb: voltage-controller current cannot be reproduced by its limiter\n"; + return 1; + } + xpiv0 = iqctl_input0 - Kvp_ * epiv0; + } + } + else + { + qv0 = qref0 / vmeas_safe0; + } + + const RealT sdip0 = Math::inside(vt0, Vdip_, Vup_); + const RealT qrate0 = q_pi_on_ * sdip0 * Math::antiwindup(Kqp_ * eq0 + xpiq0, Kqi_ * eq0, vmin, vmax); + const ScalarT vstate0{Kvp_ * epiv0 + xpiv0}; + const ScalarT vderiv0{Kvi_ * epiv0}; + const RealT vrate0 = q_on_ * sdip0 * static_cast(awband(vstate0, vderiv0, ScalarT{iqmax0})); + const RealT iqbase0 = Math::clamp(Kvp_ * epiv0 + xpiv0, -iqmax0, iqmax0); + const RealT iqcmd_check = Math::clamp(q_on_ * iqbase0 + q_off_ * qv0 + iqv0, -iqmax0, iqmax0); + const RealT ipcmd_check = Math::clamp(pord0 / vmeas_safe0, ZERO, ipmax0); + + if (!std::isfinite(imax) || !std::isfinite(ilmax0) || !std::isfinite(ilcap0) + || !std::isfinite(iqmax0) || !std::isfinite(ipmax0) + || !std::isfinite(ipraw0) || !std::isfinite(iqraw0) || !std::isfinite(pord0) + || !std::isfinite(pref0_system) || !std::isfinite(qtarget0) || !std::isfinite(qref0) + || !std::isfinite(qext0_port) || !std::isfinite(pfaref0) || !std::isfinite(eq0) + || !std::isfinite(xpiq0) || !std::isfinite(epiv0) || !std::isfinite(xpiv0) + || !std::isfinite(qv0) || !std::isfinite(qrate0) || !std::isfinite(vrate0) + || !std::isfinite(iqcmd_check) || !std::isfinite(ipcmd_check)) + { + Log::error() << "Reecb: initialization produced a nonfinite value\n"; + return 1; + } + if (std::abs(qrate0) > INITIALIZATION_TOLERANCE || std::abs(vrate0) > INITIALIZATION_TOLERANCE) + { + Log::error() << "Reecb: controller state rate is nonzero at initialization\n"; + return 1; + } + if (std::abs(iqcmd_check - iqcmd0) > INITIALIZATION_TOLERANCE + || std::abs(ipcmd_check - ipcmd0) > INITIALIZATION_TOLERANCE) + { + Log::error() << "Reecb: current-command limiter reconstruction is inexact\n"; + return 1; + } + + const bool q_adjusted = qmin != Qmin_ || qmax != Qmax_; + const bool v_adjusted = vmin != Vmin_ || vmax != Vmax_; + const bool p_adjusted = pmin != Pmin_ || pmax != Pmax_; + const bool imax_adjusted = imax != Imax_; + + Qmin_ = qmin; + Qmax_ = qmax; + Vmin_ = vmin; + Vmax_ = vmax; + Pmin_ = pmin; + Pmax_ = pmax; + Imax_ = imax; + + y[VMEAS] = vmeas0; + y[PMEAS] = pmeas0; + y[XPIQ] = xpiq0; + y[XPIV] = xpiv0; + y[QV] = qv0; + y[PORD] = pord0; + y[VT] = vt0; + y[ILMAX] = ilmax0; + + if (!Vref0_given_) + { + Vref0_ = vref0; + } + + pe_set_ = static_cast(pe0_system); + qgen_set_ = static_cast(qgen0_system); + qext_set_ = static_cast(qext0_port); + pfaref_set_ = static_cast(pfaref0); + pref_set_ = static_cast(pref0_system); + + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(qext_set_); + } + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(pfaref_set_); + } + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(pref_set_); + } + + if (q_adjusted) + { + Log::warning() << "Reecb: Qmin/Qmax adjusted to include the initial reactive power\n"; + } + if (v_adjusted) + { + Log::warning() << "Reecb: Vmin/Vmax adjusted to include the initial terminal voltage\n"; + } + if (p_adjusted) + { + Log::warning() << "Reecb: Pmin/Pmax adjusted to include the initial active-power order\n"; + } + if (imax_adjusted) + { + Log::warning() << "Reecb: Imax adjusted to include the initial current commands\n"; + } + + y_.setDataUpdated(); + yp_.setToConst(static_cast(ZERO)); + return 0; + } + + /** + * @brief Identify the differential variables + * + * The two measurement filters, two PI states, reactive-current lag, and + * active-power order carry derivatives; all other rows are algebraic. + */ + template + int Reecb::tagDifferentiable() + { + const auto VMEAS = static_cast(ReecbInternalVariables::VMEAS); + const auto PMEAS = static_cast(ReecbInternalVariables::PMEAS); + const auto XPIQ = static_cast(ReecbInternalVariables::XPIQ); + const auto XPIV = static_cast(ReecbInternalVariables::XPIV); + const auto QV = static_cast(ReecbInternalVariables::QV); + const auto PORD = static_cast(ReecbInternalVariables::PORD); + + std::fill(tag_.begin(), tag_.end(), false); + tag_[VMEAS] = true; + tag_[PMEAS] = true; + tag_[XPIQ] = true; + tag_[XPIV] = true; + tag_[QV] = true; + tag_[PORD] = true; + return 0; + } + + /** + * @brief Set one absolute tolerance for every REECB variable + * + * @param[in] rel_tol Solver relative tolerance used as the absolute floor. + */ + template + int Reecb::setAbsoluteTolerance(RealT rel_tol) + { + abs_tol_.setToConst(static_cast(rel_tol)); + return 0; + } + + /** + * @brief Evaluate the model residuals + * + * Starts from latched values, refreshes attached signals and their indices, + * refreshes terminal-bus voltage, and evaluates the internal residual. + * REECB contributes no bus residual. + */ + template + int Reecb::evaluateResidual() + { + const auto PE = static_cast(ReecbExternalVariables::PE); + const auto QGEN = static_cast(ReecbExternalVariables::QGEN); + const auto QEXT = static_cast(ReecbExternalVariables::QEXT); + const auto PFAREF = static_cast(ReecbExternalVariables::PFAREF); + const auto PREF = static_cast(ReecbExternalVariables::PREF); + + ws_[PE] = pe_set_; + ws_[QGEN] = qgen_set_; + ws_[QEXT] = qext_set_; + ws_[PFAREF] = pfaref_set_; + ws_[PREF] = pref_set_; + std::fill(ws_indices_.begin(), ws_indices_.end(), INVALID_INDEX); + + if (signals_.template isAttached()) + { + ws_[PE] = signals_.template readExternalVariable(); + ws_indices_[PE] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[QGEN] = signals_.template readExternalVariable(); + ws_indices_[QGEN] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[QEXT] = signals_.template readExternalVariable(); + ws_indices_[QEXT] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[PFAREF] = + signals_.template readExternalVariable(); + ws_indices_[PFAREF] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[PREF] = signals_.template readExternalVariable(); + ws_indices_[PREF] = + signals_.template readExternalVariableIndex(); + } + + wb_[0] = Vr(); + wb_[1] = Vi(); + + evaluateInternalResidual(y_.getData(), yp_.getData(), wb_.data(), ws_.data(), f_.getData()); + f_.setDataUpdated(); + return 0; + } + + /** + * @brief Access the REECB signal interface + * + * @return Interface used to assign current-command outputs and attach + * optional feedback and reference signals. + */ + template + auto Reecb::getSignals() + -> ComponentSignals& + { + return signals_; + } + + /** + * @brief Access the optional variable monitor + * + * @return Monitor, or nullptr when constructed without model data. + */ + template + const Model::VariableMonitorBase* Reecb::getMonitor() const + { + return monitor_.get(); + } + + /** + * @brief Evaluate the REECB internal residual + * + * The branch-free equation body preserves a fixed dependency structure; + * parameter-selected paths enter through selector masks resolved by + * setDerivedParameters(). The `ILMAX` row uses a smooth signed-square + * continuation, while a smooth magnitude supplies the limiter bounds. + * + * @param[in] y Internal variables. + * @param[in] yp Internal variable derivatives. + * @param[in] wb Terminal-bus voltage components. + * @param[in] ws External signal values in their documented port units and bases. + * @param[out] f Internal residuals. + */ + template + [[gnu::always_inline]] inline int + Reecb::evaluateInternalResidual( + const ScalarT* y, + const ScalarT* yp, + const ScalarT* wb, + const ScalarT* ws, + ScalarT* f) + { + const auto VMEAS = static_cast(ReecbInternalVariables::VMEAS); + const auto PMEAS = static_cast(ReecbInternalVariables::PMEAS); + const auto XPIQ = static_cast(ReecbInternalVariables::XPIQ); + const auto XPIV = static_cast(ReecbInternalVariables::XPIV); + const auto QV = static_cast(ReecbInternalVariables::QV); + const auto PORD = static_cast(ReecbInternalVariables::PORD); + const auto VT = static_cast(ReecbInternalVariables::VT); + const auto ILMAX = static_cast(ReecbInternalVariables::ILMAX); + const auto IQCMD = static_cast(ReecbInternalVariables::IQCMD); + const auto IPCMD = static_cast(ReecbInternalVariables::IPCMD); + + const auto PE = static_cast(ReecbExternalVariables::PE); + const auto QGEN = static_cast(ReecbExternalVariables::QGEN); + const auto QEXT = static_cast(ReecbExternalVariables::QEXT); + const auto PFAREF = static_cast(ReecbExternalVariables::PFAREF); + const auto PREF = static_cast(ReecbExternalVariables::PREF); + + const ScalarT vmeas = y[VMEAS]; + const ScalarT pmeas = y[PMEAS]; + const ScalarT xpiq = y[XPIQ]; + const ScalarT xpiv = y[XPIV]; + const ScalarT qv = y[QV]; + const ScalarT pord = y[PORD]; + const ScalarT vt = y[VT]; + const ScalarT ilmax = y[ILMAX]; + const ScalarT iqcmd_system = y[IQCMD]; + const ScalarT ipcmd_system = y[IPCMD]; + + const ScalarT vmeas_dot = yp[VMEAS]; + const ScalarT pmeas_dot = yp[PMEAS]; + const ScalarT xpiq_dot = yp[XPIQ]; + const ScalarT xpiv_dot = yp[XPIV]; + const ScalarT qv_dot = yp[QV]; + const ScalarT pord_dot = yp[PORD]; + + const ScalarT vr = wb[0]; + const ScalarT vi = wb[1]; + + const ScalarT pe = toComponentBase(ws[PE]); + const ScalarT qgen = toComponentBase(ws[QGEN]); + const ScalarT extref = ws[QEXT]; + const ScalarT pfaref = ws[PFAREF]; + const ScalarT pref = toComponentBase(ws[PREF]); + const ScalarT iqcmd = toComponentBase(iqcmd_system); + const ScalarT ipcmd = toComponentBase(ipcmd_system); + + const ScalarT vmeas_safe = Math::max(vmeas, VMEAS_MINIMUM); + const ScalarT sdip = Math::inside(vt, Vdip_, Vup_); + const ScalarT verr = Math::deadband2(Vref0_ - vmeas, dbd1_, dbd2_); + const ScalarT iqv = Math::clamp(kqv_ * verr, Iql1_, Iqh1_); + // The Volt/VAr channel is a system-base reactive power unless + // direct-voltage mode selects it as a terminal-voltage reference, + // which takes no power-base conversion. + const ScalarT qref = q_ref_on_ * (pf_on_ * pmeas * std::tan(pfaref) + pf_off_ * toComponentBase(extref)); + const ScalarT eq = Math::clamp(qref, Qmin_, Qmax_) - qgen; + const ScalarT vpiq = Math::clamp(Kqp_ * eq + xpiq, Vmin_, Vmax_); + const ScalarT epiv = q_pi_on_ * vpiq + v_ref_on_ * extref - q_on_ * vmeas; + const ScalarT fpord = (pref - pord) / Tpord_; + const ScalarT rpord = aslew(fpord, dPmin_, dPmax_); + const ScalarT ilnorm = std::sqrt(ilmax * ilmax + INITIALIZATION_TOLERANCE); + const ScalarT ilcap = (ilmax / ilnorm) * ilmax; + const ScalarT high = pq_on_ * ipcmd + pq_off_ * iqcmd; + const ScalarT iqmax = pq_on_ * ilcap + pq_off_ * Imax_; + const ScalarT ipmax = pq_on_ * Imax_ + pq_off_ * ilcap; + const ScalarT iqbase = Math::clamp(Kvp_ * epiv + xpiv, -iqmax, iqmax); + const ScalarT iqraw = q_on_ * iqbase + q_off_ * qv + iqv; + + f[VMEAS] = -vmeas_dot + (vt - vmeas) / Trv_; + f[PMEAS] = -pmeas_dot + (pe - pmeas) / Tp_; + f[XPIQ] = -xpiq_dot + q_pi_on_ * sdip * Math::antiwindup(Kqp_ * eq + xpiq, Kqi_ * eq, Vmin_, Vmax_); + f[XPIV] = -xpiv_dot + q_on_ * sdip * awband(Kvp_ * epiv + xpiv, Kvi_ * epiv, iqmax); + f[QV] = -qv_dot + q_off_ * sdip * (qref / vmeas_safe - qv) / Tiq_; + f[PORD] = -pord_dot + sdip * Math::antiwindup(pord, rpord, Pmin_, Pmax_); + f[VT] = -vt * vt + vr * vr + vi * vi; + f[ILMAX] = -ilmax * ilnorm + (Imax_ - high) * (Imax_ + high); + f[IQCMD] = -iqcmd + Math::clamp(iqraw, -iqmax, iqmax); + f[IPCMD] = -ipcmd + Math::clamp(pord / vmeas_safe, ZERO, ipmax); + + return 0; + } + + // + // Private methods + // + + /** + * @brief Smooth asymmetric slew-rate limiter + * + * @param[in] rate Unconstrained rate. + * @param[in] lower Negative rate limit. + * @param[in] upper Positive rate limit. + * @return Limited rate. + */ + template + [[gnu::always_inline]] inline scalar_type + Reecb::aslew(ScalarT rate, RealT lower, RealT upper) + { + assert(lower < ZERO && ZERO < upper); + return rate + / (ONE + Math::ramp(rate / upper - ONE) + Math::ramp(rate / lower - ONE)); + } + + /** + * @brief Smooth anti-windup derivative within a moving symmetric band + * + * Math::antiwindup over [-band, band] with a band edge that is an + * algebraic quantity, so differentiation carries the band's own + * contributions through the gate. + * + * @param[in] state Limited PI state. + * @param[in] rate Pre-limit derivative of state. + * @param[in] band Nonnegative symmetric band edge. + * @return Anti-windup-limited derivative. + * + * @todo Fold moving-limit support into Math::antiwindup in CommonMath. + */ + template + [[gnu::always_inline]] inline scalar_type + Reecb::awband(ScalarT state, ScalarT rate, ScalarT band) + { + const ScalarT above_min = Math::above(state, -band); + const ScalarT below_max = Math::below(state, band); + return (above_min * below_max + (ONE - below_max) * Math::sigmoid(-rate) + + (ONE - above_min) * Math::sigmoid(rate)) + * rate; + } + + /** + * @brief Compute the initial current-circle continuation state + * + * Solves the implemented `ILMAX` row for its nonnegative continuation + * state at a total-current limit and priority-axis current. + * + * @param[in] imax Total-current limit on the component base. + * @param[in] high Priority-axis current command on the component base. + * @return Initial `ILMAX` state on the component base, or a quiet NaN + * when the circle geometry is invalid or unrepresentable. + * @pre `imax` and `high` are finite and @f$0 \le high \le imax@f$. + */ + template + typename Reecb::RealT + Reecb::circleState(RealT imax, RealT high) + { + const RealT rhs = (imax - high) * (imax + high); + if (!std::isfinite(rhs) || rhs < ZERO) + { + return std::numeric_limits::quiet_NaN(); + } + if (rhs == ZERO) + { + return ZERO; + } + + const RealT ratio = INITIALIZATION_TOLERANCE / rhs; + return std::sqrt(rhs) + * std::sqrt(TWO + / (std::hypot(ratio, TWO) + ratio)); + } + + /** + * @brief Compute off-axis capacity from a continuation state + * + * This expression must match the `ilcap` calculation during residual + * evaluation so initialization lands on the implemented model. + * + * @param[in] ilmax Current-circle continuation state on the component base. + * @return Available off-axis current on the component base. + */ + template + typename Reecb::RealT + Reecb::capacity(RealT ilmax) + { + const RealT ilnorm = std::sqrt(ilmax * ilmax + INITIALIZATION_TOLERANCE); + return (ilmax / ilnorm) * ilmax; + } + + /** + * @brief Bisect an initial-limit interval to machine rounding + * + * The upper endpoint is returned because initialization needs the first + * not-below point; the caller separately validates that point as finite + * and feasible. A finite floating-point interval contains finitely many + * representable values, so the loop terminates when no interior midpoint + * remains. + * + * @tparam FuncT Monotone predicate type. + * @param[in] a Lower endpoint, where `below` is true. + * @param[in] b Upper endpoint, where `below` is false. + * @param[in] below Predicate returning true only when finite reconstructed + * capacity lies below the requirement. + * @return The first representable upper-side endpoint. + * @pre `a` and `b` are finite, `a < b`, and `below` is monotone. + */ + template + template + typename Reecb::RealT + Reecb::bisect(RealT a, RealT b, FuncT below) + { + while (true) + { + const RealT mid = std::midpoint(a, b); + if (mid <= a || b <= mid) + { + break; + } + + if (below(mid)) + { + a = mid; + } + else + { + b = mid; + } + } + + return b; + } + + /** + * @brief Solve the smallest initial total-current limit + * + * @param[in] lower Lower bound for the component-base total-current limit. + * @param[in] high Priority-axis current command on the component base. + * @param[in] low Required off-axis capacity on the component base. + * @return Smallest representable feasible component-base limit at or above + * `lower`, or a quiet NaN when no finite limit is found. + * @pre The arguments are finite and nonnegative, with `lower >= high` and + * `lower >= low`. + * @warning This function contains conditional branching and may be used + * during initialization, but not during residual evaluation. + */ + template + typename Reecb::RealT + Reecb::solveInitialLimit( + RealT lower, RealT high, RealT low) + { + const RealT nan = std::numeric_limits::quiet_NaN(); + + const RealT ilmax = circleState(lower, high); + const RealT ilcap = capacity(ilmax); + if (!std::isfinite(ilcap)) + { + return nan; + } + if (!(ilcap < low)) + { + return lower; + } + + const auto below = [high, low](RealT limit) + { + const RealT state = circleState(limit, high); + const RealT cap = capacity(state); + return std::isfinite(cap) && cap < low; + }; + + // Invert ilcap = ilmax^2 / sqrt(ilmax^2 + tolerance), then recover + // imax from imax^2 = high^2 + ilmax * sqrt(ilmax^2 + tolerance). + const RealT delta = std::sqrt(INITIALIZATION_TOLERANCE); + const RealT ilreq = std::sqrt(low) + * std::sqrt(HALF + * (low + std::hypot(low, TWO * delta))); + const RealT seed = std::hypot( + high, std::sqrt(ilreq) * std::sqrt(std::hypot(ilreq, delta))); + + const RealT maximum = std::numeric_limits::max(); + RealT a = lower; + RealT b = seed; + if (!std::isfinite(b)) + { + b = maximum; + } + else + { + b = std::max(b, std::nextafter(a, maximum)); + } + if (!(a < b)) + { + return nan; + } + + while (below(b)) + { + a = b; + if (b >= maximum) + { + return nan; + } + if (b > maximum / TWO) + { + b = maximum; + } + else + { + b *= TWO; + } + } + + const RealT result = bisect(a, b, below); + const RealT final_state = circleState(result, high); + const RealT final_cap = capacity(final_state); + if (!std::isfinite(result) || !std::isfinite(final_state) + || !std::isfinite(final_cap) || final_cap < low) + { + return nan; + } + return result; + } + + /** + * @brief Load one real-valued parameter + * + * Real and integer serialized values are accepted. Any other stored type + * records a loading error while preserving the existing value. + * + * @param[in] data Model parameter data. + * @param[in] parameter Parameter key to load. + * @param[in,out] target Stored parameter value. + * @param[in] name Serialized parameter name for diagnostics. + */ + template + void Reecb::loadRealParameter( + const ModelDataT& data, + ReecbParameters parameter, + RealT& target, + const char* name) + { + if (!data.parameters.contains(parameter)) + { + return; + } + + const auto& value = data.parameters.at(parameter); + 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() << "Reecb: parameter '" << name << "' must be numeric\n"; + ++parameter_error_count_; + } + } + + /** + * @brief Load one optional Boolean parameter + * + * Any non-Boolean stored type records a loading error while preserving + * the existing default. + * + * @param[in] data Model parameter data. + * @param[in] parameter Parameter key to load. + * @param[in,out] target Stored Boolean value. + * @param[in] name Serialized parameter name for diagnostics. + */ + template + void Reecb::loadBooleanParameter( + const ModelDataT& data, + ReecbParameters parameter, + bool& target, + const char* name) + { + if (!data.parameters.contains(parameter)) + { + return; + } + + const auto& value = data.parameters.at(parameter); + if (const auto* bool_value = std::get_if(&value)) + { + target = *bool_value; + } + else + { + Log::error() << "Reecb: parameter '" << name << "' must be boolean\n"; + ++parameter_error_count_; + } + } + + /** + * @brief Validate and floor one explicit controller lag + * + * Nonfinite and negative values record errors before replacement so + * verify() retains the evidence. A valid value below 1 ms is raised and + * reported through the return value, preserving an explicit Hessenberg + * residual and avoiding division by zero. + * + * @param[in,out] value Time constant to validate and floor. + * @param[in] name Parameter name for diagnostics. + * @return true only when a valid nonnegative value was raised. + */ + template + bool Reecb::floorTimeConstant( + RealT& value, const char* name) + { + if (!std::isfinite(value)) + { + Log::error() << "Reecb: " << name << " must be finite\n"; + ++parameter_error_count_; + value = TIME_CONSTANT_MINIMUM; + return false; + } + if (value < ZERO) + { + Log::error() << "Reecb: " << name << " must be non-negative\n"; + ++parameter_error_count_; + value = TIME_CONSTANT_MINIMUM; + return false; + } + + const bool raised = value < TIME_CONSTANT_MINIMUM; + value = std::max(value, TIME_CONSTANT_MINIMUM); + return raised; + } + + /** + * @brief Read parameters from model data + * + * Omitted parameters retain their documented defaults. Loading errors + * are counted for verify() rather than thrown. + * + * @param[in] data Parameters and monitored-variable selections. + */ + template + void Reecb::initializeParameters(const ModelDataT& data) + { + using Params = typename ModelDataT::Parameters; + + parameter_error_count_ = 0; + mva_given_ = data.parameters.contains(Params::mva); + Vref0_given_ = false; + + loadRealParameter(data, Params::mva, mva_base_, "mva"); + loadBooleanParameter(data, Params::PfFlag, PfFlag_, "PfFlag"); + loadBooleanParameter(data, Params::VFlag, VFlag_, "VFlag"); + loadBooleanParameter(data, Params::QFlag, QFlag_, "QFlag"); + loadBooleanParameter(data, Params::Pqflag, Pqflag_, "Pqflag"); + loadRealParameter(data, Params::Trv, Trv_, "Trv"); + loadRealParameter(data, Params::Tp, Tp_, "Tp"); + if (data.parameters.contains(Params::Vref0)) + { + loadRealParameter(data, Params::Vref0, Vref0_, "Vref0"); + Vref0_given_ = true; + } + loadRealParameter(data, Params::Vdip, Vdip_, "Vdip"); + loadRealParameter(data, Params::Vup, Vup_, "Vup"); + loadRealParameter(data, Params::dbd1, dbd1_, "dbd1"); + loadRealParameter(data, Params::dbd2, dbd2_, "dbd2"); + loadRealParameter(data, Params::kqv, kqv_, "kqv"); + loadRealParameter(data, Params::Iql1, Iql1_, "Iql1"); + loadRealParameter(data, Params::Iqh1, Iqh1_, "Iqh1"); + loadRealParameter(data, Params::Qmax, Qmax_, "Qmax"); + loadRealParameter(data, Params::Qmin, Qmin_, "Qmin"); + loadRealParameter(data, Params::Kqp, Kqp_, "Kqp"); + loadRealParameter(data, Params::Kqi, Kqi_, "Kqi"); + loadRealParameter(data, Params::Vmax, Vmax_, "Vmax"); + loadRealParameter(data, Params::Vmin, Vmin_, "Vmin"); + loadRealParameter(data, Params::Kvp, Kvp_, "Kvp"); + loadRealParameter(data, Params::Kvi, Kvi_, "Kvi"); + loadRealParameter(data, Params::Tiq, Tiq_, "Tiq"); + loadRealParameter(data, Params::Tpord, Tpord_, "Tpord"); + loadRealParameter(data, Params::dPmax, dPmax_, "dPmax"); + loadRealParameter(data, Params::dPmin, dPmin_, "dPmin"); + loadRealParameter(data, Params::Pmax, Pmax_, "Pmax"); + loadRealParameter(data, Params::Pmin, Pmin_, "Pmin"); + loadRealParameter(data, Params::Imax, Imax_, "Imax"); + + setDerivedParameters(); + } + + /** + * @brief Bind monitor selections to REECB internal states + */ + template + void Reecb::initializeMonitor() + { + using Variable = typename ModelDataT::MonitorableVariables; + + monitor_->set(Variable::iqcmd, [this] + { return y_.getData()[static_cast(ReecbInternalVariables::IQCMD)]; }); + monitor_->set(Variable::ipcmd, [this] + { return y_.getData()[static_cast(ReecbInternalVariables::IPCMD)]; }); + monitor_->set(Variable::vmeas, [this] + { return y_.getData()[static_cast(ReecbInternalVariables::VMEAS)]; }); + monitor_->set(Variable::pmeas, [this] + { return y_.getData()[static_cast(ReecbInternalVariables::PMEAS)]; }); + } + + /** + * @brief Resolve parameter-derived constants and selector masks + * + * Raises explicit controller lags in place, converts any supplied component + * rating, and resolves selector masks. Invalid lag inputs are + * recorded before replacement so verify() retains each error. + */ + template + void Reecb::setDerivedParameters() + { + bool floor_warning = false; + + floor_warning |= floorTimeConstant(Trv_, "Trv"); + floor_warning |= floorTimeConstant(Tp_, "Tp"); + floor_warning |= floorTimeConstant(Tiq_, "Tiq"); + floor_warning |= floorTimeConstant(Tpord_, "Tpord"); + + if (floor_warning) + { + Log::warning() << "Reecb: any of Trv, Tp, Tiq, or Tpord below " + << TIME_CONSTANT_MINIMUM + << " s is raised to that floor to keep the controller lags well posed\n"; + } + + va_component_base_ = mva_base_ * static_cast(1.0e6); + + if (PfFlag_ && QFlag_) + { + Log::warning() << "Reecb: PfFlag and QFlag are both enabled; " + << "this is an atypical control configuration\n"; + } + + pf_on_ = ZERO; + if (PfFlag_) + { + pf_on_ = ONE; + } + pf_off_ = ONE - pf_on_; + + q_on_ = ZERO; + if (QFlag_) + { + q_on_ = ONE; + } + q_off_ = ONE - q_on_; + + q_pi_on_ = ZERO; + if (QFlag_ && VFlag_) + { + q_pi_on_ = ONE; + } + + v_ref_on_ = ZERO; + if (QFlag_ && !VFlag_) + { + v_ref_on_ = ONE; + } + q_ref_on_ = ONE - v_ref_on_; + + pq_on_ = ZERO; + if (Pqflag_) + { + pq_on_ = ONE; + } + pq_off_ = ONE - pq_on_; + } + + /** + * @brief Evaluate log(1 - exp(-x)) accurately for positive x + * + * The small-x hyperbolic form avoids cancellation; the large-x form uses + * log1p. The algebraically equivalent branches agree in value and first + * derivative at x = log(2). + * + * @param[in] x Strictly positive argument. + * @return Numerically stable value of log(1 - exp(-x)). + * @pre `x > 0`. + * @warning This function contains conditional branching and as such can + * be used in initialization methods but not in residual evaluation. + */ + template + typename Reecb::RealT + Reecb::logOneMinusExp(RealT x) + { + static constexpr auto log_two = std::numbers::ln2_v; + + if (x < log_two) + { + return log_two - HALF * x + std::log(std::sinh(HALF * x)); + } + return std::log1p(-std::exp(-x)); + } + + /** + * @brief Recover the input that produces a requested smooth-clamp output + * + * Exact bounds use a finite offset derived from the initialization + * tolerance; collapsed bounds are reproduced directly. + * + * @param[in] output Requested output. + * @param[in] lower Lower smooth-clamp limit. + * @param[in] upper Upper smooth-clamp limit. + * @param[out] input Recovered clamp input on success. + * @return true when a finite admissible input was recovered. + * @warning This function contains conditional branching and as such can + * be used in initialization methods but not in residual evaluation. + */ + template + bool Reecb::iclamp(RealT output, RealT lower, RealT upper, RealT& input) const + { + if (!std::isfinite(output) || !std::isfinite(lower) || !std::isfinite(upper) || lower > upper + || output < lower - INITIALIZATION_TOLERANCE || output > upper + INITIALIZATION_TOLERANCE) + { + return false; + } + + output = std::clamp(output, lower, upper); + if (upper == lower) + { + input = lower; + return true; + } + + const RealT mu = Math::MU; + const RealT offset = -std::log(std::expm1(mu * HALF * INITIALIZATION_TOLERANCE)) / mu; + if (output == lower) + { + input = lower - offset; + return true; + } + if (output == upper) + { + input = upper + offset; + return true; + } + + const RealT a = mu * (output - lower); + const RealT b = mu * (upper - output); + input = lower + (a + logOneMinusExp(a) - logOneMinusExp(b)) / mu; + return std::isfinite(input); + } + + /** + * @brief Resolve the REECB component power base + * + * @return The supplied component base, or the system base when `mva` is omitted. + */ + template + typename Reecb::RealT + Reecb::componentPowerBase() const + { + if (mva_given_) + { + return va_component_base_; + } + return va_system_base_; + } + + /** + * @brief Convert a system-base power or current to REECB component base + * + * @param[in] value Quantity on the system base. + * @return The same quantity on the REECB component base. + */ + template + template + [[gnu::always_inline]] inline ValueT + Reecb::toComponentBase(ValueT value) const + { + return value * (va_system_base_ / componentPowerBase()); + } + + /** + * @brief Convert a component-base power or current to system base + * + * @param[in] value Quantity on the REECB component base. + * @return The same quantity on the system base. + */ + template + template + ValueT Reecb::toSystemBase(ValueT value) const + { + return value * (componentPowerBase() / va_system_base_); + } + + /** + * @brief Access the terminal-bus real voltage component + * + * @return Mutable reference to the bus real voltage state. + */ + template + scalar_type& Reecb::Vr() + { + return bus_->Vr(); + } + + /** + * @brief Access the terminal-bus imaginary voltage component + * + * @return Mutable reference to the bus imaginary voltage state. + */ + template + scalar_type& Reecb::Vi() + { + return bus_->Vi(); + } + } // namespace Controller + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/README.md b/GridKit/Model/PhasorDynamics/Converter/README.md index ad38ba19a..aa937a8b2 100644 --- a/GridKit/Model/PhasorDynamics/Converter/README.md +++ b/GridKit/Model/PhasorDynamics/Converter/README.md @@ -2,8 +2,9 @@ ## Introduction -Converter models represent inverter-coupled resources in the phasor dynamics model. They provide the network interface between renewable-energy control -models and the bus equations, typically through commanded active and reactive current components. +Converter models represent inverter-coupled resources in the phasor dynamics +model and provide the network interface between renewable-energy controller +models and the bus equations. ## Types diff --git a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md index 8ae69fcb1..ce4a114e5 100644 --- a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md +++ b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md @@ -152,6 +152,7 @@ are specified: [Gensal](SynchronousMachine/GENSAL/README.md) | 5th order salient-pole machine model [GenClassical](SynchronousMachine/GenClassical/README.md) | the classical machine model [Regca](Converter/REGCA/README.md) | WECC REGCA renewable generator/converter model + [Reecb](Controller/REECB/README.md) | WECC REECB renewable electrical-control model [Repca](Controller/REPCA/README.md) | the REPCA renewable plant-control model [Tgov1](Governor/Tgov1/README.md) | the TGOV1 governor model [Hygov](Governor/HYGOV/README.md) | the HYGOV hydro turbine-governor model diff --git a/GridKit/Model/PhasorDynamics/SignalSource/ConstantSignalSourceImpl.hpp b/GridKit/Model/PhasorDynamics/SignalSource/ConstantSignalSourceImpl.hpp index a541dbe00..5310d4a88 100644 --- a/GridKit/Model/PhasorDynamics/SignalSource/ConstantSignalSourceImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SignalSource/ConstantSignalSourceImpl.hpp @@ -112,10 +112,13 @@ namespace GridKit return 0; } + /** + * @brief Construct the empty Jacobian for this stateless source. + */ template int ConstantSignalSource::evaluateJacobian() { - return 0; + return this->constructCoo(); } } // namespace PhasorDynamics diff --git a/GridKit/Model/PhasorDynamics/SystemModelData.hpp b/GridKit/Model/PhasorDynamics/SystemModelData.hpp index 224c95e3c..c83d60189 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelData.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelData.hpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -46,6 +47,7 @@ namespace GridKit using BusToSignalAdapterDataT = BusToSignalAdapterData; using BusFaultDataT = BusFaultData; using RegcaDataT = Converter::RegcaData; + using ReecbDataT = Controller::ReecbData; using RepcaDataT = Controller::RepcaData; using Tgov1DataT = Governor::Tgov1Data; using Esdc1aDataT = Exciter::Esdc1aData; @@ -104,6 +106,7 @@ namespace GridKit 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 reecb; ///< REECB electrical controllers within the model std::vector repca; ///< REPCA plant controllers within the model std::vector genrou; ///< GENROU instances within the model std::vector gensal; ///< GENSAL instances within the model diff --git a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp index c1701687a..6847c769f 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp @@ -141,6 +141,12 @@ namespace GridKit raw_component.get_to(regca); sm.regca.push_back(regca); } + else if (kind == "Reecb") + { + typename SystemModelData::ReecbDataT reecb; + raw_component.get_to(reecb); + sm.reecb.push_back(reecb); + } else if (kind == "Repca") { typename SystemModelData::RepcaDataT repca; diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index 1edc3f62c..1076b3e68 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -294,6 +294,64 @@ namespace GridKit addComponent(gen); } + // Add REECB after its current-command and feedback producers because + // components initialize in insertion order. + for (const auto& reecbdata : data.reecb) + { + BusT* bus = nullptr; + if (reecbdata.buses.contains(ReecbBuses::bus)) + { + bus = getBus(reecbdata.buses.at(ReecbBuses::bus)); + } + + auto* reecb = new Reecb(bus, reecbdata); + + if (reecbdata.signal_inputs.contains(ReecbSignalInputs::pe)) + { + const IdxT pe = reecbdata.signal_inputs.at(ReecbSignalInputs::pe); + constexpr auto PE = ReecbExternalVariables::PE; + reecb->getSignals().template attachSignalNode(getSignal(pe)); + } + if (reecbdata.signal_inputs.contains(ReecbSignalInputs::qgen)) + { + const IdxT qgen = reecbdata.signal_inputs.at(ReecbSignalInputs::qgen); + constexpr auto QGEN = ReecbExternalVariables::QGEN; + reecb->getSignals().template attachSignalNode(getSignal(qgen)); + } + if (reecbdata.signal_inputs.contains(ReecbSignalInputs::qext)) + { + const IdxT qext = reecbdata.signal_inputs.at(ReecbSignalInputs::qext); + constexpr auto QEXT = ReecbExternalVariables::QEXT; + reecb->getSignals().template attachSignalNode(getSignal(qext)); + } + if (reecbdata.signal_inputs.contains(ReecbSignalInputs::pfaref)) + { + const IdxT pfaref = reecbdata.signal_inputs.at(ReecbSignalInputs::pfaref); + constexpr auto PFAREF = ReecbExternalVariables::PFAREF; + reecb->getSignals().template attachSignalNode(getSignal(pfaref)); + } + if (reecbdata.signal_inputs.contains(ReecbSignalInputs::pref)) + { + const IdxT pref = reecbdata.signal_inputs.at(ReecbSignalInputs::pref); + constexpr auto PREF = ReecbExternalVariables::PREF; + reecb->getSignals().template attachSignalNode(getSignal(pref)); + } + if (reecbdata.signal_outputs.contains(ReecbSignalOutputs::iqcmd)) + { + const IdxT iqcmd = reecbdata.signal_outputs.at(ReecbSignalOutputs::iqcmd); + constexpr auto IQCMD = ReecbInternalVariables::IQCMD; + reecb->getSignals().template assignSignalNode(getSignal(iqcmd)); + } + if (reecbdata.signal_outputs.contains(ReecbSignalOutputs::ipcmd)) + { + const IdxT ipcmd = reecbdata.signal_outputs.at(ReecbSignalOutputs::ipcmd); + constexpr auto IPCMD = ReecbInternalVariables::IPCMD; + reecb->getSignals().template assignSignalNode(getSignal(ipcmd)); + } + + addComponent(reecb); + } + // Add Tgov1 governors for (const auto& govdata : data.gov) { @@ -657,8 +715,8 @@ namespace GridKit * * @note System model composition is flat; nested systems are not supported. * - * @throws std::runtime_error if storage allocation, child binding, or - * model verification fails. + * @throws std::runtime_error if storage allocation, child binding, model + * verification, or initialization for sparse Jacobian discovery fails. */ template int SystemModel::allocate() @@ -773,21 +831,27 @@ namespace GridKit throw std::runtime_error("SystemModel allocation failed"); } - // Start variable monitors - initializeMonitor(); - startMonitor(); - - // Perform an initial Jacobian evaluation for sparse Jacobians, such that - // the dynamic solver can querry the NNZ value when it is configured. - // @todo Replace with a sparsity analysis that sets the NNZ and allocates the Jacobian - // without needing the Jacobian values. + // Sparse-pattern discovery requires an initialized operating point. A failed + // initialization aborts allocation before residual/Jacobian evaluation or + // monitor startup. + // @todo Replace with a sparsity analysis that sets the NNZ and allocates + // the Jacobian without needing the Jacobian values. if (hasJacobian()) { - initialize(); + const int status = initialize(); + if (status != 0) + { + Log::error() << "System model initialization failed with status " + << status << '\n'; + throw std::runtime_error("SystemModel allocation failed"); + } evaluateResidual(); evaluateJacobian(); } + initializeMonitor(); + startMonitor(); + allocated_ = true; return 0; } diff --git a/docs/Figures/PhasorDynamics/REECB/diagram.png b/docs/Figures/PhasorDynamics/REECB/diagram.png new file mode 100644 index 000000000..c3aaa1478 Binary files /dev/null and b/docs/Figures/PhasorDynamics/REECB/diagram.png differ diff --git a/docs/GridKit/Model/PhasorDynamics/Controller/README.md b/docs/GridKit/Model/PhasorDynamics/Controller/README.md index c7b0dd611..3962c8eb4 100644 --- a/docs/GridKit/Model/PhasorDynamics/Controller/README.md +++ b/docs/GridKit/Model/PhasorDynamics/Controller/README.md @@ -5,6 +5,7 @@ :titlesonly: :hidden: +REECB REPCA ``` diff --git a/docs/GridKit/Model/PhasorDynamics/Controller/REECB/README.md b/docs/GridKit/Model/PhasorDynamics/Controller/REECB/README.md new file mode 100644 index 000000000..6926035a0 --- /dev/null +++ b/docs/GridKit/Model/PhasorDynamics/Controller/REECB/README.md @@ -0,0 +1,6 @@ +# REECB + +```{include} ../../../../../../GridKit/Model/PhasorDynamics/Controller/REECB/README.md +:start-line: 1 +:relative-images: +``` diff --git a/tests/IntegrationTests/PhasorDynamics/PDIntegrationTests.hpp b/tests/IntegrationTests/PhasorDynamics/PDIntegrationTests.hpp index 5f58c410d..cd7011983 100644 --- a/tests/IntegrationTests/PhasorDynamics/PDIntegrationTests.hpp +++ b/tests/IntegrationTests/PhasorDynamics/PDIntegrationTests.hpp @@ -4,7 +4,14 @@ #include #include #include +#include +#include +#include +#include +#include +#include #include +#include #include #include #include @@ -716,6 +723,222 @@ namespace GridKit auto success = compare(set_data, file_data); return success.report(__func__); } + + /// A finite plant-reference pulse moves the coupled REPCA, REECB, and + /// REGCA states, after which the closed loop returns to equilibrium. + TestOutcome regcaReecbRepca() + { + using namespace GridKit::PhasorDynamics::Controller; + using namespace GridKit::PhasorDynamics::Converter; + using ReecbVar = ReecbInternalVariables; + using RepcaVar = RepcaInternalVariables; + using RegcaVar = RegcaInternalVariables; + + constexpr IdxT RENEWABLE_BUS_ID = static_cast(23); + constexpr IdxT IPCMD_SIGNAL_ID = static_cast(201); + constexpr IdxT IQCMD_SIGNAL_ID = static_cast(202); + constexpr IdxT IBRANCHR_SIGNAL_ID = static_cast(203); + constexpr IdxT IBRANCHI_SIGNAL_ID = static_cast(204); + constexpr IdxT PBRANCH_SIGNAL_ID = static_cast(205); + constexpr IdxT QBRANCH_SIGNAL_ID = static_cast(206); + constexpr IdxT QEXT_SIGNAL_ID = static_cast(207); + constexpr IdxT PEXT_SIGNAL_ID = static_cast(208); + constexpr IdxT PLANT_PREF_SIGNAL_ID = static_cast(209); + constexpr IdxT REGCA_COMPONENT_ID = static_cast(0); + constexpr IdxT REECB_COMPONENT_ID = static_cast(1); + constexpr IdxT REPCA_COMPONENT_ID = static_cast(2); + constexpr RealT COMPONENT_MVA = static_cast(50.0); + constexpr RealT INITIAL_ACTIVE_POWER = static_cast(0.4); + constexpr RealT INITIAL_REACTIVE_POWER = static_cast(0.05); + constexpr RealT REFERENCE_PULSE = static_cast(0.05); + constexpr RealT PULSE_END = static_cast(0.1); + constexpr RealT RESPONSE_TOLERANCE = static_cast(0.01); + constexpr RealT RECOVERY_HORIZON = static_cast(25.0); + constexpr RealT RECOVERY_MONITOR_STEP = static_cast(1.0 / 60.0); + constexpr RealT RECOVERY_TOLERANCE = static_cast(1.0e-6); + + TestStatus success = true; + + SystemModelDataT data; + data.va_base = static_cast(100.0e6); + + auto& bus = data.bus.emplace_back(); + bus.bus_id = RENEWABLE_BUS_ID; + bus.bus_type = BusDataT::BusType::SLACK; + bus.Vr0 = ONE; + bus.Vi0 = ZERO; + + data.signal = {{"Active Current Command", IPCMD_SIGNAL_ID}, + {"Reactive Current Command", IQCMD_SIGNAL_ID}, + {"Branch Current Real", IBRANCHR_SIGNAL_ID}, + {"Branch Current Imaginary", IBRANCHI_SIGNAL_ID}, + {"Branch Active Power", PBRANCH_SIGNAL_ID}, + {"Branch Reactive Power", QBRANCH_SIGNAL_ID}, + {"Reactive Power Command", QEXT_SIGNAL_ID}, + {"Active Power Command", PEXT_SIGNAL_ID}, + {"Plant Active Power Reference", PLANT_PREF_SIGNAL_ID}}; + + auto& converter = data.regca.emplace_back(); + converter.buses[RegcaBuses::bus] = RENEWABLE_BUS_ID; + converter.signal_inputs[RegcaSignalInputs::ipcmd] = IPCMD_SIGNAL_ID; + converter.signal_inputs[RegcaSignalInputs::iqcmd] = IQCMD_SIGNAL_ID; + converter.signal_outputs[RegcaSignalOutputs::ibranchr] = IBRANCHR_SIGNAL_ID; + converter.signal_outputs[RegcaSignalOutputs::ibranchi] = IBRANCHI_SIGNAL_ID; + converter.signal_outputs[RegcaSignalOutputs::pbranch] = PBRANCH_SIGNAL_ID; + converter.signal_outputs[RegcaSignalOutputs::qbranch] = QBRANCH_SIGNAL_ID; + converter.parameters[RegcaParameters::p0] = INITIAL_ACTIVE_POWER; + converter.parameters[RegcaParameters::q0] = INITIAL_REACTIVE_POWER; + converter.parameters[RegcaParameters::mva] = COMPONENT_MVA; + converter.parameters[RegcaParameters::Tg] = static_cast(0.02); + converter.parameters[RegcaParameters::TM] = static_cast(0.02); + converter.parameters[RegcaParameters::Rqmax] = static_cast(999.0); + converter.parameters[RegcaParameters::Rqmin] = static_cast(-999.0); + converter.parameters[RegcaParameters::Rpmax] = static_cast(999.0); + converter.parameters[RegcaParameters::sL] = true; + converter.parameters[RegcaParameters::IL1] = static_cast(1.1); + converter.parameters[RegcaParameters::VL0] = static_cast(0.4); + converter.parameters[RegcaParameters::VL1] = static_cast(0.9); + converter.parameters[RegcaParameters::VA0] = static_cast(0.4); + converter.parameters[RegcaParameters::VA1] = static_cast(0.9); + converter.parameters[RegcaParameters::Vhvmax] = static_cast(1.2); + + auto& controller = data.reecb.emplace_back(); + controller.buses[ReecbBuses::bus] = RENEWABLE_BUS_ID; + controller.signal_inputs[ReecbSignalInputs::pe] = PBRANCH_SIGNAL_ID; + controller.signal_inputs[ReecbSignalInputs::qgen] = QBRANCH_SIGNAL_ID; + controller.signal_inputs[ReecbSignalInputs::qext] = QEXT_SIGNAL_ID; + controller.signal_inputs[ReecbSignalInputs::pref] = PEXT_SIGNAL_ID; + controller.signal_outputs[ReecbSignalOutputs::ipcmd] = IPCMD_SIGNAL_ID; + controller.signal_outputs[ReecbSignalOutputs::iqcmd] = IQCMD_SIGNAL_ID; + controller.parameters[ReecbParameters::mva] = COMPONENT_MVA; + controller.parameters[ReecbParameters::Trv] = static_cast(0.02); + controller.parameters[ReecbParameters::Tp] = static_cast(0.02); + controller.parameters[ReecbParameters::Kvi] = static_cast(5.0); + controller.parameters[ReecbParameters::QFlag] = true; + controller.parameters[ReecbParameters::VFlag] = true; + + auto& plant = data.repca.emplace_back(); + plant.buses[RepcaBuses::bus] = RENEWABLE_BUS_ID; + plant.signal_inputs[RepcaSignalInputs::ir] = IBRANCHR_SIGNAL_ID; + plant.signal_inputs[RepcaSignalInputs::ii] = IBRANCHI_SIGNAL_ID; + plant.signal_inputs[RepcaSignalInputs::p] = PBRANCH_SIGNAL_ID; + plant.signal_inputs[RepcaSignalInputs::q] = QBRANCH_SIGNAL_ID; + plant.signal_inputs[RepcaSignalInputs::pref] = PLANT_PREF_SIGNAL_ID; + plant.signal_outputs[RepcaSignalOutputs::qext] = QEXT_SIGNAL_ID; + plant.signal_outputs[RepcaSignalOutputs::pext] = PEXT_SIGNAL_ID; + plant.parameters[RepcaParameters::mva] = COMPONENT_MVA; + plant.parameters[RepcaParameters::Freqflag] = true; + plant.parameters[RepcaParameters::Ddn] = ZERO; + plant.parameters[RepcaParameters::Dup] = ZERO; + plant.parameters[RepcaParameters::Tp] = static_cast(0.02); + plant.parameters[RepcaParameters::Tlag] = static_cast(0.5); + + auto& reference = data.constant_source.emplace_back(); + reference.parameters[ConstantSignalSourceParameters::Sr] = INITIAL_ACTIVE_POWER; + reference.signal_outputs[ConstantSignalSourceSignalOutputs::sr] = PLANT_PREF_SIGNAL_ID; + + SystemModel system(data); + success *= system.allocate() == 0; + + auto* regca = + dynamic_cast*>(system.getComponent(REGCA_COMPONENT_ID)); + auto* reecb = + dynamic_cast*>(system.getComponent(REECB_COMPONENT_ID)); + auto* repca = + dynamic_cast*>(system.getComponent(REPCA_COMPONENT_ID)); + if (regca == nullptr || reecb == nullptr || repca == nullptr) + { + success = false; + return success.report(__func__); + } + + AnalysisManager::Sundials::Ida ida(&system); + success *= ida.configureSimulation() == 0; + success *= ida.initializeSimulation(ZERO) == 0; + + const auto pord_index = static_cast( + reecb->getVariableIndex(static_cast(ReecbVar::PORD))); + const auto ipcmd_index = static_cast( + reecb->getVariableIndex(static_cast(ReecbVar::IPCMD))); + const auto regca_ip_index = static_cast( + regca->getVariableIndex(static_cast(RegcaVar::IP))); + const auto pext_index = static_cast( + repca->getVariableIndex(static_cast(RepcaVar::PEXT))); + + const auto* equilibrium_values = system.y().getData(); + const std::vector equilibrium( + equilibrium_values, + equilibrium_values + static_cast(system.y().getSize())); + + success *= isEqual(static_cast(system.getSignal(PLANT_PREF_SIGNAL_ID)->read()), + INITIAL_ACTIVE_POWER, + RECOVERY_TOLERANCE); + success *= isEqual(static_cast(system.getSignal(PEXT_SIGNAL_ID)->read()), + INITIAL_ACTIVE_POWER, + RECOVERY_TOLERANCE); + success *= isEqual(static_cast(system.getSignal(QEXT_SIGNAL_ID)->read()), + INITIAL_REACTIVE_POWER, + RECOVERY_TOLERANCE); + + auto* pref_signal = system.getSignal(PLANT_PREF_SIGNAL_ID); + const RealT pref0 = static_cast(pref_signal->read()); + + pref_signal->init(pref0 + REFERENCE_PULSE); + success *= ida.initializeSimulation(ZERO) == 0; + success *= ida.runSimulation(PULSE_END, RECOVERY_MONITOR_STEP) == 0; + + const auto* pulse_values = system.y().getData(); + const RealT pord_response = pulse_values[pord_index] - equilibrium[pord_index]; + const RealT ipcmd_response = pulse_values[ipcmd_index] - equilibrium[ipcmd_index]; + const RealT regca_ip_response = pulse_values[regca_ip_index] - equilibrium[regca_ip_index]; + const RealT pext_response = pulse_values[pext_index] - equilibrium[pext_index]; + + if (pext_response <= RESPONSE_TOLERANCE) + { + std::cout << "REPCA PEXT responded by only " << pext_response + << " during the reference pulse\n"; + success = false; + } + + if (pord_response <= RESPONSE_TOLERANCE) + { + std::cout << "REECB PORD responded by only " << pord_response + << " during the reference pulse\n"; + success = false; + } + if (ipcmd_response <= RESPONSE_TOLERANCE) + { + std::cout << "REECB IPCMD responded by only " << ipcmd_response + << " during the reference pulse\n"; + success = false; + } + if (regca_ip_response <= RESPONSE_TOLERANCE) + { + std::cout << "REGCA IP responded by only " << regca_ip_response + << " during the reference pulse\n"; + success = false; + } + + pref_signal->init(pref0); + success *= ida.initializeSimulation(PULSE_END) == 0; + success *= ida.runSimulation(PULSE_END + RECOVERY_HORIZON, + RECOVERY_MONITOR_STEP) + == 0; + + const auto* final_values = system.y().getData(); + for (size_t entry = 0; entry < equilibrium.size(); ++entry) + { + const RealT deviation = final_values[entry] - equilibrium[entry]; + if (!isEqual(deviation, ZERO, RECOVERY_TOLERANCE)) + { + std::cout << "State " << entry << " remains " << deviation + << " from its equilibrium after recovery\n"; + success = false; + } + } + + return success.report(__func__); + } }; } // namespace Testing } // namespace GridKit diff --git a/tests/IntegrationTests/PhasorDynamics/runPDIntegrationTests.cpp b/tests/IntegrationTests/PhasorDynamics/runPDIntegrationTests.cpp index 1f4c1fe36..5be32152b 100644 --- a/tests/IntegrationTests/PhasorDynamics/runPDIntegrationTests.cpp +++ b/tests/IntegrationTests/PhasorDynamics/runPDIntegrationTests.cpp @@ -12,6 +12,7 @@ int main() result += test.twoBusTgov1(); result += test.threeBusBasic(); result += test.threeBusClassical(); + result += test.regcaReecbRepca(); return result.summary(); } diff --git a/tests/UnitTests/PhasorDynamics/CMakeLists.txt b/tests/UnitTests/PhasorDynamics/CMakeLists.txt index 6c8721977..af3ec6f23 100644 --- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt +++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt @@ -132,6 +132,16 @@ target_link_libraries( GridKit::phasor_dynamics_bus_dependency_tracking GridKit::testing) +add_executable(test_phasor_controller_reecb runControllerReecbTests.cpp) +target_link_libraries( + test_phasor_controller_reecb + GridKit::definitions + GridKit::phasor_dynamics_controller_reecb + GridKit::phasor_dynamics_controller_reecb_dependency_tracking + GridKit::phasor_dynamics_bus + GridKit::phasor_dynamics_bus_dependency_tracking + GridKit::testing) + add_executable(test_phasor_controller_repca runControllerRepcaTests.cpp) target_link_libraries( test_phasor_controller_repca @@ -197,6 +207,7 @@ add_test(NAME PhasorDynamicsExciterEsdc1aTest COMMAND test_phasor_exciter_esdc1a 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 PhasorDynamicsControllerReecbTest COMMAND test_phasor_controller_reecb) add_test(NAME PhasorDynamicsControllerRepcaTest COMMAND test_phasor_controller_repca) add_test(NAME PhasorDynamicsStabilizerIeeestTest COMMAND test_phasor_stabilizer_ieeest) add_test(NAME PhasorDynamicsGenClassicalTest COMMAND test_phasor_gen_classical) @@ -224,6 +235,7 @@ install( test_phasor_gensal test_phasor_exciter_sexspti test_phasor_converter_regca + test_phasor_controller_reecb test_phasor_controller_repca test_phasor_stabilizer_ieeest test_phasor_gen_classical diff --git a/tests/UnitTests/PhasorDynamics/ComponentConnectionTests.hpp b/tests/UnitTests/PhasorDynamics/ComponentConnectionTests.hpp index a849c5c4e..920805837 100644 --- a/tests/UnitTests/PhasorDynamics/ComponentConnectionTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ComponentConnectionTests.hpp @@ -4,6 +4,8 @@ #include #include +#include +#include #include #include #include @@ -228,6 +230,102 @@ namespace GridKit return success.report(__func__); } + + /// REGCA initializes first and publishes the current commands it + /// resolves to the shared nodes, alongside its branch powers. REECB + /// then initializes around all four published values and must leave + /// them unchanged at a steady state. + TestOutcome regcaReecb() + { + using ConverterExternal = PhasorDynamics::Converter::RegcaExternalVariables; + using ConverterInternal = PhasorDynamics::Converter::RegcaInternalVariables; + using ConverterParams = PhasorDynamics::Converter::RegcaParameters; + using ControllerExternal = PhasorDynamics::Controller::ReecbExternalVariables; + using ControllerInternal = PhasorDynamics::Controller::ReecbInternalVariables; + using ControllerParams = PhasorDynamics::Controller::ReecbParameters; + + TestStatus success = true; + + PhasorDynamics::SystemModel system; + PhasorDynamics::BusInfinite bus( + static_cast(1.0), + static_cast(0.0)); + PhasorDynamics::SignalNode ipcmd; + PhasorDynamics::SignalNode iqcmd; + PhasorDynamics::SignalNode pe; + PhasorDynamics::SignalNode qgen; + + // The operating point is exactly representable, so the pair rests at + // an exact steady state. + PhasorDynamics::Converter::RegcaData converter_data; + converter_data.parameters[ConverterParams::p0] = static_cast(0.375); + converter_data.parameters[ConverterParams::q0] = static_cast(0.0625); + converter_data.parameters[ConverterParams::mva] = static_cast(100.0); + converter_data.parameters[ConverterParams::Tg] = static_cast(0.02); + converter_data.parameters[ConverterParams::TM] = static_cast(0.02); + converter_data.parameters[ConverterParams::Rqmax] = static_cast(999.0); + converter_data.parameters[ConverterParams::Rqmin] = static_cast(-999.0); + converter_data.parameters[ConverterParams::Rpmax] = static_cast(999.0); + converter_data.parameters[ConverterParams::sL] = true; + converter_data.parameters[ConverterParams::IL1] = static_cast(1.1); + converter_data.parameters[ConverterParams::VL0] = static_cast(0.25); + converter_data.parameters[ConverterParams::VL1] = static_cast(0.75); + converter_data.parameters[ConverterParams::VA0] = static_cast(0.25); + converter_data.parameters[ConverterParams::VA1] = static_cast(0.75); + converter_data.parameters[ConverterParams::Vhvmax] = static_cast(1.5); + + PhasorDynamics::Converter::Regca converter(&bus, converter_data); + + PhasorDynamics::Controller::ReecbData controller_data; + controller_data.parameters[ControllerParams::mva] = static_cast(100.0); + controller_data.parameters[ControllerParams::Tp] = static_cast(0.02); + controller_data.parameters[ControllerParams::QFlag] = true; + controller_data.parameters[ControllerParams::VFlag] = true; + controller_data.parameters[ControllerParams::Pqflag] = true; + controller_data.parameters[ControllerParams::Imax] = static_cast(0.625); + controller_data.parameters[ControllerParams::Vmin] = static_cast(0.5); + controller_data.parameters[ControllerParams::Vmax] = static_cast(1.5); + + PhasorDynamics::Controller::Reecb controller(&bus, controller_data); + + controller.getSignals().template assignSignalNode(&ipcmd); + controller.getSignals().template assignSignalNode(&iqcmd); + converter.getSignals().template attachSignalNode(&ipcmd); + converter.getSignals().template attachSignalNode(&iqcmd); + converter.getSignals().template assignSignalNode(&pe); + converter.getSignals().template assignSignalNode(&qgen); + controller.getSignals().template attachSignalNode(&pe); + controller.getSignals().template attachSignalNode(&qgen); + + system.addBus(&bus); + system.addComponent(&converter); + system.addComponent(&controller); + + success *= system.allocate() == 0; + success *= ipcmd.linked() && iqcmd.linked() && pe.linked() && qgen.linked(); + success *= ipcmd.getVariableIndex() + == controller.getVariableIndex( + static_cast(ControllerInternal::IPCMD)); + success *= iqcmd.getVariableIndex() + == controller.getVariableIndex( + static_cast(ControllerInternal::IQCMD)); + success *= system.initialize() == 0; + success *= system.evaluateResidual() == 0; + + // At unit terminal voltage the shared nodes carry the scheduled powers. + success *= isEqual(ipcmd.read(), static_cast(0.375), kTol); + success *= isEqual(iqcmd.read(), static_cast(0.0625), kTol); + success *= isEqual(pe.read(), static_cast(0.375), kTol); + success *= isEqual(qgen.read(), static_cast(0.0625), kTol); + + const auto* residual = controller.getResidual().getData(); + for (IdxT row = 0; row < controller.size(); ++row) + { + success *= isEqual(residual[row], static_cast(0.0), kTol); + } + + return success.report(__func__); + } }; } // namespace Testing diff --git a/tests/UnitTests/PhasorDynamics/ControllerReecbTests.hpp b/tests/UnitTests/PhasorDynamics/ControllerReecbTests.hpp new file mode 100644 index 000000000..7e2752a4d --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/ControllerReecbTests.hpp @@ -0,0 +1,2681 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + template + class ControllerReecbTests + { + public: + using ScalarT = scalar_type; + using IdxT = index_type; + using RealT = typename PhasorDynamics::Component::RealT; + + ControllerReecbTests() = default; + ~ControllerReecbTests() = default; + + static constexpr RealT kTol = + static_cast(100.0) * std::numeric_limits::epsilon(); + + /// Validate construction, row layout, defaults, parameters, buses, + /// signal links, and the time-constant floor. + TestOutcome validation() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + + PhasorDynamics::Controller::Reecb empty(&bus); + success *= (empty.size() == static_cast(index(Vars::MAXIMUM))); + success *= (empty.getMonitor() == nullptr); + + const std::array row_order{{ + Vars::VMEAS, + Vars::PMEAS, + Vars::XPIQ, + Vars::XPIV, + Vars::QV, + Vars::PORD, + Vars::VT, + Vars::ILMAX, + Vars::IQCMD, + Vars::IPCMD, + }}; + for (size_t row = 0; row < row_order.size(); ++row) + { + success *= (index(row_order[row]) == row); + } + + Fixture configured(makeData()); + success *= (configured.reecb.size() == static_cast(index(Vars::MAXIMUM))); + success *= (configured.reecb.getMonitor() != nullptr); + success *= (configured.reecb.verify() == 0); + success *= (configured.reecb.initialize() != 0); + success *= (configured.reecb.allocate() == 0); + success *= (configured.reecb.tagDifferentiable() == 0); + success *= (static_cast(configured.reecb.getResidual().getSize()) + == index(Vars::MAXIMUM)); + + for (size_t row = 0; row < index(Vars::MAXIMUM); ++row) + { + const bool expected = row <= index(Vars::PORD); + if (configured.reecb.tag()[row] != expected) + { + std::cout << "REECB differentiability tag " << row << " mismatch\n"; + success = false; + } + } + + Fixture documented_defaults(makeMinimalData()); + success *= (documented_defaults.reecb.verify() == 0); + success *= defaultsMatchDocumentedValues(); + + // Integer JSON values are accepted for real parameters; booleans are + // not numeric. + auto integer_numeric = makeData(); + integer_numeric.parameters[Params::mva] = static_cast(100); + Fixture integer_parameter(integer_numeric); + success *= (integer_parameter.reecb.verify() == 0); + success *= invalidParameterCase(Params::mva, true); + + const RealT nan = std::numeric_limits::quiet_NaN(); + const RealT infinity = std::numeric_limits::infinity(); + const std::array real_parameters{{ + Params::mva, + Params::Trv, + Params::Tp, + Params::Vref0, + Params::Vdip, + Params::Vup, + Params::dbd1, + Params::dbd2, + Params::kqv, + Params::Iql1, + Params::Iqh1, + Params::Qmax, + Params::Qmin, + Params::Kqp, + Params::Kqi, + Params::Vmax, + Params::Vmin, + Params::Kvp, + Params::Kvi, + Params::Tiq, + Params::Tpord, + Params::dPmax, + Params::dPmin, + Params::Pmax, + Params::Pmin, + Params::Imax, + }}; + for (const Params parameter : real_parameters) + { + success *= invalidParameterCase(parameter, nan); + success *= invalidParameterCase(parameter, infinity); + success *= invalidParameterCase(parameter, -infinity); + } + + success *= invalidParameterCase(Params::mva, 0.0); + success *= invalidParameterCase(Params::Trv, -0.1); + success *= invalidParameterCase(Params::Vdip, 2.0); + success *= invalidParameterCase(Params::dbd1, 0.1); + success *= invalidParameterCase(Params::dbd2, -0.1); + success *= invalidParameterCase(Params::Iql1, 2.0); + success *= invalidParameterCase(Params::Qmin, 3.0); + success *= invalidParameterCase(Params::Vmin, 2.0); + success *= invalidParameterCase(Params::dPmin, 0.0); + success *= invalidParameterCase(Params::dPmax, 0.0); + success *= invalidParameterCase(Params::Pmin, 3.0); + success *= invalidParameterCase(Params::Imax, 0.0); + + const std::array nonnegative_gains{{ + Params::kqv, + Params::Kqp, + Params::Kqi, + Params::Kvp, + Params::Kvi, + }}; + for (const Params gain : nonnegative_gains) + { + success *= invalidParameterCase(gain, -0.1); + } + + const std::array flag_parameters{{ + Params::PfFlag, + Params::VFlag, + Params::QFlag, + Params::Pqflag, + }}; + const std::array valid_flag_values{{false, true}}; + const std::array invalid_integral_flag_values{{ + static_cast(0), + static_cast(1), + static_cast(2), + }}; + const std::array invalid_real_flag_values{{ + static_cast(0.0), + static_cast(0.5), + static_cast(1.0), + nan, + infinity, + }}; + for (const Params flag : flag_parameters) + { + for (const bool value : valid_flag_values) + { + auto data = makeData(); + data.parameters[flag] = value; + Fixture model(data); + success *= (model.reecb.verify() == 0); + } + + for (const IdxT value : invalid_integral_flag_values) + { + success *= invalidParameterCase(flag, value); + } + + for (const RealT value : invalid_real_flag_values) + { + success *= invalidParameterCase(flag, value); + } + } + + PhasorDynamics::Controller::Reecb busless(nullptr, makeData()); + busless.setSystemBase(kNominalFrequency, kSystemBaseVa); + success *= (busless.verify() > 0); + + success *= unlinkedSignalRejected(); + success *= unlinkedSignalRejected(); + success *= unlinkedSignalRejected(); + success *= unlinkedSignalRejected(); + success *= unlinkedSignalRejected(); + + auto floor_data = makeData(); + floor_data.parameters[Params::Trv] = 0.0; + floor_data.parameters[Params::Tp] = 0.0; + floor_data.parameters[Params::Tiq] = 0.0; + floor_data.parameters[Params::Tpord] = 0.0; + + Fixture floored(floor_data); + success *= floored.initialize(kInitialIqcmd, kInitialIpcmd); + success *= (floored.evaluate() == 0); + success *= allResidualsWithinInitTolerance(floored.reecb); + + // Each floored lag turns a half-unit state offset into a rate of 500, + // and saturates the active-power ramp limiter. + setState(floored.reecb, {{Vars::VMEAS, 0.5}}); + success *= (floored.evaluate() == 0); + success *= residualsMatch(floored.reecb, + {{Vars::VMEAS, 500.0}}, + "floored voltage filter"); + + setState(floored.reecb, + {{Vars::VMEAS, 1.0}, + {Vars::PMEAS, 1.0}, + {Vars::QV, 1.0}, + {Vars::PORD, 1.0}}); + success *= (floored.evaluate() == 0); + success *= residualsMatch(floored.reecb, + {{Vars::PMEAS, 500.0}, + {Vars::QV, 500.0}, + {Vars::PORD, 1.0}}, + "floored time constants"); + + return success.report(__func__); + } + + /// Check initialization state, known-input preservation, unknown-reference + /// publication, command aliases, latches, monitors, and the power base. + TestOutcome initializationAndSignals() + { + TestStatus success = true; + + Fixture fixture(makeData(), 0.8, 0.6); + fixture.attachAllInputs(99.0); + fixture.input(Ext::PE) = kInitialIpcmd; + fixture.input(Ext::QGEN) = kInitialIqcmd; + success *= fixture.initialize(kInitialIqcmd, kInitialIpcmd); + success *= (fixture.evaluate() == 0); + + const std::array initial_state{{ + {Vars::VMEAS, 1.0}, + {Vars::PMEAS, 1.5}, + {Vars::XPIQ, 0.0}, + {Vars::XPIV, 0.0}, + {Vars::QV, 1.5}, + {Vars::PORD, 1.5}, + {Vars::VT, 1.0}, + {Vars::ILMAX, 2.0}, + }}; + success *= stateMatches(fixture.reecb, initial_state, "initialization"); + + success *= scalarPreserved(fixture.iqcmd(), kInitialIqcmd, "preserved iqcmd"); + success *= scalarPreserved(fixture.ipcmd(), kInitialIpcmd, "preserved ipcmd"); + success *= scalarPreserved(fixture.input(Ext::PE), kInitialIpcmd, "preserved pe"); + success *= scalarPreserved(fixture.input(Ext::QGEN), kInitialIqcmd, "preserved qgen"); + success *= scalarMatches(fixture.input(Ext::QEXT), 0.75, "published qext"); + success *= scalarMatches(fixture.input(Ext::PFAREF), 0.0, "published pfaref"); + success *= scalarMatches(fixture.input(Ext::PREF), 0.75, "published pref"); + success *= allResidualsWithinInitTolerance(fixture.reecb); + + success *= monitorMatches(fixture.reecb, + {{kInitialIqcmd, kInitialIpcmd, 1.0, 1.5}}, + "initialization"); + + constexpr RealT absolute_tolerance = 2.5e-7; + success *= (fixture.reecb.setAbsoluteTolerance(absolute_tolerance) == 0); + const auto* tolerances = fixture.reecb.absoluteTolerance().getData(); + for (size_t row = 0; row < index(Vars::MAXIMUM); ++row) + { + success *= valueUnchanged(tolerances[row], absolute_tolerance, "absolute tolerance", row); + } + + // Unassigned command outputs keep the commands in the model vector. + Fixture latched(makeData(), 1.0, 0.0, kSystemBaseVa, false); + success *= latched.initialize(kInitialIqcmd, kInitialIpcmd); + success *= (latched.evaluate() == 0); + success *= scalarPreserved(latched.iqcmd(), kInitialIqcmd, "unassigned iqcmd"); + success *= scalarPreserved(latched.ipcmd(), kInitialIpcmd, "unassigned ipcmd"); + success *= allResidualsWithinInitTolerance(latched.reecb); + + // An omitted component rating falls back to the system power base, so + // the same commands land on a different measured power. + auto system_base_data = makeData(); + system_base_data.parameters.erase(Params::mva); + Fixture system_base(system_base_data, 1.0, 0.0, static_cast(50.0e6)); + system_base.attachAllInputs(); + system_base.input(Ext::PE) = 0.75; + success *= system_base.initialize(kInitialIqcmd, 1.5); + success *= (system_base.evaluate() == 0); + success *= stateMatches(system_base.reecb, + {{Vars::PMEAS, 0.75}, + {Vars::PORD, 1.5}, + {Vars::ILMAX, 2.0}}, + "omitted component rating"); + success *= allResidualsWithinInitTolerance(system_base.reecb); + + return success.report(__func__); + } + + /// Check adjusted limits, initialization rejection, and atomicity. + TestOutcome initializationDomain() + { + TestStatus success = true; + + const auto data = makeData(); + + success *= initializationRejectedAtomically(data, 0.75, -0.1, "negative active-current command"); + success *= initializationRejectedAtomically(data, 0.75, std::numeric_limits::infinity(), "nonfinite active-current command"); + success *= initializationRejectedAtomically( + data, std::numeric_limits::infinity(), 0.75, "nonfinite reactive-current command"); + + auto pord_above = data; + pord_above.parameters[Params::Pmax] = 1.0; + Fixture adjusted_pmax(pord_above); + success *= adjusted_pmax.initialize(0.75, 0.75); + success *= (adjusted_pmax.evaluate() == 0); + success *= stateMatches(adjusted_pmax.reecb, {{Vars::PORD, 1.5}}, "adjusted Pmax"); + success *= allResidualsWithinInitTolerance(adjusted_pmax.reecb); + setState(adjusted_pmax.reecb, {{Vars::PORD, 1.25}}); + success *= (adjusted_pmax.evaluate() == 0); + success *= residualsMatch(adjusted_pmax.reecb, {{Vars::PORD, 1.0}}, "adjusted Pmax"); + + auto pord_below = data; + pord_below.parameters[Params::Pmin] = 2.0; + Fixture adjusted_pmin(pord_below); + success *= adjusted_pmin.initialize(0.75, 0.75); + success *= (adjusted_pmin.evaluate() == 0); + success *= stateMatches(adjusted_pmin.reecb, {{Vars::PORD, 1.5}}, "adjusted Pmin"); + success *= allResidualsWithinInitTolerance(adjusted_pmin.reecb); + setState(adjusted_pmin.reecb, {{Vars::PORD, 1.75}}); + success *= (adjusted_pmin.evaluate() == 0); + success *= residualsMatch(adjusted_pmin.reecb, {{Vars::PORD, -1.0}}, "adjusted Pmin"); + + auto expanded_current = data; + expanded_current.parameters[Params::Imax] = 1.0; + Fixture adjusted_imax(expanded_current); + success *= adjusted_imax.initialize(0.75, 0.75); + success *= (adjusted_imax.evaluate() == 0); + success *= stateMatches(adjusted_imax.reecb, {{Vars::ILMAX, 1.5}}, "adjusted Imax"); + success *= allResidualsWithinInitTolerance(adjusted_imax.reecb); + + auto reactive_pi = data; + reactive_pi.parameters[Params::QFlag] = true; + reactive_pi.parameters[Params::VFlag] = true; + reactive_pi.parameters[Params::Kqi] = 5.0; + const std::array q_limits{{-1.25, 1.25}}; + for (const RealT qgen : q_limits) + { + Fixture adjusted_q(reactive_pi); + adjusted_q.attachAllInputs(); + adjusted_q.input(Ext::PE) = 0.75; + adjusted_q.input(Ext::QGEN) = qgen; + success *= adjusted_q.initialize(0.75, 0.75); + success *= (adjusted_q.evaluate() == 0); + success *= allResidualsWithinInitTolerance(adjusted_q.reecb); + } + + auto voltage_pi = reactive_pi; + voltage_pi.parameters[Params::Kqi] = 0.0; + voltage_pi.parameters[Params::Kvi] = 5.0; + + struct VoltageLimitCase + { + RealT voltage; + RealT pe; + RealT qgen; + }; + + const std::array voltage_limits{{ + {0.3, 0.3, 0.3}, + {1.6, 0.96, 0.8}, + }}; + for (const auto& test_case : voltage_limits) + { + Fixture adjusted_v(voltage_pi, test_case.voltage); + adjusted_v.attachAllInputs(); + adjusted_v.input(Ext::PE) = test_case.pe; + adjusted_v.input(Ext::QGEN) = test_case.qgen; + success *= adjusted_v.initialize(0.75, 0.75); + success *= (adjusted_v.evaluate() == 0); + success *= allResidualsWithinInitTolerance(adjusted_v.reecb); + } + + // Power-factor control needs a representable angle. + auto power_factor = data; + power_factor.parameters[Params::PfFlag] = true; + + auto late_data = power_factor; + late_data.parameters[Params::Pmax] = 1.0; + Fixture late(late_data); + late.attachAllInputs(); + late.input(Ext::PE) = 0.0; + success *= late.prepare(0.75, 0.75); + if (late.reecb.initialize() == 0) + { + std::cout << "Expected REECB initialization rejection after Pmax adjustment\n"; + success = false; + } + late.input(Ext::PREF) = 0.75; + setState(late.reecb, {{Vars::PORD, 1.25}, {Vars::VT, 1.0}}); + setDerivative(late.reecb, {{Vars::PORD, 0.0}}); + success *= (late.evaluate() == 0); + success *= residualsMatch(late.reecb, {{Vars::PORD, 0.0}}, "rejected Pmax adjustment"); + success *= initializationRejectedAtomically( + power_factor, 0.75, 0.75, "unrepresentable power-factor reference", 1.0e-8, 0.75); + + success *= initializationRejectedAtomically(data, 0.75, 0.75, "zero terminal voltage", 0.75, 0.75, 0.0); + success *= initializationRejectedAtomically( + data, 0.75, 0.75, "nonfinite active-power feedback", std::numeric_limits::infinity(), 0.75); + + // Collapsed limits remain pinned at their initial output or expand to + // include it. + auto collapsed = reactive_pi; + collapsed.parameters[Params::Kvi] = 0.5; + collapsed.parameters[Params::Qmin] = 1.2; + collapsed.parameters[Params::Qmax] = 1.2; + collapsed.parameters[Params::Vmin] = 1.0; + collapsed.parameters[Params::Vmax] = 1.0; + Fixture collapsed_limits(collapsed); + collapsed_limits.attachAllInputs(); + collapsed_limits.input(Ext::PE) = 0.6; + collapsed_limits.input(Ext::QGEN) = 0.6; + success *= collapsed_limits.initialize(0.75, 0.75); + success *= (collapsed_limits.evaluate() == 0); + success *= allResidualsWithinInitTolerance(collapsed_limits.reecb); + + auto collapsed_reactive = collapsed; + collapsed_reactive.parameters[Params::Vmin] = 0.5; + collapsed_reactive.parameters[Params::Vmax] = 1.5; + collapsed_reactive.parameters[Params::Kvi] = 0.0; + Fixture expanded_reactive(collapsed_reactive); + expanded_reactive.attachAllInputs(); + expanded_reactive.input(Ext::PE) = 0.75; + expanded_reactive.input(Ext::QGEN) = 0.3; + success *= expanded_reactive.initialize(0.75, 0.75); + success *= (expanded_reactive.evaluate() == 0); + success *= allResidualsWithinInitTolerance(expanded_reactive.reecb); + + auto collapsed_voltage = collapsed; + collapsed_voltage.parameters[Params::Qmin] = -2.0; + collapsed_voltage.parameters[Params::Qmax] = 2.0; + collapsed_voltage.parameters[Params::Kqi] = 0.0; + collapsed_voltage.parameters[Params::Vmin] = 1.4; + collapsed_voltage.parameters[Params::Vmax] = 1.4; + Fixture expanded_voltage(collapsed_voltage); + expanded_voltage.attachAllInputs(); + expanded_voltage.input(Ext::PE) = 0.75; + expanded_voltage.input(Ext::QGEN) = 0.75; + success *= expanded_voltage.initialize(0.75, 0.75); + success *= (expanded_voltage.evaluate() == 0); + success *= allResidualsWithinInitTolerance(expanded_voltage.reecb); + + // Zero integral gains leave both controllers unconstrained, so any + // reactive feedback initializes. + auto zero_gains = reactive_pi; + zero_gains.parameters[Params::Kqi] = 0.0; + zero_gains.parameters[Params::Kvi] = 0.0; + Fixture unconstrained(zero_gains); + unconstrained.attachAllInputs(); + unconstrained.input(Ext::PE) = 0.75; + unconstrained.input(Ext::QGEN) = 4.0; + success *= unconstrained.initialize(0.75, 0.75); + success *= (unconstrained.evaluate() == 0); + success *= allResidualsWithinInitTolerance(unconstrained.reecb); + + // An invalid configuration is rejected before any state is written. + auto invalid_data = data; + invalid_data.parameters[Params::Imax] = 0.0; + Fixture invalid_fixture(invalid_data); + invalid_fixture.attachAllInputs(); + success *= (invalid_fixture.reecb.allocate() == 0); + poisonState(invalid_fixture, 0.75, 0.75); + const auto invalid_y = copyVector(invalid_fixture.reecb.y()); + const auto invalid_yp = copyVector(invalid_fixture.reecb.yp()); + if (invalid_fixture.reecb.initialize() == 0) + { + std::cout << "Expected REECB initialization rejection: invalid configuration\n"; + success = false; + } + success *= vectorUnchanged(invalid_fixture.reecb.y(), invalid_y, "state"); + success *= vectorUnchanged(invalid_fixture.reecb.yp(), invalid_yp, "derivative"); + + return success.report(__func__); + } + + /// The smooth-limiter inverse reproduces interior and boundary commands, + /// including points that expand the current circle. + TestOutcome initializationExactness() + { + TestStatus success = true; + + struct ExactnessCase + { + RealT ipcmd; + const char* label; + }; + + const std::array active_cases{{ + {1.0e-6, "near the lower active-current limit"}, + {0.75, "interior active-current command"}, + {1.249999, "near the upper active-current limit"}, + }}; + + // The recovered order limits are widened so the reconstruction, not + // the order limit, decides admissibility at the command endpoints. + auto exactness_data = makeData(); + exactness_data.parameters[Params::Pmin] = -1.0; + exactness_data.parameters[Params::Pmax] = 3.0; + + for (const auto& test_case : active_cases) + { + Fixture fixture(exactness_data); + success *= fixture.initialize(0.0, test_case.ipcmd); + success *= (fixture.evaluate() == 0); + success *= scalarPreserved(fixture.ipcmd(), test_case.ipcmd, test_case.label); + success *= allResidualsWithinInitTolerance(fixture.reecb); + } + + // An interior command recovers the ideal active-power order exactly. + Fixture interior(makeData()); + success *= interior.initialize(0.0, 0.75); + success *= (interior.evaluate() == 0); + success *= stateMatches(interior.reecb, + {{Vars::PORD, 1.5}}, + "interior active-power order"); + + // A command whose recovered order lands exactly on the order limit is + // admitted rather than rejected. + auto limit_data = makeData(); + limit_data.parameters[Params::Pmax] = 1.5; + Fixture at_limit(limit_data); + success *= at_limit.initialize(0.0, 0.75); + success *= (at_limit.evaluate() == 0); + success *= stateMatches(at_limit.reecb, {{Vars::PORD, 1.5}}, "order at Pmax"); + success *= allResidualsWithinInitTolerance(at_limit.reecb); + + struct AsymmetricSlewCase + { + RealT minimum; + RealT maximum; + const char* label; + }; + + const std::array asymmetric_slews{{ + {-0.001, 0.1, "narrow negative ramp limit"}, + {-0.1, 0.001, "narrow positive ramp limit"}, + }}; + + for (const auto& test_case : asymmetric_slews) + { + auto asymmetric_data = makeData(); + asymmetric_data.parameters[Params::dPmin] = test_case.minimum; + asymmetric_data.parameters[Params::dPmax] = test_case.maximum; + + Fixture asymmetric(asymmetric_data); + asymmetric.attachAllInputs(); + success *= asymmetric.initialize(0.0, 0.75); + success *= (asymmetric.evaluate() == 0); + success *= stateMatches(asymmetric.reecb, + {{Vars::PORD, 1.5}}, + test_case.label); + success *= allResidualsWithinInitTolerance(asymmetric.reecb); + success *= scalarMatches(asymmetric.input(Ext::PREF), + 0.75, + test_case.label); + } + + // The reactive command shares the inverse, at both signs. + const std::array reactive_commands{{ + static_cast(0.999999), + static_cast(-0.999999), + }}; + for (const RealT iqcmd : reactive_commands) + { + Fixture reactive(exactness_data); + success *= reactive.initialize(iqcmd, 0.75); + success *= (reactive.evaluate() == 0); + success *= scalarPreserved(reactive.iqcmd(), iqcmd, "near-limit reactive command"); + success *= allResidualsWithinInitTolerance(reactive.reecb); + } + + struct BoundaryCase + { + bool p_priority; + RealT iqcmd; + RealT ipcmd; + RealT ilmax; + const char* label; + }; + + const std::array boundary_cases{{ + {true, 0.75, 0.0, 2.5, "zero active-current command"}, + {true, 0.0, 1.25, 0.0, "zero reactive-current capacity"}, + {true, 1.0, 0.75, 2.0, "upper reactive-current command"}, + {true, -1.0, 0.75, 2.0, "lower reactive-current command"}, + {false, 0.75, 1.0, 2.0, "upper active-current command"}, + {false, 1.25, 0.0, 0.0, "zero active-current capacity"}, + {true, 0.75, 1.5, 1.5, "expanded current circle"}, + }}; + for (const auto& test_case : boundary_cases) + { + auto boundary_data = exactness_data; + boundary_data.parameters[Params::Pqflag] = test_case.p_priority; + Fixture boundary(boundary_data); + success *= boundary.initialize(test_case.iqcmd, test_case.ipcmd); + success *= (boundary.evaluate() == 0); + success *= scalarPreserved(boundary.iqcmd(), test_case.iqcmd, test_case.label); + success *= scalarPreserved(boundary.ipcmd(), test_case.ipcmd, test_case.label); + success *= stateMatches(boundary.reecb, {{Vars::ILMAX, test_case.ilmax}}, test_case.label); + success *= allResidualsWithinInitTolerance(boundary.reecb); + } + + auto separated_data = exactness_data; + separated_data.parameters[Params::Imax] = 1.0; + Fixture separated(separated_data); + success *= separated.initialize(5.0e-13, 0.5); + success *= (separated.evaluate() == 0); + success *= scalarPreserved(separated.iqcmd(), 5.0e-13, "scale-separated current command"); + success *= allResidualsWithinInitTolerance(separated.reecb); + + // A low configured Imax requires representable bisection to preserve + // a strict low-priority command. + auto capacity_data = exactness_data; + capacity_data.parameters[Params::mva] = 100.0; + capacity_data.parameters[Params::Pqflag] = true; + capacity_data.parameters[Params::Imax] = 0.1; + + const std::array, 2> capacity_cases{{ + {1.7, 0.3}, + {1.5, 0.5}, + }}; + for (const auto& [iqcmd, ipcmd] : capacity_cases) + { + Fixture capacity_fixture(capacity_data); + success *= capacity_fixture.initialize(iqcmd, ipcmd); + success *= (capacity_fixture.evaluate() == 0); + success *= scalarPreserved(capacity_fixture.iqcmd(), iqcmd, "low-priority command"); + success *= scalarPreserved(capacity_fixture.ipcmd(), ipcmd, "high-priority command"); + const RealT ilmax = static_cast(capacity_fixture.reecb.y().getData()[index(Vars::ILMAX)]); + const RealT ilcap = ilmax * ilmax + / std::sqrt(ilmax * ilmax + ReecbT::INITIALIZATION_TOLERANCE); + if (ilcap < iqcmd) + { + std::cout << "REECB low-priority capacity does not include its initial command\n"; + success = false; + } + success *= allResidualsWithinInitTolerance(capacity_fixture.reecb); + } + + auto nested_data = exactness_data; + nested_data.parameters[Params::mva] = 100.0; + nested_data.parameters[Params::Pqflag] = true; + nested_data.parameters[Params::QFlag] = true; + nested_data.parameters[Params::VFlag] = false; + nested_data.parameters[Params::Imax] = 1.0; + Fixture nested(nested_data); + nested.attachAllInputs(); + nested.input(Ext::PE) = 0.6; + nested.input(Ext::QGEN) = 0.8; + success *= nested.initialize(0.8, 0.6); + success *= (nested.evaluate() == 0); + success *= scalarPreserved(nested.iqcmd(), 0.8, "nested-clamp reactive command"); + success *= scalarPreserved(nested.ipcmd(), 0.6, "nested-clamp active command"); + success *= allResidualsWithinInitTolerance(nested.reecb); + + auto exhausted_data = exactness_data; + exhausted_data.parameters[Params::QFlag] = true; + exhausted_data.parameters[Params::kqv] = 1.0; + exhausted_data.parameters[Params::Iql1] = -0.4; + exhausted_data.parameters[Params::Iqh1] = 1.2; + exhausted_data.parameters[Params::Vref0] = 2.2; + Fixture exhausted(exhausted_data); + exhausted.attachAllInputs(); + exhausted.input(Ext::PE) = 1.25; + success *= exhausted.initialize(0.0, 1.25); + success *= (exhausted.evaluate() == 0); + success *= scalarPreserved(exhausted.iqcmd(), 0.0, "exhausted reactive-current capacity"); + success *= stateMatches(exhausted.reecb, {{Vars::ILMAX, 0.0}}, "injection does not expand current circle"); + success *= allResidualsWithinInitTolerance(exhausted.reecb); + + return success.report(__func__); + } + + /// Check every residual row against an independent numerical answer key. + /// The expected values are literals, not a second implementation of REECB. + TestOutcome residualEquations() + { + TestStatus success = true; + + Fixture fixture(makeResidualData(), kStateVr, kStateVi); + fixture.attachAllInputs(); + setAnswerKeyInputs(fixture); + success *= fixture.prepare(0.25, 0.35); + setAnswerKeyState(fixture.reecb); + success *= (fixture.evaluate() == 0); + + const std::array expected_residuals{{ + {Vars::VMEAS, 0.99}, + {Vars::PMEAS, 0.145}, + {Vars::XPIQ, 0.21}, + {Vars::XPIV, 0.13}, + {Vars::QV, -0.05}, + {Vars::PORD, 0.26}, + {Vars::VT, -0.03}, + {Vars::ILMAX, 0.32}, + {Vars::IQCMD, 0.19}, + {Vars::IPCMD, 0.05}, + }}; + + success *= (static_cast(fixture.reecb.getResidual().getSize()) + == expected_residuals.size()); + for (size_t row = 0; row < expected_residuals.size(); ++row) + { + if (index(expected_residuals[row].variable) != row) + { + std::cout << "REECB residual key position " << row << " names row " + << variableName(expected_residuals[row].variable) << '\n'; + success = false; + } + } + success *= residualsMatch(fixture.reecb, + expected_residuals, + "independent numerical answer key"); + + return success.report(__func__); + } + + /// Every selector combination initializes attached and unattached signals + /// to a zero-residual state. + TestOutcome selectorConfigurations() + { + TestStatus success = true; + + const std::array selector_values{{false, true}}; + for (const bool pf : selector_values) + { + for (const bool voltage : selector_values) + { + for (const bool reactive : selector_values) + { + for (const bool p_priority : selector_values) + { + auto data = makeData(); + data.parameters[Params::PfFlag] = pf; + data.parameters[Params::VFlag] = voltage; + data.parameters[Params::QFlag] = reactive; + data.parameters[Params::Pqflag] = p_priority; + data.parameters[Params::Kqi] = reactive && voltage ? 0.4 : 0.0; + data.parameters[Params::Kvi] = reactive ? 0.5 : 0.0; + + for (const bool attached : selector_values) + { + Fixture fixture(data); + if (attached) + { + fixture.attachAllInputs(7.0); + fixture.input(Ext::PE) = 0.75; + fixture.input(Ext::QGEN) = 0.75; + } + + success *= fixture.initialize(0.75, 0.75); + success *= (fixture.evaluate() == 0); + success *= allResidualsWithinInitTolerance(fixture.reecb); + success *= scalarPreserved(fixture.iqcmd(), 0.75, "selector iqcmd"); + success *= scalarPreserved(fixture.ipcmd(), 0.75, "selector ipcmd"); + success *= stateMatches(fixture.reecb, {{Vars::ILMAX, 2.0}}, "selector ILMAX"); + + // Exactly one reactive path carries the operating point. + const auto* y = fixture.reecb.y().getData(); + if (reactive) + { + success *= (y[index(Vars::XPIV)] != ZERO); + if (voltage) + { + success *= (y[index(Vars::XPIQ)] != ZERO); + } + } + else + { + success *= (y[index(Vars::QV)] != ZERO); + } + + if (attached) + { + success *= scalarPreserved(fixture.input(Ext::PE), 0.75, "selector pe"); + success *= scalarPreserved(fixture.input(Ext::QGEN), 0.75, "selector qgen"); + const RealT expected_qext = reactive && !voltage ? 1.0 : (pf ? 0.0 : 0.75); + const RealT expected_pfaref = pf && (!reactive || voltage) ? kUnitSlopeAngle : 0.0; + success *= scalarMatches(fixture.input(Ext::QEXT), expected_qext, "published qext"); + success *= scalarMatches(fixture.input(Ext::PFAREF), expected_pfaref, "published pfaref"); + success *= scalarMatches(fixture.input(Ext::PREF), 0.75, "published pref"); + } + } + } + } + } + } + + return success.report(__func__); + } + + /// Direct-voltage mode consumes and publishes the Volt/VAr reference + /// without power-base conversion; the reactive modes convert it. + TestOutcome voltVarReferenceBase() + { + TestStatus success = true; + + { + auto data = makeData(); + data.parameters[Params::QFlag] = true; + data.parameters[Params::VFlag] = false; + data.parameters[Params::Kvi] = 0.5; + + Fixture fixture(data); + fixture.attachAllInputs(); + success *= fixture.initialize(0.75, 0.75); + success *= scalarMatches(fixture.input(Ext::QEXT), 1.0, "published voltage reference"); + success *= (fixture.evaluate() == 0); + success *= allResidualsWithinInitTolerance(fixture.reecb); + + // A raised external voltage reference enters the V-PI rate raw. + fixture.input(Ext::QEXT) = 1.2; + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.reecb, + {{Vars::XPIV, 0.1}}, + "unconverted voltage-reference rate"); + } + + { + Fixture fixture(makeData()); + fixture.attachAllInputs(); + success *= fixture.initialize(0.75, 0.75); + success *= scalarMatches(fixture.input(Ext::QEXT), 0.75, "published system-base reactive power"); + success *= (fixture.evaluate() == 0); + success *= allResidualsWithinInitTolerance(fixture.reecb); + + // The reactive-current lag keeps the power-base conversion, so the + // same raise produces twice the component-base rate. + fixture.input(Ext::QEXT) = 0.85; + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.reecb, + {{Vars::QV, 10.0}}, + "converted reactive-reference rate"); + } + + return success.report(__func__); + } + + /// Check the reactive selector paths, the voltage-band gate, the + /// reactive limits, both anti-windup gates, and the injection curve. + TestOutcome reactiveControl() + { + TestStatus success = true; + + { + // The constant-reactive path drives the current-command lag. + auto data = makeResidualData(); + data.parameters[Params::QFlag] = false; + data.parameters[Params::kqv] = 0.0; + Fixture fixture(data); + fixture.attachAllInputs(); + fixture.input(Ext::QEXT) = 0.4; + success *= fixture.prepare(0.0, 0.2); + setControlState(fixture.reecb); + setState(fixture.reecb, {{Vars::QV, 0.1}}); + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.reecb, {{Vars::QV, 1.4}}, "constant-reactive lag"); + + setState(fixture.reecb, {{Vars::VT, 0.0}}); + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.reecb, {{Vars::QV, 0.0}}, "gated constant-reactive lag"); + } + + { + // Cascaded Volt/VAr control runs the reactive PI and bypasses the lag. + auto data = makeResidualData(); + data.parameters[Params::QFlag] = true; + data.parameters[Params::VFlag] = true; + data.parameters[Params::kqv] = 0.0; + Fixture fixture(data); + fixture.attachAllInputs(); + fixture.input(Ext::QEXT) = 0.1; + fixture.input(Ext::QGEN) = -0.05; + success *= fixture.prepare(0.0, 0.2); + setControlState(fixture.reecb); + setState(fixture.reecb, {{Vars::XPIQ, 0.82}}); + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.reecb, + {{Vars::XPIQ, 0.12}, {Vars::QV, 0.0}}, + "reactive-power integral rate"); + + setState(fixture.reecb, {{Vars::VT, 0.0}}); + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.reecb, {{Vars::XPIQ, 0.0}}, "gated reactive-power integrator"); + } + + { + // The direct-voltage reference bypasses the reactive PI. + auto data = makeResidualData(); + data.parameters[Params::QFlag] = true; + data.parameters[Params::VFlag] = false; + data.parameters[Params::kqv] = 0.0; + Fixture fixture(data); + fixture.attachAllInputs(); + fixture.input(Ext::QEXT) = 1.05; + success *= fixture.prepare(0.0, 0.2); + setControlState(fixture.reecb); + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.reecb, + {{Vars::XPIQ, 0.0}, {Vars::XPIV, 0.025}}, + "voltage-control integral rate"); + + setState(fixture.reecb, {{Vars::VT, 2.0}}); + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.reecb, {{Vars::XPIV, 0.0}}, "gated voltage-control integrator"); + } + + { + // The reactive-power reference is limited before the error forms. + auto data = makeResidualData(); + data.parameters[Params::QFlag] = true; + data.parameters[Params::VFlag] = true; + data.parameters[Params::kqv] = 0.0; + data.parameters[Params::Kqp] = 0.0; + + const std::array reference_cases{{ + {-0.6, -0.28}, + {0.05, 0.04}, + {0.6, 0.32}, + }}; + for (const auto& test_case : reference_cases) + { + Fixture fixture(data); + fixture.attachAllInputs(); + fixture.input(Ext::QEXT) = test_case.input; + success *= fixture.prepare(0.0, 0.2); + setControlState(fixture.reecb); + setState(fixture.reecb, {{Vars::XPIQ, 1.0}}); + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.reecb, + {{Vars::XPIQ, test_case.expected}}, + "reactive-power reference limit"); + } + } + + { + // Saturated probes sit beyond their limit by a margin, so a blocked + // gate contributes nothing and an admitted gate passes the full rate. + auto data = makeResidualData(); + data.parameters[Params::QFlag] = true; + data.parameters[Params::VFlag] = true; + data.parameters[Params::kqv] = 0.0; + data.parameters[Params::Kqp] = 0.0; + data.parameters[Params::Qmin] = -3.0; + data.parameters[Params::Qmax] = 3.0; + + const std::array antiwindup_cases{{ + {2.2, 1.0, 0.0}, + {2.2, -1.0, -0.8}, + {-0.2, -1.0, 0.0}, + {-0.2, 1.0, 0.8}, + {1.0, 1.0, 0.8}, + }}; + for (const auto& test_case : antiwindup_cases) + { + Fixture fixture(data); + fixture.attachAllInputs(); + fixture.input(Ext::QEXT) = test_case.reference; + success *= fixture.prepare(0.0, 0.2); + setControlState(fixture.reecb); + setState(fixture.reecb, {{Vars::XPIQ, test_case.state}}); + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.reecb, + {{Vars::XPIQ, test_case.expected}}, + "reactive-power antiwindup"); + } + } + + { + // The voltage-control integrator saturates on the reactive-current + // limit carried by the current circle. + auto data = makeResidualData(); + data.parameters[Params::QFlag] = true; + data.parameters[Params::VFlag] = false; + data.parameters[Params::kqv] = 0.0; + data.parameters[Params::Kvp] = 0.0; + data.parameters[Params::Kvi] = 1.0; + + const std::array antiwindup_cases{{ + {1.0, 1.6, 0.0}, + {1.0, 0.4, -0.6}, + {-1.0, 0.4, 0.0}, + {-1.0, 1.6, 0.6}, + {0.0, 1.6, 0.6}, + }}; + for (const auto& test_case : antiwindup_cases) + { + Fixture fixture(data); + fixture.attachAllInputs(); + fixture.input(Ext::QEXT) = test_case.reference; + success *= fixture.prepare(0.0, 0.2); + setControlState(fixture.reecb); + setState(fixture.reecb, + {{Vars::XPIV, test_case.state}, {Vars::ILMAX, 0.5}}); + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.reecb, + {{Vars::XPIV, test_case.expected}}, + "voltage-control antiwindup"); + } + } + + { + // With the reactive lag at zero the command row reads the injection + // curve directly: a deadbanded voltage error scaled and limited. + auto data = makeResidualData(); + data.parameters[Params::QFlag] = false; + data.parameters[Params::kqv] = 1.0; + data.parameters[Params::dbd1] = -0.6; + data.parameters[Params::dbd2] = 0.6; + data.parameters[Params::Iql1] = -1.2; + data.parameters[Params::Iqh1] = 1.5; + data.parameters[Params::Vref0] = 3.0; + + const std::array injection_cases{{ + {5.5, -1.2}, + {4.2, -0.6}, + {3.0, 0.0}, + {1.8, 0.6}, + {0.4, 1.5}, + }}; + for (const auto& test_case : injection_cases) + { + Fixture fixture(data); + success *= fixture.prepare(0.0, 0.2); + setControlState(fixture.reecb); + setState(fixture.reecb, + {{Vars::VMEAS, test_case.input}, + {Vars::IQCMD, 0.0}, + {Vars::ILMAX, 3.0}}); + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.reecb, + {{Vars::IQCMD, test_case.expected}}, + "reactive-current injection"); + } + } + + { + // Power-factor control resolves the reactive reference from the + // measured active power and the commanded angle. + auto data = makeResidualData(); + data.parameters[Params::PfFlag] = true; + data.parameters[Params::QFlag] = false; + data.parameters[Params::kqv] = 0.0; + Fixture fixture(data); + fixture.attachAllInputs(); + fixture.input(Ext::PFAREF) = std::atan(HALF); + success *= fixture.prepare(0.0, 0.2); + setControlState(fixture.reecb); + setState(fixture.reecb, {{Vars::PMEAS, 0.6}, {Vars::QV, 0.1}}); + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.reecb, + {{Vars::QV, 0.4}}, + "power-factor reference"); + } + + return success.report(__func__); + } + + /// Check the active-power ramp, its voltage gate and anti-windup, both + /// command limits, the priority circle, and the signed continuation. + TestOutcome activeCurrentControl() + { + TestStatus success = true; + + { + // The ramp-rate limiter bounds the active-power order rate. + const std::array rate_cases{{ + {-1.0, -0.5}, + {0.2, 0.2}, + {1.0, 0.6}, + }}; + for (const auto& test_case : rate_cases) + { + Fixture fixture(makeResidualData()); + fixture.attachAllInputs(); + fixture.input(Ext::PREF) = rampReference(0.5, test_case.input); + success *= fixture.prepare(0.0, 0.2); + setControlState(fixture.reecb); + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.reecb, + {{Vars::PORD, test_case.expected}}, + "active-power ramp limit"); + } + } + + { + // Strongly asymmetric limits retain interior rates and bound each + // direction independently. + struct AsymmetricRateCase + { + RealT minimum; + RealT maximum; + RealT input; + RealT expected; + }; + + const std::array rate_cases{{ + {-0.001, 0.1, -0.002, -0.001}, + {-0.001, 0.1, 0.05, 0.05}, + {-0.1, 0.001, -0.05, -0.05}, + {-0.1, 0.001, 0.002, 0.001}, + }}; + + for (const auto& test_case : rate_cases) + { + auto data = makeResidualData(); + data.parameters[Params::dPmin] = test_case.minimum; + data.parameters[Params::dPmax] = test_case.maximum; + + Fixture fixture(data); + fixture.attachAllInputs(); + fixture.input(Ext::PREF) = rampReference(0.5, test_case.input); + success *= fixture.prepare(0.0, 0.2); + setControlState(fixture.reecb); + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.reecb, + {{Vars::PORD, test_case.expected}}, + "asymmetric active-power ramp limit"); + } + } + + { + // The voltage band gates the active-power order outside it. + const std::array gate_cases{{ + {0.0, 0.0}, + {1.0, 0.2}, + {2.0, 0.0}, + }}; + for (const auto& test_case : gate_cases) + { + Fixture fixture(makeResidualData()); + fixture.attachAllInputs(); + fixture.input(Ext::PREF) = rampReference(0.5, 0.2); + success *= fixture.prepare(0.0, 0.2); + setControlState(fixture.reecb); + setState(fixture.reecb, {{Vars::VT, test_case.input}}); + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.reecb, + {{Vars::PORD, test_case.expected}}, + "active-power voltage gate"); + } + } + + { + // Saturated probes sit beyond their limit by a margin, so a blocked + // gate contributes nothing and an admitted gate passes the full rate. + const std::array antiwindup_cases{{ + {2.0, 1.0, 0.0}, + {2.0, -1.0, -0.5}, + {-0.4, -1.0, 0.0}, + {-0.4, 1.0, 0.6}, + {0.5, 1.0, 0.6}, + }}; + for (const auto& test_case : antiwindup_cases) + { + Fixture fixture(makeResidualData()); + fixture.attachAllInputs(); + fixture.input(Ext::PREF) = rampReference(test_case.state, test_case.reference); + success *= fixture.prepare(0.0, 0.2); + setControlState(fixture.reecb); + setState(fixture.reecb, {{Vars::PORD, test_case.state}}); + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.reecb, + {{Vars::PORD, test_case.expected}}, + "active-power antiwindup"); + } + } + + { + // P priority leaves the reactive command on the residual capacity. + const std::array reactive_limit_cases{{ + {-3.0, -2.0}, + {-1.0, -1.0}, + {0.0, 0.0}, + {1.0, 1.0}, + {3.0, 2.0}, + }}; + for (const auto& test_case : reactive_limit_cases) + { + Fixture fixture(makeData()); + success *= fixture.prepare(0.0, 0.2); + setControlState(fixture.reecb); + setState(fixture.reecb, + {{Vars::ILMAX, 2.0}, + {Vars::IQCMD, 0.0}, + {Vars::QV, test_case.input}}); + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.reecb, + {{Vars::IQCMD, test_case.expected}}, + "reactive-command limit"); + } + } + + { + // Q priority leaves the active command on the residual capacity, and + // the active command is one-sided. + auto data = makeData(); + data.parameters[Params::Pqflag] = false; + + const std::array active_limit_cases{{ + {-1.0, 0.0}, + {0.6, 0.6}, + {1.0, 1.0}, + {3.0, 2.0}, + }}; + for (const auto& test_case : active_limit_cases) + { + Fixture fixture(data); + success *= fixture.prepare(0.2, 0.0); + setControlState(fixture.reecb); + setState(fixture.reecb, + {{Vars::ILMAX, 2.0}, + {Vars::IPCMD, 0.0}, + {Vars::PORD, test_case.input}}); + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.reecb, + {{Vars::IPCMD, test_case.expected}}, + "active-command limit"); + } + } + + { + // The priority selector chooses which command consumes the circle. + const std::array, 2> priority_cases{{ + {true, 0.32}, + {false, 0.56}, + }}; + for (const auto& [p_priority, expected] : priority_cases) + { + auto data = makeResidualData(); + data.parameters[Params::Pqflag] = p_priority; + Fixture fixture(data, kStateVr, kStateVi); + fixture.attachAllInputs(); + setAnswerKeyInputs(fixture); + success *= fixture.prepare(0.25, 0.35); + setAnswerKeyState(fixture.reecb); + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.reecb, + {{Vars::ILMAX, expected}}, + p_priority ? "P-priority current circle" + : "Q-priority current circle"); + } + } + + { + // Selecting the priority command before forming the circle keeps an + // overflowing inactive-command factor from producing NaN. + const RealT maximum = std::numeric_limits::max(); + const RealT limit = maximum / 1024.0; + const std::array priorities{{false, true}}; + for (const bool p_priority : priorities) + { + auto data = makeData(); + data.parameters[Params::mva] = 100.0; + data.parameters[Params::Imax] = limit; + data.parameters[Params::Pqflag] = p_priority; + + Fixture fixture(data); + success *= fixture.prepare(0.0, 0.0); + setControlState(fixture.reecb); + RealT iqcmd = limit; + RealT ipcmd = maximum; + if (p_priority) + { + iqcmd = maximum; + ipcmd = limit; + } + setState(fixture.reecb, + {{Vars::ILMAX, 0.0}, + {Vars::IQCMD, iqcmd}, + {Vars::IPCMD, ipcmd}}); + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.reecb, + {{Vars::ILMAX, 0.0}}, + "finite selected current circle"); + success *= allResidualsFinite(fixture.reecb); + } + } + + { + // The signed-square continuation keeps a negative capacity iterate + // finite, and its magnitude still bounds the low-priority command. + auto data = makeData(); + data.parameters[Params::Imax] = 1.0; + data.parameters[Params::Pqflag] = true; + + const std::array continuation_cases{{ + {-0.5, 1.0}, + {0.5, 0.5}, + {0.0, 0.75}, + }}; + for (const auto& test_case : continuation_cases) + { + Fixture fixture(data); + success *= fixture.prepare(0.0, 0.25); + setControlState(fixture.reecb); + setState(fixture.reecb, + {{Vars::ILMAX, test_case.input}, + {Vars::IPCMD, 0.25}, + {Vars::IQCMD, 0.0}, + {Vars::QV, 1.0}}); + success *= (fixture.evaluate() == 0); + success *= residualsMatch(fixture.reecb, + {{Vars::ILMAX, test_case.expected}}, + "signed capacity continuation"); + + const RealT expected_command = + test_case.input == ZERO ? 0.0 : 0.5; + success *= residualsMatch(fixture.reecb, + {{Vars::IQCMD, expected_command}}, + "capacity magnitude bound"); + success *= allResidualsFinite(fixture.reecb); + } + } + + return success.report(__func__); + } + + /// Fixed coefficients and the complete structure pin every selector + /// path at a non-unit alpha. + TestOutcome dependencyTracking() + { + TestStatus success = true; + + const std::array selector_values{{false, true}}; + for (const bool pf : selector_values) + { + for (const bool voltage : selector_values) + { + for (const bool reactive : selector_values) + { + for (const bool p_priority : selector_values) + { + auto data = makeJacobianData(); + data.parameters[Params::PfFlag] = pf; + data.parameters[Params::VFlag] = voltage; + data.parameters[Params::QFlag] = reactive; + data.parameters[Params::Pqflag] = p_priority; + + const auto dependency = dependencyTrackingJacobian(data, kNonunitAlpha, success); + + success *= jacobianStructureMatches(dependency, "dependency tracking"); + success *= derivativeMatches(dependency, Vars::VMEAS, Vars::VMEAS, -5.7, "VMEAS diagonal"); + success *= derivativeMatches(dependency, Vars::VMEAS, Vars::VT, 5.0, "VMEAS-VT"); + success *= derivativeMatches(dependency, Vars::PMEAS, Vars::PMEAS, -3.2, "PMEAS diagonal"); + success *= derivativeMatches(dependency, Vars::PMEAS, kPeColumn, 5.0, "PMEAS-PE"); + success *= derivativeMatches(dependency, Vars::XPIQ, Vars::XPIQ, -kNonunitAlpha, "XPIQ diagonal"); + success *= derivativeMatches(dependency, Vars::XPIV, Vars::XPIV, -kNonunitAlpha, "XPIV diagonal"); + success *= derivativeMatches(dependency, Vars::QV, Vars::QV, -kNonunitAlpha - (reactive ? 0.0 : 2.0), "QV diagonal"); + success *= derivativeMatches(dependency, Vars::PORD, Vars::PORD, -4.7, "PORD diagonal"); + success *= derivativeMatches(dependency, Vars::PORD, kPrefColumn, 8.0, "PORD-PREF"); + success *= derivativeMatches(dependency, Vars::VT, Vars::VT, -2.0, "VT diagonal"); + success *= derivativeMatches(dependency, Vars::VT, kBusVrColumn, 1.8, "VT-Vr"); + success *= derivativeMatches(dependency, Vars::VT, kBusViColumn, 0.8, "VT-Vi"); + success *= derivativeMatches(dependency, Vars::ILMAX, Vars::ILMAX, -4.0, "ILMAX diagonal"); + success *= derivativeMatches(dependency, Vars::IQCMD, Vars::IQCMD, -2.0, "IQCMD diagonal"); + success *= derivativeMatches(dependency, Vars::IPCMD, Vars::IPCMD, -2.0, "IPCMD diagonal"); + success *= derivativeMatches(dependency, Vars::IPCMD, Vars::PORD, 1.0, "IPCMD-PORD"); + success *= derivativeMatches(dependency, Vars::IPCMD, Vars::VMEAS, -0.5, "IPCMD-VMEAS"); + success *= derivativeMatches(dependency, Vars::IQCMD, Vars::XPIV, reactive ? 1.0 : 0.0, "IQCMD-XPIV selector path"); + success *= derivativeMatches(dependency, Vars::IQCMD, Vars::QV, reactive ? 0.0 : 1.0, "IQCMD-QV selector path"); + success *= derivativeMatches( + dependency, Vars::XPIQ, kQgenColumn, reactive && voltage ? -0.8 : 0.0, "XPIQ-QGEN selector path"); + + // The direct-voltage coefficient carries no power-base factor, + // while the cascaded path converts the reference to component base. + RealT xpiv_qext = 0.0; + if (reactive) + { + xpiv_qext = voltage ? (pf ? 0.0 : 0.6) : 0.5; + } + success *= derivativeMatches(dependency, Vars::XPIV, kQextColumn, xpiv_qext, "XPIV-QEXT selector path"); + success *= derivativeMatches(dependency, Vars::QV, kQextColumn, !reactive && !pf ? 4.0 : 0.0, "QV-QEXT selector path"); + success *= derivativeMatches(dependency, Vars::QV, kPfarefColumn, !reactive && pf ? 1.0 : 0.0, "QV-PFAREF selector path"); + + if (p_priority) + { + success *= derivativeMatches(dependency, Vars::ILMAX, Vars::IPCMD, -1.6, "P-priority current-circle column"); + success *= derivativeMatches(dependency, Vars::ILMAX, Vars::IQCMD, 0.0, "P-priority inactive current-circle column"); + } + else + { + success *= derivativeMatches(dependency, Vars::ILMAX, Vars::IQCMD, -0.8, "Q-priority current-circle column"); + success *= derivativeMatches(dependency, Vars::ILMAX, Vars::IPCMD, 0.0, "Q-priority inactive current-circle column"); + } + } + } + } + } + + // A negative capacity iterate keeps the signed-square derivative. + for (const bool p_priority : selector_values) + { + auto data = makeJacobianData(); + data.parameters[Params::Pqflag] = p_priority; + + const auto dependency = dependencyTrackingJacobian(data, kNonunitAlpha, success, -2.0); + success *= jacobianStructureMatches(dependency, "dependency tracking"); + success *= derivativeMatches(dependency, Vars::ILMAX, Vars::ILMAX, -4.0, "negative capacity continuation"); + } + + const auto zero_capacity = dependencyTrackingJacobian(makeJacobianData(), kNonunitAlpha, success, 0.0); + success *= jacobianStructureMatches(zero_capacity, "dependency tracking"); + success *= derivativeMatches(zero_capacity, Vars::ILMAX, Vars::ILMAX, -std::sqrt(ReecbT::INITIALIZATION_TOLERANCE), "zero capacity continuation"); + + // The selector sweep zeroes the injection gain, so this configuration + // exercises the injection derivative on its own. + { + auto data = makeJacobianData(); + data.parameters[Params::QFlag] = false; + data.parameters[Params::kqv] = 1.0; + data.parameters[Params::dbd1] = -0.6; + data.parameters[Params::dbd2] = 0.6; + data.parameters[Params::Iql1] = -1.2; + data.parameters[Params::Iqh1] = 1.5; + data.parameters[Params::Vref0] = 2.2; + + const auto dependency = dependencyTrackingJacobian(data, kNonunitAlpha, success); + success *= jacobianStructureMatches(dependency, "dependency tracking"); + success *= derivativeMatches(dependency, Vars::IQCMD, Vars::VMEAS, -1.0, "IQCMD-VMEAS injection path"); + success *= derivativeMatches(dependency, Vars::IQCMD, Vars::QV, 1.0, "IQCMD-QV alongside injection"); + } + + return success.report(__func__); + } + +#ifdef GRIDKIT_ENABLE_ENZYME + /// Every selector mode and continuation probe has the expected structure + /// and agrees between Enzyme and dependency tracking. + TestOutcome jacobian() + { + TestStatus success = true; + + const std::array selector_values{{false, true}}; + const std::array continuation_states{{-2.0, 0.0}}; + const std::array slew_bases{{25.0, 100.0}}; + const std::array moving_band_bases{{25.0, 200.0}}; + for (const bool pf : selector_values) + { + for (const bool voltage : selector_values) + { + for (const bool reactive : selector_values) + { + for (const bool p_priority : selector_values) + { + auto data = makeJacobianData(); + data.parameters[Params::PfFlag] = pf; + data.parameters[Params::VFlag] = voltage; + data.parameters[Params::QFlag] = reactive; + data.parameters[Params::Pqflag] = p_priority; + + success *= jacobiansMatch( + dependencyTrackingJacobian(data, kNonunitAlpha, success), + enzymeJacobian(data, kNonunitAlpha, success)); + } + } + } + } + + for (const bool p_priority : selector_values) + { + auto data = makeJacobianData(); + data.parameters[Params::Pqflag] = p_priority; + + for (const RealT ilmax : continuation_states) + { + success *= jacobiansMatch( + dependencyTrackingJacobian(data, kNonunitAlpha, success, ilmax), + enzymeJacobian(data, kNonunitAlpha, success, ilmax)); + } + } + + // Base conversion drives the active-order rate through both asymmetric + // slew limits. + for (const RealT mva : slew_bases) + { + auto data = makeJacobianData(); + data.parameters[Params::mva] = mva; + success *= jacobiansMatch( + dependencyTrackingJacobian(data, kNonunitAlpha, success), + enzymeJacobian(data, kNonunitAlpha, success)); + } + + // These points place the voltage PI state below and above its moving band. + for (const RealT mva : moving_band_bases) + { + auto data = makeJacobianData(); + data.parameters[Params::mva] = mva; + success *= jacobiansMatch( + dependencyTrackingJacobian(data, kNonunitAlpha, success, 0.5), + enzymeJacobian(data, kNonunitAlpha, success, 0.5)); + } + + auto injection_data = makeJacobianData(); + injection_data.parameters[Params::QFlag] = false; + injection_data.parameters[Params::kqv] = 1.0; + injection_data.parameters[Params::dbd1] = -0.6; + injection_data.parameters[Params::dbd2] = 0.6; + injection_data.parameters[Params::Iql1] = -1.2; + injection_data.parameters[Params::Iqh1] = 1.5; + injection_data.parameters[Params::Vref0] = 2.2; + success *= jacobiansMatch( + dependencyTrackingJacobian(injection_data, kNonunitAlpha, success), + enzymeJacobian(injection_data, kNonunitAlpha, success)); + + return success.report(__func__); + } +#endif + + private: + using Params = PhasorDynamics::Controller::ReecbParameters; + using Vars = PhasorDynamics::Controller::ReecbInternalVariables; + using Ext = PhasorDynamics::Controller::ReecbExternalVariables; + using Mon = PhasorDynamics::Controller::ReecbMonitorableVariables; + using Data = PhasorDynamics::Controller::ReecbData; + using ReecbT = PhasorDynamics::Controller::Reecb; + + static constexpr size_t index(Vars variable) + { + return static_cast(variable); + } + + static constexpr size_t index(Ext variable) + { + return static_cast(variable); + } + + struct VariableValue + { + Vars variable; + RealT value; + }; + + struct DrivenCase + { + RealT input; + RealT expected; + }; + + struct AntiWindupCase + { + RealT state; + RealT reference; + RealT expected; + }; + + /// Owns the terminal bus, REECB, the assigned command nodes, and the + /// attached input nodes. Signal storage precedes the model so every + /// referenced node outlives REECB; copying would invalidate the model + /// and node pointers. + template + class Fixture + { + private: + std::array input_values_{}; + std::array input_indices_{}; + std::array, index(Ext::MAXIMUM)> input_nodes_{}; + + PhasorDynamics::SignalNode iqcmd_node_; + PhasorDynamics::SignalNode ipcmd_node_; + bool commands_assigned_{true}; + + public: + explicit Fixture(const Data& data, + RealT vr = 1.0, + RealT vi = 0.0, + RealT system_va_base = kSystemBaseVa, + bool assign_commands = true) + : commands_assigned_(assign_commands), + bus(static_cast(vr), static_cast(vi)), + reecb(&bus, data) + { + reecb.setSystemBase(kNominalFrequency, system_va_base); + if (commands_assigned_) + { + reecb.getSignals().template assignSignalNode(&iqcmd_node_); + reecb.getSignals().template assignSignalNode(&ipcmd_node_); + } + } + + Fixture(const Fixture&) = delete; + Fixture& operator=(const Fixture&) = delete; + + void attachAllInputs(RealT initial_value = 0.0) + { + const IdxT external_index_base = reecb.size() + bus.size(); + for (size_t port = 0; port < index(Ext::MAXIMUM); ++port) + { + input_values_[port] = static_cast(initial_value); + input_indices_[port] = external_index_base + static_cast(port); + input_nodes_[port].set(&input_values_[port], &input_indices_[port]); + } + + auto& signals = reecb.getSignals(); + signals.template attachSignalNode(&input_nodes_[index(Ext::PE)]); + signals.template attachSignalNode(&input_nodes_[index(Ext::QGEN)]); + signals.template attachSignalNode(&input_nodes_[index(Ext::QEXT)]); + signals.template attachSignalNode(&input_nodes_[index(Ext::PFAREF)]); + signals.template attachSignalNode(&input_nodes_[index(Ext::PREF)]); + } + + void setCommands(RealT iqcmd, RealT ipcmd) + { + if (commands_assigned_) + { + iqcmd_node_.init(static_cast(iqcmd)); + ipcmd_node_.init(static_cast(ipcmd)); + return; + } + auto* y = reecb.y().getData(); + y[index(Vars::IQCMD)] = static_cast(iqcmd); + y[index(Vars::IPCMD)] = static_cast(ipcmd); + reecb.y().setDataUpdated(); + } + + /// Arrange the allocation, verification, bus, and command prerequisites. + bool prepare(RealT iqcmd, RealT ipcmd) + { + const bool ready = (bus.allocate() == 0) && (reecb.allocate() == 0) + && (reecb.verify() == 0) && (bus.initialize() == 0); + if (!ready) + { + std::cout << "REECB fixture preparation failed\n"; + return false; + } + setCommands(iqcmd, ipcmd); + return true; + } + + bool initialize(RealT iqcmd, RealT ipcmd) + { + if (!prepare(iqcmd, ipcmd)) + { + return false; + } + if (reecb.initialize() != 0) + { + std::cout << "REECB initialization failed\n"; + return false; + } + return true; + } + + int evaluate() + { + return reecb.evaluateResidual(); + } + + T iqcmd() const + { + return reecb.y().getData()[index(Vars::IQCMD)]; + } + + T ipcmd() const + { + return reecb.y().getData()[index(Vars::IPCMD)]; + } + + T& input(Ext port) + { + return input_values_[index(port)]; + } + + IdxT inputIndex(Ext port) const + { + return input_indices_[index(port)]; + } + + PhasorDynamics::Bus bus; + PhasorDynamics::Controller::Reecb reecb; + }; + + static constexpr RealT kSystemBaseVa = static_cast(100.0e6); + static constexpr RealT kNominalFrequency = static_cast(60.0); + static constexpr RealT kStateVr = 0.9; + static constexpr RealT kStateVi = 0.4; + // The commands, current circle, and voltage give a well-conditioned + // interior initialization point. + static constexpr RealT kInitialIqcmd = 0.75; + static constexpr RealT kInitialIpcmd = 0.75; + static constexpr RealT kNonunitAlpha = 0.7; + + inline static const RealT kUnitSlopeAngle = std::atan(ONE); + + static constexpr size_t kBusVrColumn = index(Vars::MAXIMUM); + static constexpr size_t kBusViColumn = kBusVrColumn + 1; + static constexpr size_t kExternalColumnBase = kBusViColumn + 1; + static constexpr size_t kPeColumn = kExternalColumnBase + index(Ext::PE); + static constexpr size_t kQgenColumn = kExternalColumnBase + index(Ext::QGEN); + static constexpr size_t kQextColumn = kExternalColumnBase + index(Ext::QEXT); + static constexpr size_t kPfarefColumn = kExternalColumnBase + index(Ext::PFAREF); + static constexpr size_t kPrefColumn = kExternalColumnBase + index(Ext::PREF); + + static std::array, index(Vars::MAXIMUM)> + expectedJacobianStructure() + { + return {{ + {index(Vars::VMEAS), index(Vars::VT)}, + {index(Vars::PMEAS), kPeColumn}, + {index(Vars::XPIQ), + index(Vars::VT), + index(Vars::PMEAS), + kQgenColumn, + kQextColumn, + kPfarefColumn}, + {index(Vars::XPIV), + index(Vars::VT), + index(Vars::ILMAX), + index(Vars::VMEAS), + index(Vars::PMEAS), + index(Vars::XPIQ), + kQgenColumn, + kQextColumn, + kPfarefColumn}, + {index(Vars::QV), + index(Vars::VT), + index(Vars::VMEAS), + index(Vars::PMEAS), + kQextColumn, + kPfarefColumn}, + {index(Vars::PORD), index(Vars::VT), kPrefColumn}, + {index(Vars::VT), kBusVrColumn, kBusViColumn}, + {index(Vars::ILMAX), index(Vars::IQCMD), index(Vars::IPCMD)}, + {index(Vars::IQCMD), + index(Vars::VMEAS), + index(Vars::PMEAS), + index(Vars::XPIQ), + index(Vars::XPIV), + index(Vars::QV), + index(Vars::ILMAX), + kQgenColumn, + kQextColumn, + kPfarefColumn}, + {index(Vars::IPCMD), + index(Vars::PORD), + index(Vars::VMEAS), + index(Vars::ILMAX)}, + }}; + } + + Data makeMinimalData() const + { + Data data; + data.device_class = "Reecb"; + data.disambiguation_string = "reecb_test"; + data.monitored_variables.insert(Mon::iqcmd); + data.monitored_variables.insert(Mon::ipcmd); + data.monitored_variables.insert(Mon::vmeas); + data.monitored_variables.insert(Mon::pmeas); + return data; + } + + Data makeExplicitDefaultData() const + { + auto data = makeMinimalData(); + + // These are the documented defaults; Vref0 has no fixed default and is + // resolved from the terminal voltage. + data.parameters[Params::mva] = 100.0; + data.parameters[Params::PfFlag] = false; + data.parameters[Params::VFlag] = false; + data.parameters[Params::QFlag] = false; + data.parameters[Params::Pqflag] = false; + data.parameters[Params::Trv] = 0.02; + data.parameters[Params::Tp] = 0.0; + data.parameters[Params::Vdip] = 0.85; + data.parameters[Params::Vup] = 1.15; + data.parameters[Params::dbd1] = 0.0; + data.parameters[Params::dbd2] = 0.0; + data.parameters[Params::kqv] = 5.0; + data.parameters[Params::Iql1] = -1.1; + data.parameters[Params::Iqh1] = 1.1; + data.parameters[Params::Qmax] = 0.436; + data.parameters[Params::Qmin] = -0.436; + data.parameters[Params::Kqp] = 0.0; + data.parameters[Params::Kqi] = 0.1; + data.parameters[Params::Vmax] = 1.1; + data.parameters[Params::Vmin] = 0.9; + data.parameters[Params::Kvp] = 18.0; + data.parameters[Params::Kvi] = 5.0; + data.parameters[Params::Tiq] = 0.02; + data.parameters[Params::Tpord] = 0.02; + data.parameters[Params::dPmax] = 99.0; + data.parameters[Params::dPmin] = -99.0; + data.parameters[Params::Pmax] = 1.0; + data.parameters[Params::Pmin] = 0.0; + data.parameters[Params::Imax] = 1.3; + return data; + } + + /// The routine fixture: a half-size component base, wide bands, and + /// limits far enough from the canonical commands that every smooth + /// transition is saturated. + Data makeData() const + { + auto data = makeMinimalData(); + + data.parameters[Params::mva] = 50.0; + data.parameters[Params::PfFlag] = false; + data.parameters[Params::VFlag] = true; + data.parameters[Params::QFlag] = false; + data.parameters[Params::Pqflag] = true; + data.parameters[Params::Trv] = 0.02; + data.parameters[Params::Tp] = 0.02; + data.parameters[Params::Vref0] = 1.0; + data.parameters[Params::Vdip] = 0.5; + data.parameters[Params::Vup] = 1.5; + data.parameters[Params::dbd1] = -0.2; + data.parameters[Params::dbd2] = 0.2; + data.parameters[Params::kqv] = 0.0; + data.parameters[Params::Iql1] = -1.0; + data.parameters[Params::Iqh1] = 1.0; + data.parameters[Params::Qmax] = 2.0; + data.parameters[Params::Qmin] = -2.0; + data.parameters[Params::Kqp] = 1.0; + data.parameters[Params::Kqi] = 0.0; + data.parameters[Params::Vmax] = 1.5; + data.parameters[Params::Vmin] = 0.5; + data.parameters[Params::Kvp] = 1.0; + data.parameters[Params::Kvi] = 0.0; + data.parameters[Params::Tiq] = 0.02; + data.parameters[Params::Tpord] = 0.02; + data.parameters[Params::dPmax] = 1.0; + data.parameters[Params::dPmin] = -1.0; + data.parameters[Params::Pmax] = 2.0; + data.parameters[Params::Pmin] = 0.0; + data.parameters[Params::Imax] = 2.5; + return data; + } + + /// Distinct nonzero values for every parameter. The bands are wide + /// enough, and the lag reciprocals exact enough, for probe states to + /// clear every smooth transition on an exact decimal. + Data makeResidualData() const + { + auto data = makeData(); + + data.parameters[Params::PfFlag] = false; + data.parameters[Params::VFlag] = true; + data.parameters[Params::QFlag] = true; + data.parameters[Params::Pqflag] = true; + data.parameters[Params::Trv] = 0.2; + data.parameters[Params::Tp] = 0.4; + data.parameters[Params::Vref0] = 1.5; + data.parameters[Params::kqv] = 2.0; + data.parameters[Params::Iql1] = -0.4; + data.parameters[Params::Iqh1] = 0.5; + data.parameters[Params::Qmax] = 0.8; + data.parameters[Params::Qmin] = -0.7; + data.parameters[Params::Kqp] = 0.6; + data.parameters[Params::Kqi] = 0.4; + data.parameters[Params::Vmax] = 1.6; + data.parameters[Params::Vmin] = 0.4; + data.parameters[Params::Kvp] = 1.2; + data.parameters[Params::Kvi] = 0.5; + data.parameters[Params::Tiq] = 0.5; + data.parameters[Params::Tpord] = 0.25; + data.parameters[Params::dPmax] = 0.6; + data.parameters[Params::dPmin] = -0.5; + data.parameters[Params::Pmax] = 1.4; + data.parameters[Params::Pmin] = 0.1; + data.parameters[Params::Imax] = 1.5; + return data; + } + + /// The residual parameters with the injection disabled and the reactive + /// limits widened, so the sensitivity probe sits interior everywhere. + Data makeJacobianData() const + { + auto data = makeResidualData(); + data.parameters[Params::Vref0] = 1.0; + data.parameters[Params::kqv] = 0.0; + data.parameters[Params::Qmin] = -2.0; + data.parameters[Params::Qmax] = 2.0; + data.parameters[Params::dPmin] = -0.001; + data.parameters[Params::dPmax] = 0.1; + data.parameters[Params::Imax] = 2.5; + return data; + } + + /// The active-power reference that produces a requested pre-limit ramp + /// rate at a given order, on the residual-parameter time constant. + static constexpr RealT rampReference(RealT pord, RealT raw_rate) + { + return (pord + 0.25 * raw_rate) / 2.0; + } + + template + void setAnswerKeyInputs(Fixture& fixture) const + { + fixture.input(Ext::PE) = static_cast(0.3); + fixture.input(Ext::QGEN) = static_cast(-0.1); + fixture.input(Ext::QEXT) = static_cast(0.2); + fixture.input(Ext::PFAREF) = static_cast(0.15); + fixture.input(Ext::PREF) = static_cast(0.325); + } + + /// The rich state shared by the residual answer key and the priority + /// circle. Every smooth-transition argument keeps a saturation margin, + /// so each row carries its ideal value. + template + void setAnswerKeyState(PhasorDynamics::Controller::Reecb& reecb) const + { + setState(reecb, + {{Vars::VMEAS, 0.80}, + {Vars::PMEAS, 0.55}, + {Vars::XPIQ, 0.64}, + {Vars::XPIV, -0.05}, + {Vars::QV, 0.30}, + {Vars::PORD, 0.60}, + {Vars::VT, 1.00}, + {Vars::ILMAX, 1.20}, + {Vars::IQCMD, 0.25}, + {Vars::IPCMD, 0.35}}); + setDerivative(reecb, + {{Vars::VMEAS, 0.01}, + {Vars::PMEAS, -0.02}, + {Vars::XPIQ, 0.03}, + {Vars::XPIV, -0.03}, + {Vars::QV, 0.05}, + {Vars::PORD, -0.06}}); + } + + /// A neutral driven state for the control probes: unit voltage, cleared + /// controller states, and a rested derivative. + template + void setControlState(PhasorDynamics::Controller::Reecb& reecb) const + { + reecb.yp().setToConst(static_cast(ZERO)); + setState(reecb, + {{Vars::VMEAS, 1.0}, + {Vars::PMEAS, 0.6}, + {Vars::XPIQ, 0.0}, + {Vars::XPIV, 0.0}, + {Vars::QV, 0.0}, + {Vars::PORD, 0.5}, + {Vars::VT, 1.0}, + {Vars::ILMAX, 1.4}, + {Vars::IQCMD, 0.1}, + {Vars::IPCMD, 0.2}}); + reecb.yp().setDataUpdated(); + } + + /// The sensitivity probe: one interior operating point that keeps every + /// smooth transition saturated in every selector combination. + template + void setJacobianState(Fixture& fixture, RealT ilmax) const + { + fixture.input(Ext::PE) = static_cast(0.25); + fixture.input(Ext::QGEN) = static_cast(0.5); + fixture.input(Ext::QEXT) = static_cast(0.0); + fixture.input(Ext::PFAREF) = static_cast(0.0); + fixture.input(Ext::PREF) = static_cast(0.25); + + fixture.reecb.yp().setToConst(static_cast(ZERO)); + setState(fixture.reecb, + {{Vars::VMEAS, 1.0}, + {Vars::PMEAS, 0.5}, + {Vars::XPIQ, 1.6}, + {Vars::XPIV, 0.0}, + {Vars::QV, 0.0}, + {Vars::PORD, 0.5}, + {Vars::VT, 1.0}, + {Vars::ILMAX, ilmax}, + {Vars::IQCMD, 0.1}, + {Vars::IPCMD, 0.2}}); + fixture.reecb.yp().setDataUpdated(); + } + + /// Omitting every parameter must give exactly the model built from the + /// defaults the README documents, at rest and under load. + bool defaultsMatchDocumentedValues() const + { + Fixture implicit_defaults(makeMinimalData(), kStateVr, kStateVi); + Fixture explicit_defaults(makeExplicitDefaultData(), kStateVr, kStateVi); + implicit_defaults.attachAllInputs(); + explicit_defaults.attachAllInputs(); + + bool success = implicit_defaults.initialize(0.1, 0.2) + && explicit_defaults.initialize(0.1, 0.2); + if (!success) + { + std::cout << "REECB documented-default comparison failed to initialize\n"; + return false; + } + + if (implicit_defaults.evaluate() != 0 || explicit_defaults.evaluate() != 0) + { + success = false; + } + if (!vectorsMatch(implicit_defaults.reecb.y(), + explicit_defaults.reecb.y(), + "documented-default state")) + { + success = false; + } + if (!vectorsMatch(implicit_defaults.reecb.yp(), + explicit_defaults.reecb.yp(), + "documented-default derivative")) + { + success = false; + } + if (!vectorsMatch(implicit_defaults.reecb.getResidual(), + explicit_defaults.reecb.getResidual(), + "documented-default residual")) + { + success = false; + } + for (size_t port = 0; port < index(Ext::MAXIMUM); ++port) + { + const auto variable = static_cast(port); + if (!rowMatches(implicit_defaults.input(variable), + explicit_defaults.input(variable), + "documented-default signal", + port, + "")) + { + success = false; + } + } + + setAnswerKeyInputs(implicit_defaults); + setAnswerKeyInputs(explicit_defaults); + setAnswerKeyState(implicit_defaults.reecb); + setAnswerKeyState(explicit_defaults.reecb); + if (implicit_defaults.evaluate() != 0 || explicit_defaults.evaluate() != 0) + { + success = false; + } + if (!vectorsMatch(implicit_defaults.reecb.getResidual(), + explicit_defaults.reecb.getResidual(), + "documented-default dynamic residual")) + { + success = false; + } + return success; + } + + template + bool invalidParameterCase(Params parameter, ValueT value) const + { + auto data = makeData(); + data.parameters[parameter] = value; + Fixture fixture(data); + return fixture.reecb.verify() > 0; + } + + template + bool unlinkedSignalRejected() const + { + PhasorDynamics::SignalNode unlinked_node; + Fixture fixture(makeData()); + fixture.reecb.getSignals().template attachSignalNode(&unlinked_node); + return fixture.reecb.verify() > 0; + } + + /// Fill state and derivative with a recognizable ramp, restoring the + /// aliased commands, so any write by a rejected initialization shows. + void poisonState(Fixture& fixture, RealT iqcmd, RealT ipcmd) const + { + auto* y = fixture.reecb.y().getData(); + auto* yp = fixture.reecb.yp().getData(); + for (size_t row = 0; row < index(Vars::MAXIMUM); ++row) + { + y[row] = 0.125 + 0.01 * static_cast(row); + yp[row] = -0.25 - 0.01 * static_cast(row); + } + fixture.setCommands(iqcmd, ipcmd); + fixture.reecb.y().setDataUpdated(); + fixture.reecb.yp().setDataUpdated(); + } + + bool initializationRejectedAtomically(const Data& data, + RealT iqcmd, + RealT ipcmd, + const char* label, + RealT pe = 0.6, + RealT qgen = 0.6, + RealT voltage = 1.0) const + { + Fixture fixture(data, voltage); + fixture.attachAllInputs(77.0); + fixture.input(Ext::PE) = pe; + fixture.input(Ext::QGEN) = qgen; + if (!fixture.prepare(iqcmd, ipcmd)) + { + return false; + } + + poisonState(fixture, iqcmd, ipcmd); + + const auto y_before = copyVector(fixture.reecb.y()); + const auto yp_before = copyVector(fixture.reecb.yp()); + const auto bus_before = copyVector(fixture.bus.y()); + std::array inputs_before{}; + for (size_t port = 0; port < index(Ext::MAXIMUM); ++port) + { + inputs_before[port] = fixture.input(static_cast(port)); + } + + bool success = true; + if (fixture.reecb.initialize() == 0) + { + std::cout << "Expected REECB initialization rejection: " << label << '\n'; + success = false; + } + + if (!scalarPreserved(fixture.iqcmd(), iqcmd, "rejected iqcmd preservation")) + { + success = false; + } + if (!scalarPreserved(fixture.ipcmd(), ipcmd, "rejected ipcmd preservation")) + { + success = false; + } + if (!vectorUnchanged(fixture.reecb.y(), y_before, "state")) + { + success = false; + } + if (!vectorUnchanged(fixture.reecb.yp(), yp_before, "derivative")) + { + success = false; + } + if (!vectorUnchanged(fixture.bus.y(), bus_before, "bus state")) + { + success = false; + } + for (size_t port = 0; port < index(Ext::MAXIMUM); ++port) + { + if (!valueUnchanged(fixture.input(static_cast(port)), + inputs_before[port], + "external signal", + port)) + { + success = false; + } + } + return success; + } + + template + void setState(PhasorDynamics::Controller::Reecb& reecb, + std::initializer_list values) const + { + auto* y = reecb.y().getData(); + for (const auto& [variable, value] : values) + { + y[index(variable)] = static_cast(value); + } + reecb.y().setDataUpdated(); + } + + template + void setDerivative(PhasorDynamics::Controller::Reecb& reecb, + std::initializer_list values) const + { + auto* yp = reecb.yp().getData(); + for (const auto& [variable, value] : values) + { + yp[index(variable)] = static_cast(value); + } + reecb.yp().setDataUpdated(); + } + + static const char* variableName(Vars variable) + { + static constexpr std::array names{{ + "VMEAS", + "PMEAS", + "XPIQ", + "XPIV", + "QV", + "PORD", + "VT", + "ILMAX", + "IQCMD", + "IPCMD", + }}; + return names[index(variable)]; + } + + static bool variableMatches(RealT actual, + RealT expected, + const char* what, + Vars variable, + const char* context, + RealT tolerance = kTol) + { + if (isEqual(actual, expected, tolerance)) + { + return true; + } + std::cout << "REECB " << what << ' ' << variableName(variable); + if (context[0] != '\0') + { + std::cout << ' ' << context; + } + std::cout << " mismatch: " + << std::setprecision(std::numeric_limits::max_digits10) + << actual << " != " << expected << '\n'; + return false; + } + + static bool rowMatches(RealT actual, + RealT expected, + const char* what, + size_t row, + const char* context, + RealT tolerance = kTol) + { + if (isEqual(actual, expected, tolerance)) + { + return true; + } + std::cout << "REECB " << what << " row " << row; + if (context[0] != '\0') + { + std::cout << ' ' << context; + } + std::cout << " mismatch: " + << std::setprecision(std::numeric_limits::max_digits10) + << actual << " != " << expected << '\n'; + return false; + } + + bool scalarMatches(RealT actual, + RealT expected, + const char* label, + RealT tolerance = kTol) const + { + if (isEqual(actual, expected, tolerance)) + { + return true; + } + std::cout << "REECB " << label << " mismatch: " + << std::setprecision(std::numeric_limits::max_digits10) + << actual << " != " << expected << '\n'; + return false; + } + + /// A value retains exactly what its owner supplied, including signed + /// infinities and NaN. + static bool preserved(RealT actual, RealT expected) + { + if (std::isnan(expected)) + { + return std::isnan(actual); + } + return actual == expected; + } + + bool scalarPreserved(RealT actual, RealT expected, const char* label) const + { + if (preserved(actual, expected)) + { + return true; + } + std::cout << "REECB " << label << " changed: " + << std::setprecision(std::numeric_limits::max_digits10) + << actual << " != " << expected << '\n'; + return false; + } + + static bool valueUnchanged(RealT actual, + RealT expected, + const char* what, + size_t row) + { + if (preserved(actual, expected)) + { + return true; + } + std::cout << "REECB " << what << ' ' << row << " changed: " + << std::setprecision(std::numeric_limits::max_digits10) + << actual << " != " << expected << '\n'; + return false; + } + + template + bool rowsMatch(const VectorT& vector, + const ValuesT& values, + const char* what, + const char* context) const + { + bool success = true; + const auto* vector_values = vector.getData(); + for (const auto& [variable, expected] : values) + { + if (!variableMatches(static_cast(vector_values[index(variable)]), + expected, + what, + variable, + context)) + { + success = false; + } + } + return success; + } + + bool residualsMatch(const ReecbT& reecb, + std::initializer_list values, + const char* context = "") const + { + return rowsMatch(reecb.getResidual(), values, "residual", context); + } + + template + bool residualsMatch(const ReecbT& reecb, + const std::array& values, + const char* context = "") const + { + return rowsMatch(reecb.getResidual(), values, "residual", context); + } + + bool stateMatches(const ReecbT& reecb, + std::initializer_list values, + const char* context = "") const + { + return rowsMatch(reecb.y(), values, "state", context); + } + + template + bool stateMatches(const ReecbT& reecb, + const std::array& values, + const char* context = "") const + { + return rowsMatch(reecb.y(), values, "state", context); + } + + bool allResidualsWithinInitTolerance(const ReecbT& reecb) const + { + bool success = true; + const auto* f = reecb.getResidual().getData(); + const auto* yp = reecb.yp().getData(); + for (size_t row = 0; row < index(Vars::MAXIMUM); ++row) + { + const auto variable = static_cast(row); + if (!variableMatches(f[row], + 0.0, + "residual", + variable, + "at rest", + ReecbT::INITIALIZATION_TOLERANCE)) + { + success = false; + } + if (!valueUnchanged(yp[row], 0.0, "derivative", row)) + { + success = false; + } + } + return success; + } + + bool allResidualsFinite(const ReecbT& reecb) const + { + bool success = true; + const auto* f = reecb.getResidual().getData(); + for (size_t row = 0; row < index(Vars::MAXIMUM); ++row) + { + if (!std::isfinite(f[row])) + { + std::cout << "REECB residual " << variableName(static_cast(row)) + << " is not finite\n"; + success = false; + } + } + return success; + } + + bool monitorMatches(const ReecbT& reecb, + const std::array& expected, + const char* context) const + { + RealT time = 0.0; + Model::VariableMonitorController monitor(time); + monitor.addMonitor(reecb.getMonitor()); + std::stringstream output; + monitor.addSink({Model::VariableMonitorFormat::CSV}, output); + monitor.start(); + monitor.print(); + monitor.stop(); + + std::string header; + std::string values_line; + std::getline(output, header); + std::getline(output, values_line); + + bool success = header == "t,Reecb_reecb_test_iqcmd,Reecb_reecb_test_ipcmd," + "Reecb_reecb_test_vmeas,Reecb_reecb_test_pmeas"; + + const auto values = Tokenizer(values_line, ',')(); + if (values.size() != expected.size() + 1) + { + std::cout << "REECB monitor emitted " << values.size() + << " values instead of " << expected.size() + 1 << '\n'; + return false; + } + + for (size_t i = 0; i < expected.size(); ++i) + { + if (!rowMatches(values[i + 1], expected[i], "monitor", i, context)) + { + success = false; + } + } + return success; + } + + template + std::vector copyVector(const VectorT& vector) const + { + const auto* values = vector.getData(); + std::vector snapshot(static_cast(vector.getSize())); + for (size_t row = 0; row < snapshot.size(); ++row) + { + snapshot[row] = static_cast(values[row]); + } + return snapshot; + } + + template + bool vectorUnchanged(const VectorT& vector, + const std::vector& snapshot, + const char* what) const + { + bool success = true; + const auto* values = vector.getData(); + for (size_t row = 0; row < snapshot.size(); ++row) + { + if (!valueUnchanged(static_cast(values[row]), snapshot[row], what, row)) + { + success = false; + } + } + return success; + } + + template + bool vectorsMatch(const LeftVectorT& left, + const RightVectorT& right, + const char* what) const + { + if (left.getSize() != right.getSize()) + { + std::cout << "REECB " << what << " size mismatch\n"; + return false; + } + bool success = true; + const auto* left_values = left.getData(); + const auto* right_values = right.getData(); + for (size_t row = 0; row < static_cast(left.getSize()); ++row) + { + if (!rowMatches(static_cast(left_values[row]), + static_cast(right_values[row]), + what, + row, + "")) + { + success = false; + } + } + return success; + } + + void numberVariables(Fixture& fixture, RealT alpha) const + { + auto* y = fixture.reecb.y().getData(); + auto* yp = fixture.reecb.yp().getData(); + auto* bus_y = fixture.bus.y().getData(); + + for (size_t row = 0; row < index(Vars::MAXIMUM); ++row) + { + y[row].setVariableNumber(row); + yp[row].setVariableNumber(row); + yp[row].scaleDependencies(alpha); + } + for (size_t row = 0; row < static_cast(fixture.bus.size()); ++row) + { + bus_y[row].setVariableNumber(kBusVrColumn + row); + } + for (size_t port = 0; port < index(Ext::MAXIMUM); ++port) + { + const auto variable = static_cast(port); + fixture.input(variable).setVariableNumber(fixture.inputIndex(variable)); + } + + fixture.reecb.y().setDataUpdated(); + fixture.reecb.yp().setDataUpdated(); + fixture.bus.y().setDataUpdated(); + } + + std::vector + dependencyTrackingJacobian(const Data& data, + RealT alpha, + TestStatus& success, + RealT ilmax = 2.0) const + { + using DepVar = DependencyTracking::Variable; + + Fixture fixture(data, kStateVr, kStateVi); + fixture.attachAllInputs(); + success *= fixture.prepare(0.0, 0.2); + setJacobianState(fixture, ilmax); + numberVariables(fixture, alpha); + success *= (fixture.evaluate() == 0); + + std::vector rows(index(Vars::MAXIMUM)); + const auto* f = fixture.reecb.getResidual().getData(); + for (size_t row = 0; row < rows.size(); ++row) + { + rows[row] = f[row].getDependencies(); + } + return rows; + } + + bool derivativeMatches( + const std::vector& jacobian, + Vars row, + Vars column, + RealT expected, + const char* label) const + { + return derivativeMatches(jacobian, row, index(column), expected, label); + } + + bool derivativeMatches( + const std::vector& jacobian, + Vars row, + size_t column, + RealT expected, + const char* label) const + { + const auto& dependencies = jacobian[index(row)]; + const auto entry = dependencies.find(column); + if (entry == dependencies.end()) + { + std::cout << "REECB Jacobian " << label + << " missing column " << column << '\n'; + return false; + } + + const RealT actual = entry->second; + if (isEqual(actual, expected, kTol)) + { + return true; + } + std::cout << "REECB Jacobian " << label << " mismatch: " + << std::setprecision(std::numeric_limits::max_digits10) + << actual << " != " << expected << '\n'; + return false; + } + + bool jacobianStructureMatches( + const std::vector& jacobian, + const char* source) const + { + const auto expected = expectedJacobianStructure(); + if (jacobian.size() != expected.size()) + { + std::cout << "REECB " << source + << " Jacobian row-count mismatch\n"; + return false; + } + + bool success = true; + for (size_t row = 0; row < expected.size(); ++row) + { + if (jacobian[row].size() != expected[row].size()) + { + std::cout << "REECB " << source << " Jacobian row " << row + << " column-count mismatch: " << jacobian[row].size() + << " != " << expected[row].size() << '\n'; + success = false; + } + + for (const size_t column : expected[row]) + { + if (!jacobian[row].contains(column)) + { + std::cout << "REECB " << source << " Jacobian row " << row + << " missing column " << column << '\n'; + success = false; + } + } + } + return success; + } + +#ifdef GRIDKIT_ENABLE_ENZYME + std::vector + enzymeJacobian(const Data& data, + RealT alpha, + TestStatus& success, + RealT ilmax = 2.0) const + { + Fixture fixture(data, kStateVr, kStateVi); + fixture.attachAllInputs(); + success *= fixture.prepare(0.0, 0.2); + + for (IdxT row = 0; row < fixture.bus.size(); ++row) + { + fixture.bus.setVariableIndex(row, fixture.reecb.size() + row); + } + + setJacobianState(fixture, ilmax); + fixture.reecb.updateTime(0.0, alpha); + success *= (fixture.evaluate() == 0); + success *= (fixture.reecb.evaluateJacobian() == 0); + success *= (fixture.reecb.constructCsr() == 0); + return MapFromCsr(fixture.reecb.getCsrJacobian()); + } + + bool jacobiansMatch( + const std::vector& dependency, + const std::vector& enzyme) const + { + bool success = true; + if (!jacobianStructureMatches(dependency, "dependency tracking")) + { + success = false; + } + if (!jacobianStructureMatches(enzyme, "Enzyme")) + { + success = false; + } + + if (dependency.size() != enzyme.size()) + { + std::cout << "REECB Jacobian row-count mismatch\n"; + return false; + } + + for (size_t row = 0; row < dependency.size(); ++row) + { + if (!isEqual(dependency[row], enzyme[row], kTol)) + { + std::cout << "REECB Jacobian row " << row + << " mismatch between dependency tracking and Enzyme\n"; + success = false; + } + } + return success; + } +#endif + }; + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp b/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp index 0feb2a8c7..18592dacc 100644 --- a/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp +++ b/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp @@ -1,3 +1,4 @@ +#include #include #include @@ -296,6 +297,43 @@ namespace GridKit return success.report(__func__); } + /// REECB through the production data path. + TestOutcome reecb() + { + using Data = PhasorDynamics::Controller::ReecbData; + using Buses = typename Data::Buses; + using Params = typename Data::Parameters; + using Vars = PhasorDynamics::Controller::ReecbInternalVariables; + + constexpr IdxT bus_id = static_cast(1); + + TestStatus success = true; + + PhasorDynamics::SystemModelData data; + data.bus.resize(1); + data.bus[0].bus_id = bus_id; + data.bus[0].bus_type = PhasorDynamics::BusData::BusType::SLACK; + data.bus[0].Vr0 = static_cast(1.0); + data.bus[0].Vi0 = static_cast(0.0); + + Data reecb_data; + reecb_data.buses[Buses::bus] = bus_id; + reecb_data.parameters[Params::Tp] = static_cast(0.02); + reecb_data.parameters[Params::Pmin] = static_cast(-1.0); + data.reecb.push_back(reecb_data); + + 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() == static_cast(Vars::MAXIMUM); + + return success.report(__func__); + } + TestOutcome genrou() { TestStatus success = true; diff --git a/tests/UnitTests/PhasorDynamics/SystemTests.hpp b/tests/UnitTests/PhasorDynamics/SystemTests.hpp index 968ad1462..02eba339d 100644 --- a/tests/UnitTests/PhasorDynamics/SystemTests.hpp +++ b/tests/UnitTests/PhasorDynamics/SystemTests.hpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include @@ -13,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -34,7 +36,68 @@ namespace GridKit class SystemTests { private: - using RealT = typename PhasorDynamics::Component::RealT; + using ComponentT = PhasorDynamics::Component; + using RealT = typename ComponentT::RealT; + + class InitializationFailureComponent final : public ComponentT + { + public: + InitializationFailureComponent() + { + this->size_ = static_cast(1); + } + + int setGridKitComponentID(IdxT component_id) override final + { + this->gridkit_component_id_ = component_id; + return 0; + } + + int allocate() override final + { + if (!this->allocated_) + { + this->allocateVectors(this->size_); + } + + const auto size = static_cast(this->size_); + this->tag_.assign(size, false); + this->variable_indices_.resize(size); + this->residual_indices_.resize(size); + this->allocated_ = true; + return 0; + } + + int verify() const override final + { + return 0; + } + + int initialize() override final + { + return 1; + } + + int tagDifferentiable() override final + { + return 0; + } + + int setAbsoluteTolerance(RealT) override final + { + return 0; + } + + int evaluateResidual() override final + { + return 0; + } + + int evaluateJacobian() override final + { + return this->constructCoo(); + } + }; public: SystemTests() = default; @@ -341,6 +404,35 @@ namespace GridKit return status.report(__func__); } + /// SystemModel propagates a statically valid component's initialization error. + TestOutcome componentInitializationError() + { + TestStatus success = true; + + PhasorDynamics::SystemModel system; + InitializationFailureComponent component; + system.addComponent(&component); + + success *= system.verify() == 0; + + const auto previous_verbosity = Log::verbosity(); + Log::setVerbosity(Log::Verbosity::NONE); + + if (system.hasJacobian()) + { + success *= throws([&]() + { system.allocate(); }); + } + else + { + success *= system.allocate() == 0; + success *= system.initialize() != 0; + } + + Log::setVerbosity(previous_verbosity); + return success.report(__func__); + } + #ifdef GRIDKIT_ENABLE_ENZYME TestOutcome jacobian() { diff --git a/tests/UnitTests/PhasorDynamics/runComponentConnectionTests.cpp b/tests/UnitTests/PhasorDynamics/runComponentConnectionTests.cpp index b8c7f36f5..ce2d37dbb 100644 --- a/tests/UnitTests/PhasorDynamics/runComponentConnectionTests.cpp +++ b/tests/UnitTests/PhasorDynamics/runComponentConnectionTests.cpp @@ -10,6 +10,7 @@ int main() result += test.genrouEsdc1a(); result += test.genrouHygov(); result += test.regcaRepca(); + result += test.regcaReecb(); return result.summary(); } diff --git a/tests/UnitTests/PhasorDynamics/runControllerReecbTests.cpp b/tests/UnitTests/PhasorDynamics/runControllerReecbTests.cpp new file mode 100644 index 000000000..09522fd69 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/runControllerReecbTests.cpp @@ -0,0 +1,31 @@ +#include + +#include "ControllerReecbTests.hpp" + +int main() +{ + using Log = GridKit::Utilities::Logger; + + const auto previous_verbosity = Log::verbosity(); + Log::setVerbosity(Log::Verbosity::NONE); + + GridKit::Testing::TestingResults result; + GridKit::Testing::ControllerReecbTests test; + + result += test.validation(); + result += test.initializationAndSignals(); + result += test.initializationDomain(); + result += test.initializationExactness(); + result += test.residualEquations(); + result += test.selectorConfigurations(); + result += test.voltVarReferenceBase(); + result += test.reactiveControl(); + result += test.activeCurrentControl(); + result += test.dependencyTracking(); +#ifdef GRIDKIT_ENABLE_ENZYME + result += test.jacobian(); +#endif + + Log::setVerbosity(previous_verbosity); + return result.summary(); +} diff --git a/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp b/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp index 65ec715a8..7fe1dcbe0 100644 --- a/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp +++ b/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp @@ -17,6 +17,7 @@ int main() result += test.loadZIP(); result += test.regca(); result += test.repca(); + result += test.reecb(); result += test.genrou(); result += test.genClassical(); result += test.tgov1(); diff --git a/tests/UnitTests/PhasorDynamics/runSystemTests.cpp b/tests/UnitTests/PhasorDynamics/runSystemTests.cpp index f1dd5c778..71cbe85a3 100644 --- a/tests/UnitTests/PhasorDynamics/runSystemTests.cpp +++ b/tests/UnitTests/PhasorDynamics/runSystemTests.cpp @@ -17,6 +17,7 @@ int main() #endif result += test.allocationError(); + result += test.componentInitializationError(); result += test.signalError(); return result.summary(); diff --git a/tests/UnitTests/Utilities/CaseFormatTests.hpp b/tests/UnitTests/Utilities/CaseFormatTests.hpp index e4f46bcc4..c0fef2641 100644 --- a/tests/UnitTests/Utilities/CaseFormatTests.hpp +++ b/tests/UnitTests/Utilities/CaseFormatTests.hpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include #include @@ -39,6 +40,7 @@ namespace GridKit using BusData = BusData; using BusType = typename BusData::BusType; using RegcaData = Converter::RegcaData; + using ReecbData = Controller::ReecbData; const char data[] = R"({ @@ -73,8 +75,9 @@ namespace GridKit "Tqop":0.75, "Xd":2.1, "Xdp":0.2, "Xdpp":0.18, "Xq":0.5, "Xqp": 0.0, "Xqpp":0.18, "Xl":0.15, "S10":0.0, "S12":0.0}, "mon": ["delta", "omega"] }, { "class": "Gensal", "ports": {"bus":1}, "id": "2", "params": {"p0":1.0, "q0":0.05013, "H":3.0, "D":0.0, "Ra":0.0, "Tdop":7.0, "Tdopp":0.04, "Tqopp":0.05, "Xd":2.1, "Xdp":0.2, "Xdpp":0.18, "Xq":0.5, "Xl":0.15, "S10":0.0, "S12":0.0}, "mon": ["delta", "omega"] }, - { "class": "BusFault", "ports": {"bus":1}, "id": "1", "params": {"state0": false, "R":0.0, "X":1e-3} }, - { "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"] } + { "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"] }, + { "class": "Reecb", "ports": {"bus":1}, "id": "REE1", "params": {"mva":50.0, "Pqflag":true}, "mon": ["iqcmd", "ipcmd", "vmeas", "pmeas"] }, + { "class": "BusFault", "ports": {"bus":1}, "id": "1", "params": {"state0": false, "R":0.0, "X":1e-3} } ] })"; @@ -101,6 +104,7 @@ namespace GridKit success *= result.genrou.size() == 1; success *= result.gensal.size() == 1; success *= result.regca.size() == 1; + success *= result.reecb.size() == 1; success *= result.loadz.size() == 0; success *= result.bus[0].bus_id == 1; @@ -176,6 +180,18 @@ namespace GridKit success *= result.regca[0].monitored_variables.contains(RegcaData::MonitorableVariables::ii); success *= result.regca[0].monitored_variables.contains(RegcaData::MonitorableVariables::p); success *= result.regca[0].monitored_variables.contains(RegcaData::MonitorableVariables::q); + success *= std::get(result.reecb[0].parameters[ReecbData::Parameters::mva]) == 50.0; + success *= std::get(result.reecb[0].parameters[ReecbData::Parameters::Pqflag]); + success *= result.reecb[0].buses[ReecbData::Buses::bus] == 1; + success *= result.reecb[0].disambiguation_string == "REE1"; + success *= result.reecb[0].monitored_variables.contains( + ReecbData::MonitorableVariables::iqcmd); + success *= result.reecb[0].monitored_variables.contains( + ReecbData::MonitorableVariables::ipcmd); + success *= result.reecb[0].monitored_variables.contains( + ReecbData::MonitorableVariables::vmeas); + success *= result.reecb[0].monitored_variables.contains( + ReecbData::MonitorableVariables::pmeas); success *= std::get(result.bus_fault[0].parameters[BusFaultParameters::R]) == 0.0; success *= std::get(result.bus_fault[0].parameters[BusFaultParameters::X]) == 1e-3; @@ -233,7 +249,14 @@ namespace GridKit { "signal_id": 18, "name": "Reactive Power Reference"}, { "signal_id": 19, "name": "Frequency Reference"}, { "signal_id": 20, "name": "Reactive Power Command"}, - { "signal_id": 21, "name": "Active Power Command"} + { "signal_id": 21, "name": "Active Power Command"}, + { "signal_id": 22, "name": "Electrical Power"}, + { "signal_id": 23, "name": "Reactive Power"}, + { "signal_id": 24, "name": "Reactive Reference"}, + { "signal_id": 25, "name": "Power Factor Reference"}, + { "signal_id": 26, "name": "Active Power Reference"}, + { "signal_id": 27, "name": "Reactive Current Command"}, + { "signal_id": 28, "name": "Active Current Command"} ], "devices": [ { "class": "Branch", "ports": {"bus1":1, "bus2":2}, "id": "BR1", "params": {"R":0.0, "X":0.1, "G":0.0, "B":0.0, "tap":1.05, "phase":0.1} }, @@ -247,6 +270,7 @@ namespace GridKit { "class": "Repca", "ports": {"bus":1, "ir":11, "ii":12, "p":13, "q":14, "freq":15, "vref":16, "pref":17, "qref":18, "freqref":19, "qext":20, "pext":21}, "id": "PC1", "params": {"mva":50, "VcompFlag":false, "RefFlag":true, "Freqflag":true, "Tfltr":0.2, "Vfrz":0.65, "Rc":0.02, "Xc":0.03, "Kc":0.4, "dbdlow":-0.02, "dbdupper":0.03, "emax":0.8, "emin":-0.7, "Kp":2.0, "Ki":3.0, "Qmax":0.9, "Qmin":-0.8, "Tft":0.2, "Tfv":1.5, "Tp":0.4, "fdbd1":-0.01, "fdbd2":0.015, "Ddn":2.0, "Dup":1.0, "femax":0.6, "femin":-0.5, "Kpg":1.7, "Kig":1.8, "Pmax":1.2, "Pmin":0.1, "Tlag":0.5}, "mon": ["qext", "pext", "vmeas", "qmeas", "pmeas"] }, { "class": "Ieeet1", "ports": {"bus":1, "speed": 1, "efd":3}, "id": "DV3", "params": {"Tr":0.0, "Ka":50.0, "Ta":0.04, "Ke":-0.06, "Te":0.6, "Kf":0.09, "Tf":1.46, "Vrmin":-1.0, "Vrmax":1.0, "E1":2.8, "E2":3.373, "Se1":0.04, "Se2":0.33, "Ispdlim":0.0}}, { "class": "SexsPti", "ports": {"bus":1, "efd":3}, "id": "DV4", "params": {"Ta":0.1, "Tb":0.5, "Te":0.8, "K":10.0, "Efdmax":5.0, "Efdmin":-5.0}}, + { "class": "Reecb", "ports": {"bus":1, "pe":22, "qgen":23, "qext":24, "pfaref":25, "pref":26, "iqcmd":27, "ipcmd":28}, "id": "EC1", "params": {"mva":100.0, "PfFlag":false, "VFlag":true, "QFlag":true, "Pqflag":true, "Trv":0.02, "Tp":0.05, "Vref0":1.0, "Vdip":0.85, "Vup":1.15, "dbd1":-0.01, "dbd2":0.01, "kqv":5.0, "Iql1":-1.1, "Iqh1":1.1, "Qmax":0.436, "Qmin":-0.436, "Kqp":0.1, "Kqi":0.2, "Vmax":1.1, "Vmin":0.9, "Kvp":18.0, "Kvi":5.0, "Tiq":0.02, "Tpord":0.02, "dPmax":99.0, "dPmin":-99.0, "Pmax":1.0, "Pmin":0.0, "Imax":1.3}, "mon": ["iqcmd", "ipcmd", "vmeas", "pmeas"]}, { "class": "BusFault", "ports": {"bus":1}, "id": "1", "params": {"state0": false, "R":0.0, "X":1e-3} } ] })"; @@ -274,7 +298,8 @@ namespace GridKit success *= result.loadz.size() == 0; success *= result.exciter.size() == 1; success *= result.sexspti.size() == 1; - success *= result.signal.size() == 20; + success *= result.reecb.size() == 1; + success *= result.signal.size() == 27; success *= result.bus[0].bus_id == 1; success *= result.bus[0].bus_type == BusType::DEFAULT; @@ -310,6 +335,8 @@ namespace GridKit success *= result.signal[7].name == "Governor Load Reference"; success *= result.signal[8].signal_id == 9; success *= result.signal[8].name == "Governor Auxiliary Power"; + success *= result.signal[26].signal_id == 28; + success *= result.signal[26].name == "Active Current Command"; success *= std::get(result.branch[0].parameters[BranchParameters::R]) == 0.0; success *= std::get(result.branch[0].parameters[BranchParameters::X]) == 0.1; @@ -513,6 +540,56 @@ namespace GridKit success *= result.sexspti[0].signal_outputs[Exciter::SexsPtiSignalOutputs::efd] == 3; success *= result.sexspti[0].disambiguation_string == "DV4"; + using ReecbData = Controller::ReecbData; + using ReecbParams = ReecbData::Parameters; + success *= std::get(result.reecb[0].parameters[ReecbParams::mva]) == 100.0; + success *= !std::get(result.reecb[0].parameters[ReecbParams::PfFlag]); + success *= std::get(result.reecb[0].parameters[ReecbParams::VFlag]); + success *= std::get(result.reecb[0].parameters[ReecbParams::QFlag]); + success *= std::get(result.reecb[0].parameters[ReecbParams::Pqflag]); + success *= std::get(result.reecb[0].parameters[ReecbParams::Trv]) == 0.02; + success *= std::get(result.reecb[0].parameters[ReecbParams::Tp]) == 0.05; + success *= std::get(result.reecb[0].parameters[ReecbParams::Vref0]) == 1.0; + success *= std::get(result.reecb[0].parameters[ReecbParams::Vdip]) == 0.85; + success *= std::get(result.reecb[0].parameters[ReecbParams::Vup]) == 1.15; + success *= std::get(result.reecb[0].parameters[ReecbParams::dbd1]) == -0.01; + success *= std::get(result.reecb[0].parameters[ReecbParams::dbd2]) == 0.01; + success *= std::get(result.reecb[0].parameters[ReecbParams::kqv]) == 5.0; + success *= std::get(result.reecb[0].parameters[ReecbParams::Iql1]) == -1.1; + success *= std::get(result.reecb[0].parameters[ReecbParams::Iqh1]) == 1.1; + success *= std::get(result.reecb[0].parameters[ReecbParams::Qmax]) == 0.436; + success *= std::get(result.reecb[0].parameters[ReecbParams::Qmin]) == -0.436; + success *= std::get(result.reecb[0].parameters[ReecbParams::Kqp]) == 0.1; + success *= std::get(result.reecb[0].parameters[ReecbParams::Kqi]) == 0.2; + success *= std::get(result.reecb[0].parameters[ReecbParams::Vmax]) == 1.1; + success *= std::get(result.reecb[0].parameters[ReecbParams::Vmin]) == 0.9; + success *= std::get(result.reecb[0].parameters[ReecbParams::Kvp]) == 18.0; + success *= std::get(result.reecb[0].parameters[ReecbParams::Kvi]) == 5.0; + success *= std::get(result.reecb[0].parameters[ReecbParams::Tiq]) == 0.02; + success *= std::get(result.reecb[0].parameters[ReecbParams::Tpord]) == 0.02; + success *= std::get(result.reecb[0].parameters[ReecbParams::dPmax]) == 99.0; + success *= std::get(result.reecb[0].parameters[ReecbParams::dPmin]) == -99.0; + success *= std::get(result.reecb[0].parameters[ReecbParams::Pmax]) == 1.0; + success *= std::get(result.reecb[0].parameters[ReecbParams::Pmin]) == 0.0; + success *= std::get(result.reecb[0].parameters[ReecbParams::Imax]) == 1.3; + success *= result.reecb[0].buses[ReecbData::Buses::bus] == 1; + success *= result.reecb[0].signal_inputs[ReecbData::SignalInputs::pe] == 22; + success *= result.reecb[0].signal_inputs[ReecbData::SignalInputs::qgen] == 23; + success *= result.reecb[0].signal_inputs[ReecbData::SignalInputs::qext] == 24; + success *= result.reecb[0].signal_inputs[ReecbData::SignalInputs::pfaref] == 25; + success *= result.reecb[0].signal_inputs[ReecbData::SignalInputs::pref] == 26; + success *= result.reecb[0].signal_outputs[ReecbData::SignalOutputs::iqcmd] == 27; + success *= result.reecb[0].signal_outputs[ReecbData::SignalOutputs::ipcmd] == 28; + success *= result.reecb[0].disambiguation_string == "EC1"; + success *= result.reecb[0].monitored_variables.contains( + ReecbData::MonitorableVariables::iqcmd); + success *= result.reecb[0].monitored_variables.contains( + ReecbData::MonitorableVariables::ipcmd); + success *= result.reecb[0].monitored_variables.contains( + ReecbData::MonitorableVariables::vmeas); + success *= result.reecb[0].monitored_variables.contains( + ReecbData::MonitorableVariables::pmeas); + success *= std::get(result.bus_fault[0].parameters[BusFaultParameters::R]) == 0.0; success *= std::get(result.bus_fault[0].parameters[BusFaultParameters::X]) == 1e-3; success *= !std::get(result.bus_fault[0].parameters[BusFaultParameters::state0]);