diff --git a/CHANGELOG.md b/CHANGELOG.md index ec2882196..e552d2cf3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -77,6 +77,7 @@ - Added `REGCA` converter model implementation for PhasorDynamics. - 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. ## v0.1 diff --git a/GridKit/Model/PhasorDynamics/CMakeLists.txt b/GridKit/Model/PhasorDynamics/CMakeLists.txt index fda957537..d2b9a20dc 100644 --- a/GridKit/Model/PhasorDynamics/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/CMakeLists.txt @@ -35,6 +35,7 @@ add_subdirectory(Branch) add_subdirectory(Bus) add_subdirectory(BusFault) add_subdirectory(BusToSignalAdapter) +add_subdirectory(Controller) add_subdirectory(Converter) add_subdirectory(Exciter) add_subdirectory(Governor) diff --git a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp index 51c2d78e5..21b7210ff 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 new file mode 100644 index 000000000..be3287ed4 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Controller/CMakeLists.txt @@ -0,0 +1,6 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +add_subdirectory(REPCA) diff --git a/GridKit/Model/PhasorDynamics/Controller/README.md b/GridKit/Model/PhasorDynamics/Controller/README.md new file mode 100644 index 000000000..ce5f57486 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Controller/README.md @@ -0,0 +1,10 @@ +# **Controller Models** + +## Introduction + +Controller models produce command signals for dynamic device models without +directly contributing to the network equations. + +## Types + +- Renewable Energy Plant Control Model REPCA (See [REPCA](REPCA/README.md)) diff --git a/GridKit/Model/PhasorDynamics/Controller/REPCA/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Controller/REPCA/CMakeLists.txt new file mode 100644 index 000000000..f0e2b345c --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Controller/REPCA/CMakeLists.txt @@ -0,0 +1,54 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +set(_install_headers Repca.hpp RepcaData.hpp) + +if(GRIDKIT_ENABLE_ENZYME) + gridkit_add_library( + phasor_dynamics_controller_repca + SOURCES RepcaEnzyme.cpp + HEADERS ${_install_headers} + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal + PRIVATE + ClangEnzymeFlags + COMPILE_OPTIONS + PRIVATE + -mllvm + -enzyme-auto-sparsity=1 + -fno-math-errno) +else() + gridkit_add_library( + phasor_dynamics_controller_repca + SOURCES Repca.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_repca_dependency_tracking + SOURCES RepcaDependencyTracking.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_repca) +target_link_libraries( + phasor_dynamics_components_dependency_tracking + INTERFACE GridKit::phasor_dynamics_controller_repca_dependency_tracking) diff --git a/GridKit/Model/PhasorDynamics/Controller/REPCA/README.md b/GridKit/Model/PhasorDynamics/Controller/REPCA/README.md new file mode 100644 index 000000000..0e5df068a --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Controller/REPCA/README.md @@ -0,0 +1,357 @@ +# **Renewable Energy Plant Control Model (REPCA)** + +REPCA is a WECC renewable energy plant control model for inverter-coupled +resources. + +## Notes + +- Fig. 1 hard nonlinearities use the linked CommonMath smooth approximations; + transition-point values may differ. + +> [!NOTE] +> `freq` is optional because the regulated bus does not yet expose a frequency +> signal. With frequency control and nonzero droop enabled, omitting it holds +> frequency at 1.0 p.u. and logs a warning. `freq` and `freqref` use absolute +> per-unit frequency. + +> [!WARNING] +> GridKit does not yet apply the associated generator's Governor Response Limits +> modes `Down Only` and `Fixed` to REPCA. + +## Block Diagram + +![REPCA plant-control block diagram](../../../../../docs/Figures/PhasorDynamics/REPCA/diagram.png) + +Figure 1: REPCA plant-control model. Figure courtesy of +[PowerWorld REPC_A model reference](https://www.powerworld.com/WebHelp/Content/TransientModels_PDF/Generator/Others/Plant%20Controller%20REPC_A.pdf). + +## Model Parameters + +All parameters are optional; omitted keys use the defaults below. + +Symbol | Units | JSON | Description | Default | Note +------------------------------------|-----------|-------------|---------------------------------------------------------|---------------|------ +$S^\mathrm{base}$ | [MVA] | `mva` | REPCA component power base | 100.0 | Set to the associated converter base +$s_\mathrm{comp}$ | [boolean] | `VcompFlag` | Voltage-compensation selector | `true` | `true` = line-drop compensation, `false` = reactive droop +$s_\mathrm{ref}$ | [boolean] | `RefFlag` | Reactive-loop reference selector | `true` | `true` = voltage control, `false` = reactive-power control +$s_\mathrm{freq}$ | [boolean] | `Freqflag` | Active-power output selector | `false` | `true` = command enabled, `false` = zero output +$T_\mathrm{fltr}$ | [sec] | `Tfltr` | Voltage and reactive-power filter time constant | 0.05 | +$V^\mathrm{frz}$ | [p.u.] | `Vfrz` | Reactive-power PI freeze-voltage threshold | 0.7 | +$R_c$ | [p.u.] | `Rc` | Line-drop compensation resistance | 0.0 | Component base +$X_c$ | [p.u.] | `Xc` | Line-drop compensation reactance | 0.0 | Component base +$K_c$ | [p.u.] | `Kc` | Reactive-current compensation gain | 1.0 | +$D_\mathrm{bd1}$ | [p.u.] | `dbdlow` | Lower reactive-loop deadband threshold | 0.0 | +$D_\mathrm{bd2}$ | [p.u.] | `dbdupper` | Upper reactive-loop deadband threshold | 0.0 | +$e^{\max}$ | [p.u.] | `emax` | Maximum reactive-loop error limit | 1.0 | +$e^{\min}$ | [p.u.] | `emin` | Minimum reactive-loop error limit | -1.0 | +$K_\mathrm{p}$ | [p.u.] | `Kp` | Reactive-power controller proportional gain | 10.0 | +$K_\mathrm{i}$ | [p.u./s] | `Ki` | Reactive-power controller integral gain | 10.0 | +$Q^{\max}$ | [p.u.] | `Qmax` | Maximum reactive-power command | 1.0 | Component base +$Q^{\min}$ | [p.u.] | `Qmin` | Minimum reactive-power command | -1.0 | Component base +$T_\mathrm{ft}$ | [sec] | `Tft` | Reactive-command lead time constant | 0.0 | +$T_\mathrm{fv}$ | [sec] | `Tfv` | Reactive-command lag time constant | 3.0 | +$T_\mathrm{p}$ | [sec] | `Tp` | Active-power measurement filter time constant | 0.0 | +$D_\mathrm{bd1}^{f}$ | [p.u.] | `fdbd1` | Lower frequency-error deadband threshold | 0.0 | +$D_\mathrm{bd2}^{f}$ | [p.u.] | `fdbd2` | Upper frequency-error deadband threshold | 0.0 | +$D_\mathrm{dn}$ | [p.u./p.u.] | `Ddn` | Down-regulation (overfrequency) gain | 20.0 | +$D_\mathrm{up}$ | [p.u./p.u.] | `Dup` | Up-regulation (underfrequency) gain | 0.0 | +$e_P^{\max}$ | [p.u.] | `femax` | Maximum active-power error limit | 1.0 | +$e_P^{\min}$ | [p.u.] | `femin` | Minimum active-power error limit | -1.0 | +$K_\mathrm{pg}$ | [p.u.] | `Kpg` | Active-power controller proportional gain | 10.0 | +$K_\mathrm{ig}$ | [p.u./s] | `Kig` | Active-power controller integral gain | 10.0 | +$P^{\max}$ | [p.u.] | `Pmax` | Maximum active-power command | 2.0 | Component base +$P^{\min}$ | [p.u.] | `Pmin` | Minimum active-power command | 0.0 | Component base +$T_\mathrm{lag}$ | [sec] | `Tlag` | Active-power command lag time constant | 3.0 | + +### Parameter Validation + +All real parameters must be finite. Invalid parameter sets are rejected by: + +```math +\begin{aligned} + S^\mathrm{base} &> 0 \\ + D_\mathrm{bd1} &\le 0 \le D_\mathrm{bd2} \\ + e^{\min} &\le 0 \le e^{\max} \\ + Q^{\min} &\le Q^{\max} \\ + D_\mathrm{bd1}^{f} &\le 0 \le D_\mathrm{bd2}^{f} \\ + D_\mathrm{dn},D_\mathrm{up} &\ge 0 \\ + e_P^{\min} &\le 0 \le e_P^{\max} \\ + P^{\min} &\le P^{\max} +\end{aligned} +``` + +The power bases and both conversion ratios must also be finite and positive; +`verify()` also enforces Model Ports. + +### Model Derived Parameters + +Let $\epsilon_T=10^{-3}\ \mathrm{s}$. Smaller time constants are raised to that +floor with a warning: + +```math +\begin{aligned} + T_x &\leftarrow \max(T_x,\epsilon_T), \quad x\in\{\mathrm{fltr},\mathrm{fv},\mathrm{p},\mathrm{lag}\} \\ + s_\mathrm{comp}^\mathrm{off} &= 1 - s_\mathrm{comp} \\ + s_\mathrm{ref}^\mathrm{off} &= 1 - s_\mathrm{ref} \\ + k_\mathrm{base} &= \dfrac{S^\mathrm{sys}}{S^\mathrm{base}} \\ + I_\mathrm{r}^\mathrm{comp} &= k_\mathrm{base}I_\mathrm{r}, \quad + I_\mathrm{i}^\mathrm{comp} = k_\mathrm{base}I_\mathrm{i} +\end{aligned} +``` + +$S^\mathrm{sys}$ is the system power base; $k_\mathrm{base}$ converts system-base +power and current signals to component base. GridKit does not expose `PUflag`; +these ports are always on system base. + +## Model Ports + +Name | Port | Init | Description +------------|--------|---------|------------ +`bus` | Bus | Known | Regulated-bus voltage +`ir` | Input | Known | Branch-current real component on system base +`ii` | Input | Known | Branch-current imaginary component on system base +`p` | Input | Known | Branch active power on system base +`q` | Input | Known | Branch reactive power on system base +`freq` | Input | Known | Optional absolute frequency; defaults to 1.0 p.u.[^frequency-measurement] +`vref` | Input | Unknown | Voltage-control reference +`pref` | Input | Unknown | Plant active-power reference on system base +`qref` | Input | Unknown | Reactive-power reference on system base +`freqref` | Input | Unknown | Absolute per-unit frequency reference +`qext` | Output | Known | Reactive-power command on system base +`pext` | Output | Known | Active-power command on system base + +The bus and measurement inputs are required; attached inputs must be linked. +Optional references are published when attached and latched otherwise. The +`qext` and `pext` assignments are optional; both remain monitorable. +Initialization preserves `qext` and preserves `pext` only when +$s_\mathrm{freq}=1$; otherwise `pext` becomes zero. + +## Model Variables + +### Internal Variables + +#### Differential + +Symbol | Units | Description | Note +------------------------|--------|-------------------------------------|------ +$V^\mathrm{meas}$ | [p.u.] | Filtered regulated voltage | State 1 in Fig. 1 +$Q^\mathrm{meas}$ | [p.u.] | Filtered reactive-power signal | State 2 in Fig. 1; component base +$x_Q^\mathrm{PI}$ | [p.u.] | Reactive-power PI controller state | State 3 in Fig. 1; component base +$x_Q^\mathrm{lag}$ | [p.u.] | Reactive-command lead-lag state | State 4 in Fig. 1; component base +$P^\mathrm{meas}$ | [p.u.] | Filtered active-power signal | State 5 in Fig. 1; component base +$x_P^\mathrm{PI}$ | [p.u.] | Active-power PI controller state | State 6 in Fig. 1; component base +$P^\mathrm{ref}$ | [p.u.] | Active-power command lag state | State 7 in Fig. 1; component base + +#### Algebraic + +Symbol | Units | Description | Note +--------------------------|--------|-------------------------------------|------ +$V$ | [p.u.] | Regulated-bus voltage magnitude | +$V^\mathrm{ldc}$ | [p.u.] | Line-drop compensated voltage magnitude | +$V^\mathrm{droop}$ | [p.u.] | Reactive-droop-compensated voltage | +$V^\mathrm{ctrl}$ | [p.u.] | Selected voltage-measurement input | +$s_\mathrm{frz}$ | [-] | Smooth reactive-power PI voltage-enable gate | +$e_\mathrm{RQ}$ | [p.u.] | Selected reactive-loop error | +$e_\mathrm{RQ}^\mathrm{db}$ | [p.u.] | Deadbanded reactive-loop error | +$e_\mathrm{RQ}^\mathrm{lim}$ | [p.u.] | Limited reactive-loop error | +$Q^\mathrm{PI}$ | [p.u.] | Reactive-power PI output | Component base +$Q^\mathrm{ext}$ | [p.u.] | Reactive-power command output | System base +$e_f$ | [p.u.] | Frequency error after deadband | +$e_P$ | [p.u.] | Active-power control error | Component base +$e_P^\mathrm{lim}$ | [p.u.] | Limited active-power control error | Component base +$P^\mathrm{PI}$ | [p.u.] | Active-power PI output | Component base +$P^\mathrm{ext}$ | [p.u.] | Active-power command output | System base + +### External Variables + +#### Differential + +None. + +#### Algebraic + +Symbol | Units | Init | Description | Note +-------------------------------|--------|---------|-----------------------------------|------ +$V_\mathrm{r}$ | [p.u.] | Known | Regulated-bus voltage, real component | Bus input +$V_\mathrm{i}$ | [p.u.] | Known | Regulated-bus voltage, imaginary component | Bus input +$I_\mathrm{r}$ | [p.u.] | Known | Branch-current real component | Signal port `ir`; system base +$I_\mathrm{i}$ | [p.u.] | Known | Branch-current imaginary component | Signal port `ii`; system base +$P$ | [p.u.] | Known | Branch active power | Signal port `p`; system base +$Q$ | [p.u.] | Known | Branch reactive power | Signal port `q`; system base +$f$ | [p.u.] | Known | Frequency input | Signal port `freq` +$V^\mathrm{ref}$ | [p.u.] | Unknown | Voltage-control reference | Optional signal port `vref` +$P_\mathrm{plant}^\mathrm{ref}$ | [p.u.] | Unknown | Plant active-power reference | Optional signal port `pref`; system base +$Q^\mathrm{ref}$ | [p.u.] | Unknown | Reactive-power reference | Optional signal port `qref`; system base +$f^\mathrm{ref}$ | [p.u.] | Unknown | Frequency reference | Optional signal port `freqref` + +## Model Equations + +### Internal Equations + +#### Differential + +```math +\begin{aligned} + 0 &= -\dot{V}^\mathrm{meas} + \dfrac{1}{T_\mathrm{fltr}} (V^\mathrm{ctrl} - V^\mathrm{meas}) \\ + 0 &= -\dot{Q}^\mathrm{meas} + \dfrac{1}{T_\mathrm{fltr}} (k_\mathrm{base}Q - Q^\mathrm{meas}) \\ + 0 &= -\dot{x}_Q^\mathrm{PI} + s_\mathrm{frz}\, \text{antiwindup}(Q^\mathrm{PI}, K_\mathrm{i}e_\mathrm{RQ}^\mathrm{lim};\,Q^{\min}, Q^{\max}) \\ + 0 &= -\dot{x}_Q^\mathrm{lag} + \dfrac{1}{T_\mathrm{fv}} (Q^\mathrm{PI} - x_Q^\mathrm{lag}) \\ + 0 &= -\dot{P}^\mathrm{meas} + \dfrac{1}{T_\mathrm{p}} (k_\mathrm{base}P - P^\mathrm{meas}) \\ + 0 &= -\dot{x}_P^\mathrm{PI} + \text{antiwindup}(P^\mathrm{PI}, K_\mathrm{ig}e_P^\mathrm{lim};\,P^{\min}, P^{\max}) \\ + 0 &= -\dot{P}^\mathrm{ref} + \dfrac{1}{T_\mathrm{lag}} (P^\mathrm{PI} - P^\mathrm{ref}). +\end{aligned} +``` + +CommonMath defines the [`antiwindup`](../../../../CommonMath.md#antiwindup) +target and smooth approximation. + +#### Algebraic + +```math +\begin{aligned} + 0 &= -V^2 + V_\mathrm{r}^2 + V_\mathrm{i}^2 \\ + 0 &= -(V^\mathrm{ldc})^2 + (V_\mathrm{r} - R_c I_\mathrm{r}^\mathrm{comp} + X_c I_\mathrm{i}^\mathrm{comp})^2 + (V_\mathrm{i} - R_c I_\mathrm{i}^\mathrm{comp} - X_c I_\mathrm{r}^\mathrm{comp})^2 \\ + 0 &= -V^\mathrm{droop} + V + K_c k_\mathrm{base}Q \\ + 0 &= -V^\mathrm{ctrl} + s_\mathrm{comp}V^\mathrm{ldc} + s_\mathrm{comp}^\mathrm{off}V^\mathrm{droop} \\ + 0 &= -s_\mathrm{frz} + \text{above}(V;\,V^\mathrm{frz}) \\ + 0 &= -e_\mathrm{RQ} + s_\mathrm{ref}(V^\mathrm{ref} - V^\mathrm{meas}) + s_\mathrm{ref}^\mathrm{off} (k_\mathrm{base}Q^\mathrm{ref} - Q^\mathrm{meas}) \\ + 0 &= -e_\mathrm{RQ}^\mathrm{db} + \text{deadband2}(e_\mathrm{RQ};\,D_\mathrm{bd1},D_\mathrm{bd2}) \\ + 0 &= -e_\mathrm{RQ}^\mathrm{lim} + \text{clamp}(e_\mathrm{RQ}^\mathrm{db};\,e^{\min},e^{\max}) \\ + 0 &= -Q^\mathrm{PI} + \text{clamp}(K_\mathrm{p}e_\mathrm{RQ}^\mathrm{lim}+x_Q^\mathrm{PI};\,Q^{\min},Q^{\max}) \\ + 0 &= -T_\mathrm{fv} (k_\mathrm{base}Q^\mathrm{ext}-x_Q^\mathrm{lag}) + T_\mathrm{ft} (Q^\mathrm{PI}-x_Q^\mathrm{lag}) \\ + 0 &= -e_f + \text{deadband2}(f^\mathrm{ref}-f;\,D_\mathrm{bd1}^{f},D_\mathrm{bd2}^{f}) \\ + 0 &= -e_P + k_\mathrm{base}P_\mathrm{plant}^\mathrm{ref} - P^\mathrm{meas} + \text{droop}(e_f;D_\mathrm{dn},D_\mathrm{up}) \\ + 0 &= -e_P^\mathrm{lim} + \text{clamp}(e_P;\,e_P^{\min},e_P^{\max}) \\ + 0 &= -P^\mathrm{PI} + \text{clamp}(K_\mathrm{pg}e_P^\mathrm{lim}+x_P^\mathrm{PI};\,P^{\min},P^{\max}) \\ + 0 &= -k_\mathrm{base}P^\mathrm{ext} + s_\mathrm{freq}P^\mathrm{ref}. +\end{aligned} +``` + +CommonMath defines the [derived limiter functions](../../../../CommonMath.md#derived-functions) +used above; Appendix A defines `droop`. + +### External Equations + +None. + +## Initialization + +REPCA reconstructs a steady operating point; arbitrary-state restart is unsupported. + +### Input Initialization + +```math +\begin{aligned} + V_\mathrm{r}, V_\mathrm{i} &\leftarrow \text{regulated-bus voltage} \\ + I_\mathrm{r}, I_\mathrm{i} &\leftarrow \text{branch current} \\ + P, Q &\leftarrow \text{branch power} \\ + f &\leftarrow \text{frequency input} \\ + Q^\mathrm{ext} &\leftarrow \text{known reactive-power command on system base} \\ + P^\mathrm{ext} &\leftarrow \text{known active-power command on system base}, \quad s_\mathrm{freq}=1 +\end{aligned} +``` + +### Internal Initialization + +Initialization-only $\text{clamp}^{-1}$ and $\text{deadband2}^{-1}$ recover +finite inputs that reproduce the requested smooth-block outputs. + +```math +\begin{aligned} + V &\leftarrow \sqrt{V_\mathrm{r}^2 + V_\mathrm{i}^2} \\ + V^\mathrm{ldc} &\leftarrow \sqrt{(V_\mathrm{r}-R_c I_\mathrm{r}^\mathrm{comp}+X_c I_\mathrm{i}^\mathrm{comp})^2 + (V_\mathrm{i}-R_c I_\mathrm{i}^\mathrm{comp}-X_c I_\mathrm{r}^\mathrm{comp})^2} \\ + V^\mathrm{droop} &\leftarrow V + K_c k_\mathrm{base}Q \\ + V^\mathrm{ctrl} &\leftarrow s_\mathrm{comp}V^\mathrm{ldc} + s_\mathrm{comp}^\mathrm{off}V^\mathrm{droop} \\ + V^\mathrm{meas} &\leftarrow V^\mathrm{ctrl} \\ + Q^\mathrm{meas} &\leftarrow k_\mathrm{base}Q \\ + P^\mathrm{meas} &\leftarrow k_\mathrm{base}P \\ + s_\mathrm{frz} &\leftarrow \text{above}(V;\,V^\mathrm{frz}) \\ + e_\mathrm{RQ}^\mathrm{db} &\leftarrow \text{clamp}^{-1}(0;\,e^{\min},e^{\max}) \\ + e_\mathrm{RQ} &\leftarrow \text{deadband2}^{-1}(e_\mathrm{RQ}^\mathrm{db};\,D_\mathrm{bd1},D_\mathrm{bd2}) \\ + e_\mathrm{RQ}^\mathrm{lim} &\leftarrow 0 \\ + Q^\mathrm{PI} &\leftarrow k_\mathrm{base}Q^\mathrm{ext} \\ + Q^{\min} &\leftarrow \min(Q^{\min},Q^\mathrm{PI}),\quad Q^{\max}\leftarrow \max(Q^{\max},Q^\mathrm{PI}) \\ + x_Q^\mathrm{lag} &\leftarrow Q^\mathrm{PI} \\ + u_Q^\mathrm{PI} &\leftarrow \text{clamp}^{-1}(Q^\mathrm{PI};\,Q^{\min},Q^{\max}) \\ + x_Q^\mathrm{PI} &\leftarrow u_Q^\mathrm{PI} - K_\mathrm{p}e_\mathrm{RQ}^\mathrm{lim} \\ + \Delta f_0 &\leftarrow \text{deadband2}^{-1}(0;\,D_\mathrm{bd1}^{f},D_\mathrm{bd2}^{f}) \\ + e_f &\leftarrow 0 \\ + P^\mathrm{freq} &\leftarrow \text{droop}(e_f;D_\mathrm{dn},D_\mathrm{up}) = 0 \\ + e_P &\leftarrow \text{clamp}^{-1}(0;\,e_P^{\min},e_P^{\max}) \\ + e_P^\mathrm{lim} &\leftarrow 0 \\ + P^\mathrm{ref} &\leftarrow \begin{cases} + k_\mathrm{base}P^\mathrm{ext} & s_\mathrm{freq}=1 \\ + P^\mathrm{meas} & s_\mathrm{freq}=0 + \end{cases} \\ + P^\mathrm{PI} &\leftarrow P^\mathrm{ref} \\ + P^{\min} &\leftarrow \min(P^{\min},P^\mathrm{PI}),\quad P^{\max}\leftarrow \max(P^{\max},P^\mathrm{PI}) \\ + u_P^\mathrm{PI} &\leftarrow \text{clamp}^{-1}(P^\mathrm{PI};\,P^{\min},P^{\max}) \\ + x_P^\mathrm{PI} &\leftarrow u_P^\mathrm{PI} - K_\mathrm{pg}e_P^\mathrm{lim} \\ + P^\mathrm{ext} &\leftarrow \dfrac{s_\mathrm{freq}}{k_\mathrm{base}}P^\mathrm{ref} \\ + \dot{V}^\mathrm{meas},\dot{Q}^\mathrm{meas}, \dot{x}_Q^\mathrm{PI},\dot{x}_Q^\mathrm{lag}, \dot{P}^\mathrm{meas},\dot{x}_P^\mathrm{PI},\dot{P}^\mathrm{ref} &\leftarrow 0. +\end{aligned} +``` + +If an initial PI output falls outside its configured limits, REPCA expands the +limits to include it and logs a warning. + +Initialization rejects an operating point if: + +- a required input or derived value is not finite; +- the gated reactive-power or active-power PI state rate is nonfinite or does + not vanish within the implementation tolerance. + +Initialization is atomic; candidates are validated before state or signal writes. + +### Output Initialization + +```math +\begin{aligned} + V^\mathrm{ref} &\leftarrow V^\mathrm{meas} + s_\mathrm{ref}e_\mathrm{RQ} \\ + Q^\mathrm{ref} &\leftarrow Q + \dfrac{s_\mathrm{ref}^\mathrm{off}e_\mathrm{RQ}}{k_\mathrm{base}} \\ + P_\mathrm{plant}^\mathrm{ref} &\leftarrow P + \dfrac{e_P-P^\mathrm{freq}}{k_\mathrm{base}} \\ + f^\mathrm{ref} &\leftarrow f + \Delta f_0. +\end{aligned} +``` + +## Monitorable Outputs + +Output | Units | Description | Note +----------------|--------|-------------------------------------|------ +`qext` | [p.u.] | Reactive-power command output | $Q^\mathrm{ext}$; system base +`pext` | [p.u.] | Active-power command output | $P^\mathrm{ext}$; system base +`vmeas` | [p.u.] | Filtered regulated voltage | $V^\mathrm{meas}$ +`qmeas` | [p.u.] | Filtered reactive-power signal | $Q^\mathrm{meas}$; component base +`pmeas` | [p.u.] | Filtered active-power signal | $P^\mathrm{meas}$; component base + +## Testing + +- `validation()` checks defaults, parameter domains, signal contracts, and time floors. +- `initializationAndSignals()` checks reconstruction, bases, signals, monitors, + tags, and selectors. +- `initializationDomain()` checks adjusted and collapsed limits, nonfinite + values, and atomicity. +- `residualEquations()` checks every residual against a fixed answer key. +- `reactiveControl()` checks compensation and reference modes, voltage freeze, + deadbands, smooth limits, anti-windup, and lead-lag behavior. +- `activePowerControl()` checks frequency selection, deadband, droop, smooth + limits, anti-windup, and the output lag. +- `derivatives()` checks differential-row derivative signs. +- `dependencyTracking()` checks fixed numerical and structural Jacobian oracles. +- `jacobian()` checks fixed numerical and structural oracles, plus Enzyme + agreement to $10^{-9}$ when enabled. + +## Appendix A: `droop` + +```math +\text{droop}(e;D_\mathrm{dn},D_\mathrm{up}) +=e\left[D_\mathrm{dn}+(D_\mathrm{up}-D_\mathrm{dn})\sigma(e)\right], +``` + +where $\sigma$ is GridKit's smooth +[`sigmoid`](../../../../CommonMath.md#primitives). The response preserves +$\text{droop}(0;D_\mathrm{dn},D_\mathrm{up})=0$. + +[^frequency-measurement]: Background for phase-derived, filtered frequency + measurement: [PSCAD Frequency/Phase/Magnitude Meter](https://www.pscad.com/webhelp-pscad-v5.1.0-ol/Master_Library_Models/Meters/Frequency_Phase_Magnitude_Meter.htm) and + [Ting et al., *Evaluating Methods for Measuring Grid Frequency in Low-Inertia Power Systems*](https://research-hub.nrel.gov/en/publications/evaluating-methods-for-measuring-grid-frequency-in-low-inertia-po-3/). diff --git a/GridKit/Model/PhasorDynamics/Controller/REPCA/Repca.cpp b/GridKit/Model/PhasorDynamics/Controller/REPCA/Repca.cpp new file mode 100644 index 000000000..76fc4ce09 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Controller/REPCA/Repca.cpp @@ -0,0 +1,30 @@ +/** + * @file Repca.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Non-Enzyme instantiation for the REPCA plant-control model. + */ + +#include "RepcaImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Controller + { + /** + * @brief Report that a separate Jacobian is unavailable in the plain-real build. + */ + template + int Repca::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Repca...\n"; + Log::misc() << "Jacobian evaluation is not implemented!\n"; + return 0; + } + + template class Repca; + template class Repca; + } // namespace Controller + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Controller/REPCA/Repca.hpp b/GridKit/Model/PhasorDynamics/Controller/REPCA/Repca.hpp new file mode 100644 index 000000000..73213d15d --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Controller/REPCA/Repca.hpp @@ -0,0 +1,222 @@ +/** + * @file Repca.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Declaration of the REPCA phasor-dynamics plant-control model. + */ + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + template + class BusBase; + + template + class SignalNode; + + namespace Controller + { + /// Internal variables of `Repca`. + enum class RepcaInternalVariables : size_t + { + VMEAS, ///< \f$V^\mathrm{meas}\f$ Differential filtered regulated voltage [p.u.] + QMEAS, ///< \f$Q^\mathrm{meas}\f$ Differential filtered reactive power on component base [p.u.] + XQPI, ///< \f$x_Q^\mathrm{PI}\f$ Differential reactive-power PI state on component base [p.u.] + XQLAG, ///< \f$x_Q^\mathrm{lag}\f$ Differential reactive-command lead-lag state on component base [p.u.] + PMEAS, ///< \f$P^\mathrm{meas}\f$ Differential filtered active power on component base [p.u.] + XPPI, ///< \f$x_P^\mathrm{PI}\f$ Differential active-power PI state on component base [p.u.] + PREF, ///< \f$P^\mathrm{ref}\f$ Differential active-power command lag state on component base [p.u.] + V, ///< \f$V\f$ Algebraic regulated-bus voltage magnitude [p.u.] + VLDC, ///< \f$V^\mathrm{ldc}\f$ Algebraic line-drop compensated voltage magnitude [p.u.] + VDROOP, ///< \f$V^\mathrm{droop}\f$ Algebraic reactive-droop-compensated voltage [p.u.] + VCTRL, ///< \f$V^\mathrm{ctrl}\f$ Algebraic selected voltage-measurement input [p.u.] + SFRZ, ///< \f$s_\mathrm{frz}\f$ Algebraic reactive-power PI voltage-enable gate [-] + ERQ, ///< \f$e_\mathrm{RQ}\f$ Algebraic selected reactive-loop error [p.u.] + ERQDB, ///< \f$e_\mathrm{RQ}^\mathrm{db}\f$ Algebraic deadbanded reactive-loop error [p.u.] + ERQLIM, ///< \f$e_\mathrm{RQ}^\mathrm{lim}\f$ Algebraic limited reactive-loop error [p.u.] + QPI, ///< \f$Q^\mathrm{PI}\f$ Algebraic reactive-power PI output on component base [p.u.] + QEXT, ///< \f$Q^\mathrm{ext}\f$ Algebraic reactive-power command on system base [p.u.] + EF, ///< \f$e_f\f$ Algebraic frequency error after deadband [p.u.] + EP, ///< \f$e_P\f$ Algebraic active-power control error on component base [p.u.] + EPLIM, ///< \f$e_P^\mathrm{lim}\f$ Algebraic limited active-power control error on component base [p.u.] + PPI, ///< \f$P^\mathrm{PI}\f$ Algebraic active-power PI output on component base [p.u.] + PEXT, ///< \f$P^\mathrm{ext}\f$ Algebraic active-power command on system base [p.u.] + MAXIMUM ///< Number of internal variables + }; + + /// External variables of `Repca`. + enum class RepcaExternalVariables : size_t + { + IR, ///< \f$I_\mathrm{r}\f$ Required branch-current real component on system base [p.u.] + II, ///< \f$I_\mathrm{i}\f$ Required branch-current imaginary component on system base [p.u.] + P, ///< \f$P\f$ Required branch active power on system base [p.u.] + Q, ///< \f$Q\f$ Required branch reactive power on system base [p.u.] + FREQ, ///< \f$f\f$ Optional absolute frequency input [p.u.] + VREF, ///< \f$V^\mathrm{ref}\f$ Optional voltage-control reference [p.u.] + PREF, ///< \f$P_\mathrm{plant}^\mathrm{ref}\f$ Optional plant active-power reference on system base [p.u.] + QREF, ///< \f$Q^\mathrm{ref}\f$ Optional reactive-power reference on system base [p.u.] + FREQREF, ///< \f$f^\mathrm{ref}\f$ Optional absolute frequency reference [p.u.] + MAXIMUM ///< Number of external variables + }; + + /** + * @class Repca + * @brief WECC renewable plant controller with reactive-power and + * active-power control paths. + * + * @tparam scalar_type Plain real or differentiable scalar type. + * @tparam index_type Integer index type. + */ + template + class Repca : public Component + { + using Component::abs_tol_; + using Component::allocated_; + using Component::alpha_; + using Component::f_; + using Component::gridkit_component_id_; + 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 SignalT = SignalNode; + using ModelDataT = RepcaData; + using MonitorT = Model::VariableMonitor; + using InternalVariablesT = RepcaInternalVariables; + using ExternalVariablesT = RepcaExternalVariables; + + static constexpr RealT INITIALIZATION_TOLERANCE = static_cast(1.0e-12); + + Repca(BusT* bus); + Repca(BusT* bus, const ModelDataT& data); + ~Repca(); + + 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 frequency-droop response. + static __attribute__((always_inline)) inline ScalarT droop(ScalarT error, RealT down, RealT up); + + void initializeParameters(const ModelDataT& data); + void initializeMonitor(); + void setDerivedParameters(); + + bool invertClamp(ScalarT output, RealT lower, RealT upper, ScalarT& input) const; + + bool invertDeadband(ScalarT output, RealT lower, RealT upper, ScalarT& input) const; + + static RealT logOneMinusExp(RealT x); + + [[gnu::always_inline]] inline ScalarT toComponentBase(ScalarT value) const; + ScalarT toSystemBase(ScalarT value) const; + + ScalarT& Vr(); + ScalarT& Vi(); + + static constexpr RealT TIME_CONSTANT_MINIMUM = static_cast(1.0e-3); + + static constexpr RealT INITIALIZATION_LIMIT_OFFSET = static_cast(0.1); + + BusT* bus_{nullptr}; + + RealT mva_base_{static_cast(100.0)}; + bool VcompFlag_{true}; + bool RefFlag_{true}; + bool Freqflag_{false}; + RealT Tfltr_{static_cast(0.05)}; + RealT Vfrz_{static_cast(0.7)}; + RealT Rc_{ZERO}; + RealT Xc_{ZERO}; + RealT Kc_{ONE}; + RealT dbdlow_{ZERO}; + RealT dbdupper_{ZERO}; + RealT emax_{ONE}; + RealT emin_{-ONE}; + RealT Kp_{static_cast(10.0)}; + RealT Ki_{static_cast(10.0)}; + RealT Qmax_{ONE}; + RealT Qmin_{-ONE}; + RealT Tft_{ZERO}; + RealT Tfv_{static_cast(3.0)}; + RealT Tp_{ZERO}; + RealT fdbd1_{ZERO}; + RealT fdbd2_{ZERO}; + RealT Ddn_{static_cast(20.0)}; + RealT Dup_{ZERO}; + RealT femax_{ONE}; + RealT femin_{-ONE}; + RealT Kpg_{static_cast(10.0)}; + RealT Kig_{static_cast(10.0)}; + RealT Pmax_{static_cast(2.0)}; + RealT Pmin_{ZERO}; + RealT Tlag_{static_cast(3.0)}; + + IdxT parameter_error_count_{0}; + RealT va_component_base_{ZERO}; + RealT vcomp_on_{ONE}; + RealT vcomp_off_{ZERO}; + RealT ref_on_{ONE}; + RealT ref_off_{ZERO}; + RealT freq_on_{ZERO}; + + ScalarT freqref_set_{ONE}; + ScalarT vref_set_{ONE}; + ScalarT qref_set_{ZERO}; + ScalarT pref_set_{ZERO}; + + ComponentSignals signals_; + std::unique_ptr monitor_; + + std::vector ws_; + std::vector ws_indices_; + }; + } // namespace Controller + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaData.hpp b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaData.hpp new file mode 100644 index 000000000..90e37bdf0 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaData.hpp @@ -0,0 +1,120 @@ +/** + * @file RepcaData.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Modeling data for the REPCA plant-control model. + */ + +#pragma once + +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Controller + { + /// Parameter keys for `Repca`. + enum class RepcaParameters + { + mva, ///< \f$S^\mathrm{base}\f$ Component power base [MVA] + VcompFlag, ///< \f$s_\mathrm{comp}\f$ Voltage-compensation selector [boolean] + RefFlag, ///< \f$s_\mathrm{ref}\f$ Reactive-loop reference selector [boolean] + Freqflag, ///< \f$s_\mathrm{freq}\f$ Active-power output selector [boolean] + Tfltr, ///< \f$T_\mathrm{fltr}\f$ Voltage and reactive-power filter time constant [sec] + Vfrz, ///< \f$V^\mathrm{frz}\f$ Reactive-power PI freeze threshold [p.u.] + Rc, ///< \f$R_c\f$ Line-drop resistance on component base [p.u.] + Xc, ///< \f$X_c\f$ Line-drop reactance on component base [p.u.] + Kc, ///< \f$K_c\f$ Reactive-current compensation gain [p.u.] + dbdlow, ///< \f$D_\mathrm{bd1}\f$ Lower reactive-loop deadband threshold [p.u.] + dbdupper, ///< \f$D_\mathrm{bd2}\f$ Upper reactive-loop deadband threshold [p.u.] + emax, ///< \f$e^{\max}\f$ Maximum reactive-loop error [p.u.] + emin, ///< \f$e^{\min}\f$ Minimum reactive-loop error [p.u.] + Kp, ///< \f$K_\mathrm{p}\f$ Reactive-power proportional gain [p.u.] + Ki, ///< \f$K_\mathrm{i}\f$ Reactive-power integral gain [p.u./s] + Qmax, ///< \f$Q^{\max}\f$ Maximum reactive-power command on component base [p.u.] + Qmin, ///< \f$Q^{\min}\f$ Minimum reactive-power command on component base [p.u.] + Tft, ///< \f$T_\mathrm{ft}\f$ Reactive-command lead time constant [sec] + Tfv, ///< \f$T_\mathrm{fv}\f$ Reactive-command lag time constant [sec] + Tp, ///< \f$T_\mathrm{p}\f$ Active-power measurement filter time constant [sec] + fdbd1, ///< \f$D_\mathrm{bd1}^{f}\f$ Lower frequency-error deadband threshold [p.u.] + fdbd2, ///< \f$D_\mathrm{bd2}^{f}\f$ Upper frequency-error deadband threshold [p.u.] + Ddn, ///< \f$D_\mathrm{dn}\f$ Down-regulation (overfrequency) gain [p.u./p.u.] + Dup, ///< \f$D_\mathrm{up}\f$ Up-regulation (underfrequency) gain [p.u./p.u.] + femax, ///< \f$e_P^{\max}\f$ Maximum active-power error [p.u.] + femin, ///< \f$e_P^{\min}\f$ Minimum active-power error [p.u.] + Kpg, ///< \f$K_\mathrm{pg}\f$ Active-power proportional gain [p.u.] + Kig, ///< \f$K_\mathrm{ig}\f$ Active-power integral gain [p.u./s] + Pmax, ///< \f$P^{\max}\f$ Maximum active-power command on component base [p.u.] + Pmin, ///< \f$P^{\min}\f$ Minimum active-power command on component base [p.u.] + Tlag ///< \f$T_\mathrm{lag}\f$ Active-power command lag time constant [sec] + }; + + /// Buses for `Repca`. + enum class RepcaBuses : size_t + { + bus, ///< \f$V_\mathrm{r},V_\mathrm{i}\f$ Required Known regulated-bus voltage [p.u.] + SIZE ///< Number of REPCA bus ports + }; + + /// Signal inputs for the `Repca`. + enum class RepcaSignalInputs : size_t + { + ir, ///< \f$I_\mathrm{r}\f$ Required Known branch-current real input on system base [p.u.] + ii, ///< \f$I_\mathrm{i}\f$ Required Known branch-current imaginary input on system base [p.u.] + p, ///< \f$P\f$ Required Known branch active-power input on system base [p.u.] + q, ///< \f$Q\f$ Required Known branch reactive-power input on system base [p.u.] + freq, ///< \f$f\f$ Optional Known frequency input [p.u.] + vref, ///< \f$V^\mathrm{ref}\f$ Optional Unknown voltage-reference input [p.u.] + pref, ///< \f$P_\mathrm{plant}^\mathrm{ref}\f$ Optional Unknown plant active-power reference on system base [p.u.] + qref, ///< \f$Q^\mathrm{ref}\f$ Optional Unknown reactive-power reference on system base [p.u.] + freqref, ///< \f$f^\mathrm{ref}\f$ Optional Unknown frequency-reference input [p.u.] + SIZE ///< Number of REPCA signal-input ports + }; + + /// Signal outputs for `Repca`. + enum class RepcaSignalOutputs : size_t + { + qext, ///< \f$Q^\mathrm{ext}\f$ Optional Known reactive-power command output on system base [p.u.] + pext, ///< \f$P^\mathrm{ext}\f$ Optional Known active-power command output on system base [p.u.] + SIZE ///< Number of REPCA signal-output ports + }; + + /// Variables available through the monitor interface. + enum class RepcaMonitorableVariables + { + qext, ///< \f$Q^\mathrm{ext}\f$ Reactive-power command output on system base [p.u.] + pext, ///< \f$P^\mathrm{ext}\f$ Active-power command output on system base [p.u.] + vmeas, ///< \f$V^\mathrm{meas}\f$ Filtered regulated voltage [p.u.] + qmeas, ///< \f$Q^\mathrm{meas}\f$ Filtered reactive-power signal on component base [p.u.] + pmeas ///< \f$P^\mathrm{meas}\f$ Filtered active-power signal on component base [p.u.] + }; + + /** + * @brief Model data for REPCA parameters, bus and signal ports, and monitored variables. + * + * @tparam real_type Real parameter value type. + * @tparam index_type Integer index type. + * + * @see Repca + */ + template + struct RepcaData : public ComponentData + { + RepcaData() = default; + + using Parameters = RepcaParameters; + using Buses = RepcaBuses; + using SignalInputs = RepcaSignalInputs; + using SignalOutputs = RepcaSignalOutputs; + using MonitorableVariables = RepcaMonitorableVariables; + }; + } // namespace Controller + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaDependencyTracking.cpp new file mode 100644 index 000000000..e3d8095a7 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaDependencyTracking.cpp @@ -0,0 +1,31 @@ +/** + * @file RepcaDependencyTracking.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Dependency-tracking instantiations for the REPCA plant-control model. + */ + +#include "RepcaImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Controller + { + /** + * @brief Report that DependencyTracking exposes structure through the + * residual rather than a separately assembled Jacobian. + */ + template + int Repca::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Repca...\n"; + Log::misc() << "Jacobian evaluation is not implemented!\n"; + return 0; + } + + template class Repca; + template class Repca; + } // namespace Controller + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaEnzyme.cpp b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaEnzyme.cpp new file mode 100644 index 000000000..0f6b1451e --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaEnzyme.cpp @@ -0,0 +1,122 @@ +/** + * @file RepcaEnzyme.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Enzyme sparse Jacobian for the REPCA plant-control model. + */ + +#include + +#include "RepcaImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Controller + { + /** + * @brief Assemble the sparse REPCA Jacobian with Enzyme. + * + * Differentiates the internal residual with respect to state, derivative, + * regulated-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 Repca::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Repca...\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::Repca; + 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 Repca; + template class Repca; + } // namespace Controller + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaImpl.hpp b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaImpl.hpp new file mode 100644 index 000000000..74a939ee1 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Controller/REPCA/RepcaImpl.hpp @@ -0,0 +1,1240 @@ +/** + * @file RepcaImpl.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Definition of the REPCA phasor-dynamics plant-control model. + */ + +#pragma once + +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Controller + { + /// Logger used for REPCA diagnostics. + using Log = ::GridKit::Utilities::Logger; + + /** + * @brief Construct REPCA with its documented parameter defaults + * + * The regulated bus is retained, the model is sized, and every + * parameter keeps its documented default. No monitor or signal + * connection is created. + * + * @param[in] bus Regulated bus measured by the controller. + */ + template + Repca::Repca(BusT* bus) + : bus_(bus) + { + size_ = static_cast(RepcaInternalVariables::MAXIMUM); + setDerivedParameters(); + } + + /** + * @brief Construct REPCA from model data + * + * @param[in] bus Regulated bus measured by the controller. + * @param[in] data Model parameters and monitor selections. + */ + template + Repca::Repca(BusT* bus, const ModelDataT& data) + : bus_(bus), + monitor_(std::make_unique(data)) + { + initializeParameters(data); + initializeMonitor(); + size_ = static_cast(RepcaInternalVariables::MAXIMUM); + } + + /** + * @brief Destroy the plant controller and its optional variable monitor. + */ + template + Repca::~Repca() + { + } + + /** + * @brief Set the component ID + * + * @param[in] component_id Identifier assigned by the system model. + */ + template + int Repca::setGridKitComponentID(IdxT component_id) + { + gridkit_component_id_ = component_id; + return 0; + } + + /** + * @brief Allocate model vectors and wire assigned command outputs + * + * Sizes the state, residual, bus, and signal-interface buffers, initializes + * identity index maps, and points assigned `qext` and `pext` nodes at + * the internal system-base states that REPCA publishes. Repeated + * allocation reuses the existing model vectors and signal links. + */ + template + int Repca::allocate() + { + const auto QEXT = static_cast(RepcaInternalVariables::QEXT); + const auto PEXT = static_cast(RepcaInternalVariables::PEXT); + + 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(RepcaExternalVariables::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[QEXT], + &(this->getVariableIndex(static_cast(QEXT)))); + } + + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y[PEXT], + &(this->getVariableIndex(static_cast(PEXT)))); + } + + allocated_ = true; + return 0; + } + + /** + * @brief Validate the REPCA configuration + * + * Checks parameter-loading errors, static parameter relationships, + * system/component bases and both conversion ratios, the regulated + * bus, required measurement signals, and attached optional reference + * signals. Command-output assignment is optional. + * + * @return Number of configuration errors; zero when valid. + */ + template + int Repca::verify() const + { + int ret = static_cast(parameter_error_count_); + + auto check = [&](bool condition, const char* message) + { + if (!condition) + { + Log::error() << "Repca: " << message << '\n'; + ret += 1; + } + }; + + check(bus_ != nullptr, "regulated bus is required"); + + const bool valid_component_base = std::isfinite(mva_base_) + && mva_base_ > ZERO + && std::isfinite(va_component_base_) + && va_component_base_ > ZERO; + const bool valid_system_base = std::isfinite(va_system_base_) + && va_system_base_ > ZERO; + check(valid_component_base, + "mva must define a finite positive component power base"); + 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_ / va_component_base_; + const RealT component_to_system = va_component_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(dbdlow_ <= ZERO && ZERO <= dbdupper_, + "dbdlow <= 0 <= dbdupper is required"); + check(emin_ <= ZERO && ZERO <= emax_, + "emin <= 0 <= emax is required"); + check(Qmin_ <= Qmax_, "Qmin must be less than or equal to Qmax"); + check(fdbd1_ <= ZERO && ZERO <= fdbd2_, + "fdbd1 <= 0 <= fdbd2 is required"); + check(Ddn_ >= ZERO, "Ddn must be non-negative"); + check(Dup_ >= ZERO, "Dup must be non-negative"); + check(femin_ <= ZERO && ZERO <= femax_, + "femin <= 0 <= femax is required"); + check(Pmin_ <= Pmax_, "Pmin must be less than or equal to Pmax"); + + auto check_required_signal = [&](const char* name) + { + if (!signals_.template isAttached()) + { + Log::error() << "Repca: " << name << " signal is required\n"; + ret += 1; + } + else if (!signals_.template isLinked()) + { + Log::error() << "Repca: " << name << " signal attached with no linked source\n"; + ret += 1; + } + }; + + check_required_signal.template operator()("ir"); + check_required_signal.template operator()("ii"); + check_required_signal.template operator()("p"); + check_required_signal.template operator()("q"); + + auto check_optional_signal = [&](const char* name) + { + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "Repca: " << name << " signal attached with no linked source\n"; + ret += 1; + } + }; + + check_optional_signal.template operator()("vref"); + check_optional_signal.template operator()("pref"); + check_optional_signal.template operator()("qref"); + check_optional_signal.template operator()("freq"); + check_optional_signal.template operator()("freqref"); + + return ret; + } + + /** + * @brief Initialize REPCA from the initial plant commands + * + * Reads the required bus and measurement ports, preserves the initial + * system-base `qext` value and, when frequency control is enabled, `pext`, + * reconstructs the steady controller state, and publishes the resolved + * optional references. + * + * @pre allocate() has completed. + * @pre verify() reports a valid parameter and port configuration. + * @pre `qext` and, when frequency control is enabled, `pext` contain the + * initial plant commands. + * + * @post On failure, states, derivatives, limits, latches, and signal + * values are unchanged. + * + * @return 0 on success; nonzero when allocation, configuration, initial- + * value, candidate, or steady-state checks fail. + */ + template + int Repca::initialize() + { + const auto VMEAS = static_cast(RepcaInternalVariables::VMEAS); + const auto QMEAS = static_cast(RepcaInternalVariables::QMEAS); + const auto XQPI = static_cast(RepcaInternalVariables::XQPI); + const auto XQLAG = static_cast(RepcaInternalVariables::XQLAG); + const auto PMEAS = static_cast(RepcaInternalVariables::PMEAS); + const auto XPPI = static_cast(RepcaInternalVariables::XPPI); + const auto PREF = static_cast(RepcaInternalVariables::PREF); + const auto V = static_cast(RepcaInternalVariables::V); + const auto VLDC = static_cast(RepcaInternalVariables::VLDC); + const auto VDROOP = static_cast(RepcaInternalVariables::VDROOP); + const auto VCTRL = static_cast(RepcaInternalVariables::VCTRL); + const auto SFRZ = static_cast(RepcaInternalVariables::SFRZ); + const auto ERQ = static_cast(RepcaInternalVariables::ERQ); + const auto ERQDB = static_cast(RepcaInternalVariables::ERQDB); + const auto ERQLIM = static_cast(RepcaInternalVariables::ERQLIM); + const auto QPI = static_cast(RepcaInternalVariables::QPI); + const auto QEXT = static_cast(RepcaInternalVariables::QEXT); + const auto EF = static_cast(RepcaInternalVariables::EF); + const auto EP = static_cast(RepcaInternalVariables::EP); + const auto EPLIM = static_cast(RepcaInternalVariables::EPLIM); + const auto PPI = static_cast(RepcaInternalVariables::PPI); + const auto PEXT = static_cast(RepcaInternalVariables::PEXT); + + if (!allocated_) + { + Log::error() << "Repca: allocate must complete before initialize\n"; + return 1; + } + + if (verify() > 0) + { + Log::error() << "Repca: cannot initialize with invalid configuration\n"; + return 1; + } + + auto* y = y_.getData(); + + const ScalarT qext0_system = y[QEXT]; + const ScalarT pext0_system = y[PEXT]; + const ScalarT qext0 = toComponentBase(qext0_system); + const ScalarT pext0 = toComponentBase(pext0_system); + + const ScalarT vr = Vr(); + const ScalarT vi = Vi(); + const ScalarT ir_system = + signals_.template readExternalVariable(); + const ScalarT ii_system = + signals_.template readExternalVariable(); + const ScalarT p_system = + signals_.template readExternalVariable(); + const ScalarT q_system = + signals_.template readExternalVariable(); + const ScalarT ir = toComponentBase(ir_system); + const ScalarT ii = toComponentBase(ii_system); + ScalarT freq = static_cast(ONE); + if (signals_.template isAttached()) + { + freq = signals_.template readExternalVariable(); + } + + auto is_finite = [](ScalarT value) + { + return std::isfinite(static_cast(value)); + }; + if (!is_finite(vr) || !is_finite(vi) || !is_finite(ir_system) + || !is_finite(ii_system) || !is_finite(p_system) + || !is_finite(q_system) || !is_finite(freq) || !is_finite(qext0) + || (Freqflag_ && !is_finite(pext0))) + { + Log::error() << "Repca: initial bus, signal, and command values must be finite\n"; + return 1; + } + + const ScalarT p = toComponentBase(p_system); + const ScalarT q = toComponentBase(q_system); + + const ScalarT vldc_r = vr - Rc_ * ir + Xc_ * ii; + const ScalarT vldc_i = vi - Rc_ * ii - Xc_ * ir; + + const ScalarT v0 = std::sqrt(vr * vr + vi * vi); + const ScalarT vldc0 = std::sqrt(vldc_r * vldc_r + vldc_i * vldc_i); + const ScalarT vdroop0 = v0 + Kc_ * q; + const ScalarT vctrl0 = vcomp_on_ * vldc0 + vcomp_off_ * vdroop0; + const ScalarT vmeas0 = vctrl0; + const ScalarT qmeas0 = q; + const ScalarT pmeas0 = p; + const ScalarT sfrz0 = Math::above(v0, Vfrz_); + + const ScalarT zero = static_cast(ZERO); + ScalarT erq0{}; + ScalarT erqdb0{}; + if (!invertClamp(zero, emin_, emax_, erqdb0) + || !invertDeadband(erqdb0, dbdlow_, dbdupper_, erq0)) + { + Log::error() << "Repca: reactive error blocks have no finite steady input\n"; + return 1; + } + const ScalarT erqlim0 = zero; + const ScalarT qpi0 = qext0; + const ScalarT xqlag0 = qpi0; + + const RealT qmin = std::min(Qmin_, static_cast(qpi0)); + const RealT qmax = std::max(Qmax_, static_cast(qpi0)); + + ScalarT qpi_input0{}; + if (!invertClamp(qpi0, qmin, qmax, qpi_input0)) + { + Log::error() << "Repca: reactive-power PI limiter has no finite steady input\n"; + return 1; + } + const ScalarT xqpi0 = qpi_input0 - Kp_ * erqlim0; + const ScalarT q_aw_rate0 = Math::antiwindup(qpi0, Ki_ * erqlim0, qmin, qmax); + const ScalarT xqpi_rate0 = sfrz0 * q_aw_rate0; + if (!is_finite(q_aw_rate0) || !is_finite(xqpi_rate0) + || std::abs(static_cast(xqpi_rate0)) > INITIALIZATION_TOLERANCE) + { + Log::error() << "Repca: reactive-power PI state rate is nonzero at initialization\n"; + return 1; + } + + ScalarT frequency_error0{}; + ScalarT ep0{}; + if (!invertDeadband(zero, fdbd1_, fdbd2_, frequency_error0) + || !invertClamp(zero, femin_, femax_, ep0)) + { + Log::error() << "Repca: active error blocks have no finite steady input\n"; + return 1; + } + const ScalarT ef0 = zero; + const ScalarT pfreq0 = droop(ef0, Ddn_, Dup_); + const ScalarT eplim0 = zero; + const ScalarT pref0 = Freqflag_ ? pext0 : pmeas0; + const ScalarT ppi0 = pref0; + + const RealT pmin = std::min(Pmin_, static_cast(ppi0)); + const RealT pmax = std::max(Pmax_, static_cast(ppi0)); + + ScalarT ppi_input0{}; + if (!invertClamp(ppi0, pmin, pmax, ppi_input0)) + { + Log::error() << "Repca: active-power PI limiter has no finite steady input\n"; + return 1; + } + const ScalarT xppi0 = ppi_input0 - Kpg_ * eplim0; + const ScalarT p_aw_rate0 = Math::antiwindup(ppi0, Kig_ * eplim0, pmin, pmax); + if (!is_finite(p_aw_rate0) + || std::abs(static_cast(p_aw_rate0)) > INITIALIZATION_TOLERANCE) + { + Log::error() << "Repca: active-power PI antiwindup rate is nonzero at initialization\n"; + return 1; + } + + const ScalarT pext_output0 = Freqflag_ ? pext0_system : zero; + const ScalarT freqref0 = freq + frequency_error0; + ScalarT vref0 = vmeas0; + ScalarT qref0_system = q_system; + if (RefFlag_) + { + vref0 += erq0; + } + else + { + qref0_system = toSystemBase(qmeas0 + erq0); + } + const ScalarT pref0_system = p_system + toSystemBase(ep0 - pfreq0); + + const bool candidates_are_finite = + is_finite(qext0_system) + && (!Freqflag_ || is_finite(pext0_system)) + && is_finite(qext0) + && (!Freqflag_ || is_finite(pext0)) + && is_finite(p) + && is_finite(q) + && is_finite(vldc_r) + && is_finite(vldc_i) + && is_finite(v0) + && is_finite(vldc0) + && is_finite(vdroop0) + && is_finite(vctrl0) + && is_finite(vmeas0) + && is_finite(qmeas0) + && is_finite(pmeas0) + && is_finite(sfrz0) + && is_finite(erq0) + && is_finite(erqdb0) + && is_finite(erqlim0) + && is_finite(qpi0) + && is_finite(xqlag0) + && is_finite(qpi_input0) + && is_finite(xqpi0) + && is_finite(q_aw_rate0) + && is_finite(xqpi_rate0) + && is_finite(ef0) + && is_finite(pfreq0) + && is_finite(ep0) + && is_finite(eplim0) + && is_finite(pref0) + && is_finite(ppi0) + && is_finite(ppi_input0) + && is_finite(xppi0) + && is_finite(p_aw_rate0) + && is_finite(pext_output0) + && is_finite(freqref0) + && is_finite(vref0) + && is_finite(qref0_system) + && is_finite(pref0_system); + if (!candidates_are_finite) + { + Log::error() << "Repca: derived initial values must be finite\n"; + return 1; + } + + y[VMEAS] = vmeas0; + y[QMEAS] = qmeas0; + y[XQPI] = xqpi0; + y[XQLAG] = xqlag0; + y[PMEAS] = pmeas0; + y[XPPI] = xppi0; + y[PREF] = pref0; + y[V] = v0; + y[VLDC] = vldc0; + y[VDROOP] = vdroop0; + y[VCTRL] = vctrl0; + y[SFRZ] = sfrz0; + y[ERQ] = erq0; + y[ERQDB] = erqdb0; + y[ERQLIM] = erqlim0; + y[QPI] = qpi0; + y[QEXT] = qext0_system; + y[EF] = ef0; + y[EP] = ep0; + y[EPLIM] = eplim0; + y[PPI] = ppi0; + y[PEXT] = pext_output0; + + const bool q_adjusted = qmin != Qmin_ || qmax != Qmax_; + const bool p_adjusted = pmin != Pmin_ || pmax != Pmax_; + Qmin_ = qmin; + Qmax_ = qmax; + Pmin_ = pmin; + Pmax_ = pmax; + + freqref_set_ = freqref0; + vref_set_ = vref0; + qref_set_ = qref0_system; + pref_set_ = pref0_system; + + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(vref_set_); + } + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable( + pref_set_); + } + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(qref_set_); + } + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable( + freqref_set_); + } + + if (q_adjusted) + { + Log::warning() << "Repca: initial reactive PI output is outside [Qmin, Qmax]; " + "the limits are adjusted to include the initialized value\n"; + } + if (p_adjusted) + { + Log::warning() << "Repca: initial active PI output is outside [Pmin, Pmax]; " + "the limits are adjusted to include the initialized value\n"; + } + if (Freqflag_ + && (Ddn_ != ZERO || Dup_ != ZERO) + && !signals_.template isAttached()) + { + Log::warning() << "Repca: Freqflag is enabled without a freq signal; " + "frequency remains nominal and cannot track bus-frequency deviations\n"; + } + + y_.setDataUpdated(); + yp_.setToConst(static_cast(ZERO)); + return 0; + } + + /** + * @brief Identify the differential variables + * + * The voltage and power filters, reactive PI and lead-lag states, and + * active PI and command-lag states carry derivatives; every other + * internal variable is algebraic. + */ + template + int Repca::tagDifferentiable() + { + const auto VMEAS = static_cast(RepcaInternalVariables::VMEAS); + const auto QMEAS = static_cast(RepcaInternalVariables::QMEAS); + const auto XQPI = static_cast(RepcaInternalVariables::XQPI); + const auto XQLAG = static_cast(RepcaInternalVariables::XQLAG); + const auto PMEAS = static_cast(RepcaInternalVariables::PMEAS); + const auto XPPI = static_cast(RepcaInternalVariables::XPPI); + const auto PREF = static_cast(RepcaInternalVariables::PREF); + + std::fill(tag_.begin(), tag_.end(), false); + tag_[VMEAS] = true; + tag_[QMEAS] = true; + tag_[XQPI] = true; + tag_[XQLAG] = true; + tag_[PMEAS] = true; + tag_[XPPI] = true; + tag_[PREF] = true; + return 0; + } + + /** + * @brief Compute the absolute tolerance for each variable in the model + * + * REPCA variables are dimensionless per-unit signals and controller + * states of the same order, so they share the relative tolerance as + * their absolute floor. + * + * @param[in] rel_tol Solver relative tolerance. + */ + template + int Repca::setAbsoluteTolerance(RealT rel_tol) + { + abs_tol_.setToConst(static_cast(rel_tol)); + return 0; + } + + /** + * @brief Evaluate the REPCA-owned residual rows + * + * Refreshes required bus and measurement inputs, starts optional + * references from the values latched by initialize(), then overwrites + * them from attached signals. REPCA contributes no bus residual. + */ + template + int Repca::evaluateResidual() + { + const auto IR = static_cast(RepcaExternalVariables::IR); + const auto II = static_cast(RepcaExternalVariables::II); + const auto P = static_cast(RepcaExternalVariables::P); + const auto Q = static_cast(RepcaExternalVariables::Q); + const auto FREQ = static_cast(RepcaExternalVariables::FREQ); + const auto VREF = static_cast(RepcaExternalVariables::VREF); + const auto PREF = static_cast(RepcaExternalVariables::PREF); + const auto QREF = static_cast(RepcaExternalVariables::QREF); + const auto FREQREF = static_cast(RepcaExternalVariables::FREQREF); + + ws_[VREF] = vref_set_; + ws_[PREF] = pref_set_; + ws_[QREF] = qref_set_; + ws_[FREQ] = static_cast(ONE); + ws_[FREQREF] = freqref_set_; + std::fill(ws_indices_.begin(), ws_indices_.end(), INVALID_INDEX); + + ws_[IR] = + signals_.template readExternalVariable(); + ws_indices_[IR] = + signals_.template readExternalVariableIndex(); + ws_[II] = + signals_.template readExternalVariable(); + ws_indices_[II] = + signals_.template readExternalVariableIndex(); + ws_[P] = + signals_.template readExternalVariable(); + ws_indices_[P] = + signals_.template readExternalVariableIndex(); + ws_[Q] = + signals_.template readExternalVariable(); + ws_indices_[Q] = + signals_.template readExternalVariableIndex(); + if (signals_.template isAttached()) + { + ws_[FREQ] = + signals_.template readExternalVariable(); + ws_indices_[FREQ] = + signals_.template readExternalVariableIndex(); + } + + if (signals_.template isAttached()) + { + ws_[VREF] = + signals_.template readExternalVariable(); + ws_indices_[VREF] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[PREF] = + signals_.template readExternalVariable(); + ws_indices_[PREF] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[QREF] = + signals_.template readExternalVariable(); + ws_indices_[QREF] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[FREQREF] = + signals_.template readExternalVariable(); + ws_indices_[FREQREF] = + signals_.template readExternalVariableIndex(); + } + + wb_[0] = Vr(); + wb_[1] = Vi(); + + const auto* y = y_.getData(); + const auto* yp = yp_.getData(); + auto* f = f_.getData(); + + evaluateInternalResidual(y, yp, wb_.data(), ws_.data(), f); + f_.setDataUpdated(); + return 0; + } + + /** + * @brief Access the REPCA signal interface + * + * @return Interface used to assign optional plant-command outputs and + * attach required measurements and optional references. + */ + template + auto Repca::getSignals() + -> ComponentSignals& + { + return signals_; + } + + /** + * @brief Access the configured monitor + * + * @return Monitor for this model, or nullptr when constructed without data. + */ + template + const Model::VariableMonitorBase* Repca::getMonitor() const + { + return monitor_.get(); + } + + /** + * @brief Evaluate the REPCA internal residual + * + * Evaluates seven differential and fifteen algebraic rows in enum order. + * Precomputed mode masks keep the differentiated path branch- and + * loop-free with a fixed dependency structure. + * + * @param[in] y Internal variables in `RepcaInternalVariables` order and + * on the bases documented by their enums. + * @param[in] yp Internal derivatives in the same enum order and bases. + * @param[in] wb Regulated-bus real and imaginary voltage components. + * @param[in] ws External signals in `RepcaExternalVariables` order. + * @param[out] f Caller-provided residual output buffer in + * `RepcaInternalVariables` order. + */ + template + [[gnu::always_inline]] inline int + Repca::evaluateInternalResidual( + const ScalarT* y, + const ScalarT* yp, + const ScalarT* wb, + const ScalarT* ws, + ScalarT* f) + { + const auto VMEAS = static_cast(RepcaInternalVariables::VMEAS); + const auto QMEAS = static_cast(RepcaInternalVariables::QMEAS); + const auto XQPI = static_cast(RepcaInternalVariables::XQPI); + const auto XQLAG = static_cast(RepcaInternalVariables::XQLAG); + const auto PMEAS = static_cast(RepcaInternalVariables::PMEAS); + const auto XPPI = static_cast(RepcaInternalVariables::XPPI); + const auto PREF_STATE = static_cast(RepcaInternalVariables::PREF); + const auto V = static_cast(RepcaInternalVariables::V); + const auto VLDC = static_cast(RepcaInternalVariables::VLDC); + const auto VDROOP = static_cast(RepcaInternalVariables::VDROOP); + const auto VCTRL = static_cast(RepcaInternalVariables::VCTRL); + const auto SFRZ = static_cast(RepcaInternalVariables::SFRZ); + const auto ERQ = static_cast(RepcaInternalVariables::ERQ); + const auto ERQDB = static_cast(RepcaInternalVariables::ERQDB); + const auto ERQLIM = static_cast(RepcaInternalVariables::ERQLIM); + const auto QPI = static_cast(RepcaInternalVariables::QPI); + const auto QEXT = static_cast(RepcaInternalVariables::QEXT); + const auto EF = static_cast(RepcaInternalVariables::EF); + const auto EP = static_cast(RepcaInternalVariables::EP); + const auto EPLIM = static_cast(RepcaInternalVariables::EPLIM); + const auto PPI = static_cast(RepcaInternalVariables::PPI); + const auto PEXT = static_cast(RepcaInternalVariables::PEXT); + + const auto IR = static_cast(RepcaExternalVariables::IR); + const auto II = static_cast(RepcaExternalVariables::II); + const auto P = static_cast(RepcaExternalVariables::P); + const auto Q = static_cast(RepcaExternalVariables::Q); + const auto FREQ = static_cast(RepcaExternalVariables::FREQ); + const auto VREF = static_cast(RepcaExternalVariables::VREF); + const auto PREF_INPUT = static_cast(RepcaExternalVariables::PREF); + const auto QREF = static_cast(RepcaExternalVariables::QREF); + const auto FREQREF = static_cast(RepcaExternalVariables::FREQREF); + + const ScalarT vmeas = y[VMEAS]; + const ScalarT qmeas = y[QMEAS]; + const ScalarT xqpi = y[XQPI]; + const ScalarT xqlag = y[XQLAG]; + const ScalarT pmeas = y[PMEAS]; + const ScalarT xppi = y[XPPI]; + const ScalarT pref = y[PREF_STATE]; + const ScalarT v = y[V]; + const ScalarT vldc = y[VLDC]; + const ScalarT vdroop = y[VDROOP]; + const ScalarT vctrl = y[VCTRL]; + const ScalarT sfrz = y[SFRZ]; + const ScalarT erq = y[ERQ]; + const ScalarT erqdb = y[ERQDB]; + const ScalarT erqlim = y[ERQLIM]; + const ScalarT qpi = y[QPI]; + const ScalarT qext = toComponentBase(y[QEXT]); + const ScalarT ef = y[EF]; + const ScalarT ep = y[EP]; + const ScalarT eplim = y[EPLIM]; + const ScalarT ppi = y[PPI]; + const ScalarT pext = toComponentBase(y[PEXT]); + + const ScalarT vmeas_dot = yp[VMEAS]; + const ScalarT qmeas_dot = yp[QMEAS]; + const ScalarT xqpi_dot = yp[XQPI]; + const ScalarT xqlag_dot = yp[XQLAG]; + const ScalarT pmeas_dot = yp[PMEAS]; + const ScalarT xppi_dot = yp[XPPI]; + const ScalarT pref_dot = yp[PREF_STATE]; + + const ScalarT vr = wb[0]; + const ScalarT vi = wb[1]; + + const ScalarT ir = toComponentBase(ws[IR]); + const ScalarT ii = toComponentBase(ws[II]); + const ScalarT p = toComponentBase(ws[P]); + const ScalarT q = toComponentBase(ws[Q]); + const ScalarT freq = ws[FREQ]; + const ScalarT freqref = ws[FREQREF]; + const ScalarT vref = ws[VREF]; + const ScalarT qref = toComponentBase(ws[QREF]); + const ScalarT pref_in = toComponentBase(ws[PREF_INPUT]); + + const ScalarT vldc_r = vr - Rc_ * ir + Xc_ * ii; + const ScalarT vldc_i = vi - Rc_ * ii - Xc_ * ir; + const ScalarT pfreq = droop(ef, Ddn_, Dup_); + + f[VMEAS] = -vmeas_dot + (vctrl - vmeas) / Tfltr_; + f[QMEAS] = -qmeas_dot + (q - qmeas) / Tfltr_; + f[XQPI] = -xqpi_dot + sfrz * Math::antiwindup(qpi, Ki_ * erqlim, Qmin_, Qmax_); + f[XQLAG] = -xqlag_dot + (qpi - xqlag) / Tfv_; + f[PMEAS] = -pmeas_dot + (p - pmeas) / Tp_; + f[XPPI] = -xppi_dot + Math::antiwindup(ppi, Kig_ * eplim, Pmin_, Pmax_); + f[PREF_STATE] = -pref_dot + (ppi - pref) / Tlag_; + + f[V] = -v * v + vr * vr + vi * vi; + f[VLDC] = -vldc * vldc + vldc_r * vldc_r + vldc_i * vldc_i; + f[VDROOP] = -vdroop + v + Kc_ * q; + f[VCTRL] = -vctrl + vcomp_on_ * vldc + vcomp_off_ * vdroop; + f[SFRZ] = -sfrz + Math::above(v, Vfrz_); + f[ERQ] = -erq + ref_on_ * (vref - vmeas) + ref_off_ * (qref - qmeas); + f[ERQDB] = -erqdb + Math::deadband2(erq, dbdlow_, dbdupper_); + f[ERQLIM] = -erqlim + Math::clamp(erqdb, emin_, emax_); + f[QPI] = -qpi + Math::clamp(Kp_ * erqlim + xqpi, Qmin_, Qmax_); + f[QEXT] = -Tfv_ * (qext - xqlag) + Tft_ * (qpi - xqlag); + + f[EF] = -ef + Math::deadband2(freqref - freq, fdbd1_, fdbd2_); + f[EP] = -ep + pref_in - pmeas + pfreq; + f[EPLIM] = -eplim + Math::clamp(ep, femin_, femax_); + f[PPI] = -ppi + Math::clamp(Kpg_ * eplim + xppi, Pmin_, Pmax_); + f[PEXT] = -pext + freq_on_ * pref; + + return 0; + } + + // + // Private methods + // + + /** + * @brief Smooth asymmetric frequency-droop response + * + * @param[in] error Deadbanded frequency error. + * @param[in] down Down-regulation (overfrequency) gain. + * @param[in] up Up-regulation (underfrequency) gain. + * @return Active-power frequency response. + */ + template + __attribute__((always_inline)) inline scalar_type + Repca::droop(ScalarT error, RealT down, RealT up) + { + return error * (down + (up - down) * Math::sigmoid(error)); + } + + /** + * @brief Read optional parameters from model data + * + * Omitted parameters retain their documented defaults. Real parameters + * accept real and integer values and must be finite; selectors require + * Boolean values. Loading errors are counted for verify() rather than + * thrown. + * + * @param[in] data Parameters and monitored-variable selections. + */ + template + void Repca::initializeParameters(const ModelDataT& data) + { + using Params = typename ModelDataT::Parameters; + + parameter_error_count_ = 0; + + auto load_real = [&](auto key, RealT& target, const char* name) + { + if (!data.parameters.contains(key)) + { + return; + } + + const auto& value = data.parameters.at(key); + RealT parsed_value{}; + if (const auto* real_value = std::get_if(&value)) + { + parsed_value = *real_value; + } + else if (const auto* index_value = std::get_if(&value)) + { + parsed_value = static_cast(*index_value); + } + else + { + Log::error() << "Repca: parameter '" << name << "' must be numeric\n"; + ++parameter_error_count_; + return; + } + + if (!std::isfinite(parsed_value)) + { + Log::error() << "Repca: parameter '" << name << "' must be finite\n"; + ++parameter_error_count_; + return; + } + + target = parsed_value; + }; + + auto load_switch = [&](auto key, bool& target, const char* name) + { + if (!data.parameters.contains(key)) + { + return; + } + + const auto& value = data.parameters.at(key); + if (const auto* bool_value = std::get_if(&value)) + { + target = *bool_value; + } + else + { + Log::error() << "Repca: parameter '" << name << "' must be boolean\n"; + ++parameter_error_count_; + } + }; + + load_real(Params::mva, mva_base_, "mva"); + load_switch(Params::VcompFlag, VcompFlag_, "VcompFlag"); + load_switch(Params::RefFlag, RefFlag_, "RefFlag"); + load_switch(Params::Freqflag, Freqflag_, "Freqflag"); + load_real(Params::Tfltr, Tfltr_, "Tfltr"); + load_real(Params::Vfrz, Vfrz_, "Vfrz"); + load_real(Params::Rc, Rc_, "Rc"); + load_real(Params::Xc, Xc_, "Xc"); + load_real(Params::Kc, Kc_, "Kc"); + load_real(Params::dbdlow, dbdlow_, "dbdlow"); + load_real(Params::dbdupper, dbdupper_, "dbdupper"); + load_real(Params::emax, emax_, "emax"); + load_real(Params::emin, emin_, "emin"); + load_real(Params::Kp, Kp_, "Kp"); + load_real(Params::Ki, Ki_, "Ki"); + load_real(Params::Qmax, Qmax_, "Qmax"); + load_real(Params::Qmin, Qmin_, "Qmin"); + load_real(Params::Tft, Tft_, "Tft"); + load_real(Params::Tfv, Tfv_, "Tfv"); + load_real(Params::Tp, Tp_, "Tp"); + load_real(Params::fdbd1, fdbd1_, "fdbd1"); + load_real(Params::fdbd2, fdbd2_, "fdbd2"); + load_real(Params::Ddn, Ddn_, "Ddn"); + load_real(Params::Dup, Dup_, "Dup"); + load_real(Params::femax, femax_, "femax"); + load_real(Params::femin, femin_, "femin"); + load_real(Params::Kpg, Kpg_, "Kpg"); + load_real(Params::Kig, Kig_, "Kig"); + load_real(Params::Pmax, Pmax_, "Pmax"); + load_real(Params::Pmin, Pmin_, "Pmin"); + load_real(Params::Tlag, Tlag_, "Tlag"); + + setDerivedParameters(); + } + + /** + * @brief Bind monitor selections to REPCA internal states + */ + template + void Repca::initializeMonitor() + { + using Variable = typename ModelDataT::MonitorableVariables; + + monitor_->set(Variable::qext, [this] + { return y_.getData()[static_cast(RepcaInternalVariables::QEXT)]; }); + monitor_->set(Variable::pext, [this] + { return y_.getData()[static_cast(RepcaInternalVariables::PEXT)]; }); + monitor_->set(Variable::vmeas, [this] + { return y_.getData()[static_cast(RepcaInternalVariables::VMEAS)]; }); + monitor_->set(Variable::qmeas, [this] + { return y_.getData()[static_cast(RepcaInternalVariables::QMEAS)]; }); + monitor_->set(Variable::pmeas, [this] + { return y_.getData()[static_cast(RepcaInternalVariables::PMEAS)]; }); + } + + /** + * @brief Resolve parameter-derived constants + * + * Raises the explicit controller lags in place, computes the component + * power base, and resolves selector masks. + */ + template + void Repca::setDerivedParameters() + { + // The lags are raised to the floor below, so negative values must be + // rejected here while the value as read is still available. + auto check_non_negative = [&](RealT value, const char* name) + { + if (value < ZERO) + { + Log::error() << "Repca: " << name << " must be non-negative\n"; + ++parameter_error_count_; + } + }; + + check_non_negative(Tfltr_, "Tfltr"); + check_non_negative(Tft_, "Tft"); + check_non_negative(Tfv_, "Tfv"); + check_non_negative(Tp_, "Tp"); + check_non_negative(Tlag_, "Tlag"); + + if (Tfltr_ < TIME_CONSTANT_MINIMUM || Tfv_ < TIME_CONSTANT_MINIMUM + || Tp_ < TIME_CONSTANT_MINIMUM || Tlag_ < TIME_CONSTANT_MINIMUM) + { + Log::warning() << "Repca: Tfltr, Tfv, Tp, and Tlag below " + << TIME_CONSTANT_MINIMUM + << " s are raised to that floor to keep the controller lags well posed\n"; + } + + Tfltr_ = std::max(Tfltr_, TIME_CONSTANT_MINIMUM); + Tfv_ = std::max(Tfv_, TIME_CONSTANT_MINIMUM); + Tp_ = std::max(Tp_, TIME_CONSTANT_MINIMUM); + Tlag_ = std::max(Tlag_, TIME_CONSTANT_MINIMUM); + + va_component_base_ = mva_base_ * static_cast(1.0e6); + + vcomp_on_ = VcompFlag_ ? ONE : ZERO; + vcomp_off_ = ONE - vcomp_on_; + ref_on_ = RefFlag_ ? ONE : ZERO; + ref_off_ = ONE - ref_on_; + freq_on_ = Freqflag_ ? ONE : ZERO; + } + + /** + * @brief Recover an input for a requested smooth-clamp output + * + * This initialization-only helper inverts the CommonMath smooth clamp, + * including collapsed limits. Exact-bound requests use a 0.1 offset, + * leaving less than 2e-13 at CommonMath MU = 240. + * + * @param[in] output Requested smooth-clamp output. + * @param[in] lower Lower clamp limit. + * @param[in] upper Upper clamp limit. + * @param[out] input Recovered clamp input on success. + * @return true when a finite admissible input was recovered. + * @warning Contains conditional branching and is for initialization only; + * do not use it during residual evaluation. + */ + template + bool Repca::invertClamp(ScalarT output, RealT lower, RealT upper, ScalarT& input) const + { + const RealT value = static_cast(output); + + if (!std::isfinite(value) + || !std::isfinite(lower) + || !std::isfinite(upper) + || lower > upper + || value < lower + || value > upper) + { + return false; + } + + const RealT width = upper - lower; + if (width <= INITIALIZATION_TOLERANCE) + { + input = static_cast(lower); + return true; + } + + const RealT distance_from_lower = value - lower; + const RealT distance_from_upper = upper - value; + if (distance_from_lower <= INITIALIZATION_TOLERANCE) + { + input = static_cast(lower - INITIALIZATION_LIMIT_OFFSET); + return true; + } + if (distance_from_upper <= INITIALIZATION_TOLERANCE) + { + input = static_cast(upper + INITIALIZATION_LIMIT_OFFSET); + return true; + } + + const RealT mu = Math::MU; + const RealT scaled_lower_distance = mu * distance_from_lower; + const RealT scaled_upper_distance = mu * distance_from_upper; + const RealT log_lower = logOneMinusExp(scaled_lower_distance); + const RealT log_upper = logOneMinusExp(scaled_upper_distance); + const RealT correction = (scaled_lower_distance + log_lower - log_upper) / mu; + + input = static_cast(lower + correction); + return std::isfinite(static_cast(input)); + } + + /** + * @brief Recover an input for a requested smooth deadband output + * + * @param[in] output Requested deadband output. + * @param[in] lower Lower deadband threshold. + * @param[in] upper Upper deadband threshold. + * @param[out] input Recovered deadband input on success. + * @return true when a finite input was recovered. + * @warning Contains conditional branching and is for initialization only; + * do not use it during residual evaluation. + */ + template + bool Repca::invertDeadband(ScalarT output, RealT lower, RealT upper, ScalarT& input) const + { + const RealT value = static_cast(output); + + if (!std::isfinite(value) + || !std::isfinite(lower) + || !std::isfinite(upper) + || lower > upper) + { + return false; + } + + const RealT midpoint = HALF * lower + HALF * upper; + if (std::abs(value) <= INITIALIZATION_TOLERANCE) + { + input = static_cast(midpoint); + return true; + } + + RealT lower_input = midpoint; + RealT upper_input = midpoint; + if (value < ZERO) + { + lower_input = lower + value; + } + else + { + upper_input = upper + value; + } + + const RealT lower_output = Math::deadband2(lower_input, lower, upper); + const RealT upper_output = Math::deadband2(upper_input, lower, upper); + if (!std::isfinite(lower_output) + || !std::isfinite(upper_output) + || lower_output - value > INITIALIZATION_TOLERANCE + || value - upper_output > INITIALIZATION_TOLERANCE) + { + return false; + } + + for (std::size_t iteration = 0; iteration < 128; ++iteration) + { + const RealT mid = lower_input + HALF * (upper_input - lower_input); + if (Math::deadband2(mid, lower, upper) < value) + { + lower_input = mid; + } + else + { + upper_input = mid; + } + } + + const RealT result = lower_input + HALF * (upper_input - lower_input); + input = static_cast(result); + return std::isfinite(result) + && std::abs(Math::deadband2(result, lower, upper) - value) + <= INITIALIZATION_TOLERANCE; + } + + /** + * @brief Evaluate log(1 - exp(-x)) accurately for positive x + * + * The small-x form avoids cancellation in `1 - exp(-x)`; the large-x + * form uses `log1p`. The two algebraically equivalent forms agree in + * value and first derivative at x = log(2). + * + * @param[in] x Positive argument. + * @return Numerically stable value of log(1 - exp(-x)). + */ + template + typename Repca::RealT + Repca::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 Convert a system-base power or current quantity to REPCA component base + * + * @param[in] value Quantity on the system base. + * @return The same quantity on the REPCA component base. + */ + template + [[gnu::always_inline]] inline scalar_type + Repca::toComponentBase(scalar_type value) const + { + return value * (va_system_base_ / va_component_base_); + } + + /** + * @brief Convert a component-base power quantity to system base + * + * @param[in] value Quantity on the REPCA component base. + * @return The same quantity on the system base. + */ + template + scalar_type Repca::toSystemBase(scalar_type value) const + { + return value * (va_component_base_ / va_system_base_); + } + + /** + * @brief Access the regulated-bus real voltage component + * + * @return Mutable reference to the bus real voltage state. + */ + template + scalar_type& Repca::Vr() + { + return bus_->Vr(); + } + + /** + * @brief Access the regulated-bus imaginary voltage component + * + * @return Mutable reference to the bus imaginary voltage state. + */ + template + scalar_type& Repca::Vi() + { + return bus_->Vi(); + } + + } // namespace Controller + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md index 7a8fa915f..8ae69fcb1 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 + [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 [Ieeet1](Exciter/IEEET1/README.md) | the IEEET1 exciter model diff --git a/GridKit/Model/PhasorDynamics/SystemModelData.hpp b/GridKit/Model/PhasorDynamics/SystemModelData.hpp index e82fc965e..224c95e3c 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelData.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelData.hpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -45,6 +46,7 @@ namespace GridKit using BusToSignalAdapterDataT = BusToSignalAdapterData; using BusFaultDataT = BusFaultData; using RegcaDataT = Converter::RegcaData; + using RepcaDataT = Controller::RepcaData; using Tgov1DataT = Governor::Tgov1Data; using Esdc1aDataT = Exciter::Esdc1aData; using HygovDataT = Governor::HygovData; @@ -102,6 +104,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 repca; ///< REPCA plant controllers within the model std::vector genrou; ///< GENROU instances within the model std::vector gensal; ///< GENSAL instances within the model std::vector genclassical; ///< Classical generator instances within the model diff --git a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp index 328053288..c1701687a 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 == "Repca") + { + typename SystemModelData::RepcaDataT repca; + raw_component.get_to(repca); + sm.repca.push_back(repca); + } else if (kind == "Tgov1") { typename SystemModelData::Tgov1DataT gov; diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index ceaa8ed6d..1edc3f62c 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -42,6 +42,7 @@ namespace GridKit using namespace Governor; using namespace Exciter; using namespace Stabilizer; + using namespace Controller; using namespace Converter; owns_components_ = true; @@ -482,6 +483,88 @@ namespace GridKit addComponent(stabilizer); } + // Add REPCA plant controllers after the signal producers they read at + // initialization + for (const auto& repcadata : data.repca) + { + BusT* bus = nullptr; + if (repcadata.buses.contains(RepcaBuses::bus)) + { + bus = getBus(repcadata.buses.at(RepcaBuses::bus)); + } + + auto* repca = new Repca(bus, repcadata); + + if (repcadata.signal_inputs.contains(RepcaSignalInputs::ir)) + { + const IdxT ir = repcadata.signal_inputs.at(RepcaSignalInputs::ir); + constexpr auto IR = RepcaExternalVariables::IR; + repca->getSignals().template attachSignalNode(getSignal(ir)); + } + if (repcadata.signal_inputs.contains(RepcaSignalInputs::ii)) + { + const IdxT ii = repcadata.signal_inputs.at(RepcaSignalInputs::ii); + constexpr auto II = RepcaExternalVariables::II; + repca->getSignals().template attachSignalNode(getSignal(ii)); + } + if (repcadata.signal_inputs.contains(RepcaSignalInputs::p)) + { + const IdxT p = repcadata.signal_inputs.at(RepcaSignalInputs::p); + constexpr auto P = RepcaExternalVariables::P; + repca->getSignals().template attachSignalNode

(getSignal(p)); + } + if (repcadata.signal_inputs.contains(RepcaSignalInputs::q)) + { + const IdxT q = repcadata.signal_inputs.at(RepcaSignalInputs::q); + constexpr auto Q = RepcaExternalVariables::Q; + repca->getSignals().template attachSignalNode(getSignal(q)); + } + if (repcadata.signal_inputs.contains(RepcaSignalInputs::freq)) + { + const IdxT freq = repcadata.signal_inputs.at(RepcaSignalInputs::freq); + constexpr auto FREQ = RepcaExternalVariables::FREQ; + repca->getSignals().template attachSignalNode(getSignal(freq)); + } + if (repcadata.signal_inputs.contains(RepcaSignalInputs::vref)) + { + const IdxT vref = repcadata.signal_inputs.at(RepcaSignalInputs::vref); + constexpr auto VREF = RepcaExternalVariables::VREF; + repca->getSignals().template attachSignalNode(getSignal(vref)); + } + if (repcadata.signal_inputs.contains(RepcaSignalInputs::pref)) + { + const IdxT pref = repcadata.signal_inputs.at(RepcaSignalInputs::pref); + constexpr auto PREF = RepcaExternalVariables::PREF; + repca->getSignals().template attachSignalNode(getSignal(pref)); + } + if (repcadata.signal_inputs.contains(RepcaSignalInputs::qref)) + { + const IdxT qref = repcadata.signal_inputs.at(RepcaSignalInputs::qref); + constexpr auto QREF = RepcaExternalVariables::QREF; + repca->getSignals().template attachSignalNode(getSignal(qref)); + } + if (repcadata.signal_inputs.contains(RepcaSignalInputs::freqref)) + { + const IdxT freqref = repcadata.signal_inputs.at(RepcaSignalInputs::freqref); + constexpr auto FREQREF = RepcaExternalVariables::FREQREF; + repca->getSignals().template attachSignalNode(getSignal(freqref)); + } + if (repcadata.signal_outputs.contains(RepcaSignalOutputs::qext)) + { + const IdxT qext = repcadata.signal_outputs.at(RepcaSignalOutputs::qext); + constexpr auto QEXT = RepcaInternalVariables::QEXT; + repca->getSignals().template assignSignalNode(getSignal(qext)); + } + if (repcadata.signal_outputs.contains(RepcaSignalOutputs::pext)) + { + const IdxT pext = repcadata.signal_outputs.at(RepcaSignalOutputs::pext); + constexpr auto PEXT = RepcaInternalVariables::PEXT; + repca->getSignals().template assignSignalNode(getSignal(pext)); + } + + addComponent(repca); + } + // Add constant signal sources for (const auto& srcdata : data.constant_source) { diff --git a/docs/Figures/PhasorDynamics/REPCA/diagram.png b/docs/Figures/PhasorDynamics/REPCA/diagram.png new file mode 100644 index 000000000..08fef3c15 Binary files /dev/null and b/docs/Figures/PhasorDynamics/REPCA/diagram.png differ diff --git a/docs/GridKit/Model/PhasorDynamics/Controller/README.md b/docs/GridKit/Model/PhasorDynamics/Controller/README.md new file mode 100644 index 000000000..c7b0dd611 --- /dev/null +++ b/docs/GridKit/Model/PhasorDynamics/Controller/README.md @@ -0,0 +1,14 @@ +# Controller + +```{toctree} +:maxdepth: 4 +:titlesonly: +:hidden: + +REPCA +``` + +```{include} ../../../../../GridKit/Model/PhasorDynamics/Controller/README.md +:start-line: 1 +:relative-images: +``` diff --git a/docs/GridKit/Model/PhasorDynamics/Controller/REPCA/README.md b/docs/GridKit/Model/PhasorDynamics/Controller/REPCA/README.md new file mode 100644 index 000000000..153acf4bc --- /dev/null +++ b/docs/GridKit/Model/PhasorDynamics/Controller/REPCA/README.md @@ -0,0 +1,6 @@ +# REPCA + +```{include} ../../../../../../GridKit/Model/PhasorDynamics/Controller/REPCA/README.md +:start-line: 1 +:relative-images: +``` diff --git a/docs/GridKit/Model/PhasorDynamics/README.md b/docs/GridKit/Model/PhasorDynamics/README.md index e5209139a..e867035a4 100644 --- a/docs/GridKit/Model/PhasorDynamics/README.md +++ b/docs/GridKit/Model/PhasorDynamics/README.md @@ -10,6 +10,7 @@ Branch Bus BusFault BusToSignalAdapter +Controller Converter Exciter Governor diff --git a/tests/UnitTests/PhasorDynamics/CMakeLists.txt b/tests/UnitTests/PhasorDynamics/CMakeLists.txt index 122e925ee..6c8721977 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_repca runControllerRepcaTests.cpp) +target_link_libraries( + test_phasor_controller_repca + GridKit::definitions + GridKit::phasor_dynamics_controller_repca + GridKit::phasor_dynamics_controller_repca_dependency_tracking + GridKit::phasor_dynamics_bus + GridKit::phasor_dynamics_bus_dependency_tracking + GridKit::testing) + add_executable(test_phasor_stabilizer_ieeest runStabilizerIeeestTests.cpp) target_link_libraries( test_phasor_stabilizer_ieeest @@ -187,6 +197,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 PhasorDynamicsControllerRepcaTest COMMAND test_phasor_controller_repca) add_test(NAME PhasorDynamicsStabilizerIeeestTest COMMAND test_phasor_stabilizer_ieeest) add_test(NAME PhasorDynamicsGenClassicalTest COMMAND test_phasor_gen_classical) add_test(NAME PhasorDynamicsLoadZTest COMMAND test_phasor_loadz) @@ -213,6 +224,7 @@ install( test_phasor_gensal test_phasor_exciter_sexspti test_phasor_converter_regca + test_phasor_controller_repca test_phasor_stabilizer_ieeest test_phasor_gen_classical test_phasor_system diff --git a/tests/UnitTests/PhasorDynamics/ComponentConnectionTests.hpp b/tests/UnitTests/PhasorDynamics/ComponentConnectionTests.hpp index 3b50a9f0b..a849c5c4e 100644 --- a/tests/UnitTests/PhasorDynamics/ComponentConnectionTests.hpp +++ b/tests/UnitTests/PhasorDynamics/ComponentConnectionTests.hpp @@ -4,6 +4,10 @@ #include #include +#include +#include +#include +#include #include #include #include @@ -131,6 +135,99 @@ namespace GridKit return success.report(__func__); } + + /// REGCA initializes first and publishes its branch current and power + /// to the four shared nodes. REPCA then initializes around those + /// measurements and must hold a steady state without a frequency input. + TestOutcome regcaRepca() + { + using ConverterInternal = PhasorDynamics::Converter::RegcaInternalVariables; + using ConverterParams = PhasorDynamics::Converter::RegcaParameters; + using PlantInternal = PhasorDynamics::Controller::RepcaInternalVariables; + using PlantExternal = PhasorDynamics::Controller::RepcaExternalVariables; + using PlantParams = PhasorDynamics::Controller::RepcaParameters; + + TestStatus success = true; + + PhasorDynamics::SystemModel system; + PhasorDynamics::BusInfinite bus( + static_cast(1.0), + static_cast(0.0)); + PhasorDynamics::SignalNode ir; + PhasorDynamics::SignalNode ii; + PhasorDynamics::SignalNode p; + PhasorDynamics::SignalNode q; + + PhasorDynamics::Converter::RegcaData converter_data; + converter_data.parameters[ConverterParams::p0] = static_cast(0.4); + converter_data.parameters[ConverterParams::q0] = static_cast(0.1); + 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.4); + converter_data.parameters[ConverterParams::VL1] = static_cast(0.9); + converter_data.parameters[ConverterParams::VA0] = static_cast(0.4); + converter_data.parameters[ConverterParams::VA1] = static_cast(0.9); + converter_data.parameters[ConverterParams::Vhvmax] = static_cast(1.2); + + PhasorDynamics::Controller::RepcaData plant_data; + plant_data.parameters[PlantParams::mva] = static_cast(50.0); + plant_data.parameters[PlantParams::Tp] = static_cast(0.05); + + PhasorDynamics::Converter::Regca converter(&bus, converter_data); + PhasorDynamics::Controller::Repca plant(&bus, plant_data); + + auto& converter_signals = converter.getSignals(); + converter_signals.template assignSignalNode(&ir); + converter_signals.template assignSignalNode(&ii); + converter_signals.template assignSignalNode(&p); + converter_signals.template assignSignalNode(&q); + + auto& plant_signals = plant.getSignals(); + plant_signals.template attachSignalNode(&ir); + plant_signals.template attachSignalNode(&ii); + plant_signals.template attachSignalNode(&p); + plant_signals.template attachSignalNode(&q); + + system.addBus(&bus); + system.addComponent(&converter); + system.addComponent(&plant); + + success *= system.allocate() == 0; + success *= ir.linked(); + success *= ii.linked(); + success *= p.linked(); + success *= q.linked(); + success *= system.initialize() == 0; + success *= system.evaluateResidual() == 0; + + success *= isEqual(ir.read(), static_cast(0.4), kTol); + success *= isEqual(ii.read(), static_cast(-0.1), kTol); + success *= isEqual(p.read(), static_cast(0.4), kTol); + success *= isEqual(q.read(), static_cast(0.1), kTol); + + const auto* state = plant.y().getData(); + + success *= isEqual(state[static_cast(PlantInternal::PMEAS)], + static_cast(0.8), + kTol); + success *= isEqual(state[static_cast(PlantInternal::QMEAS)], + static_cast(0.2), + kTol); + + const auto* residual = plant.getResidual().getData(); + for (IdxT row = 0; row < plant.size(); ++row) + { + success *= isEqual(residual[row], static_cast(0.0), kTol); + } + + return success.report(__func__); + } }; } // namespace Testing diff --git a/tests/UnitTests/PhasorDynamics/ControllerRepcaTests.hpp b/tests/UnitTests/PhasorDynamics/ControllerRepcaTests.hpp new file mode 100644 index 000000000..8728a23eb --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/ControllerRepcaTests.hpp @@ -0,0 +1,2361 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + using Log = ::GridKit::Utilities::Logger; + + template + class ControllerRepcaTests + { + public: + using ScalarT = scalar_type; + using IdxT = index_type; + using RealT = typename PhasorDynamics::Component::RealT; + + ControllerRepcaTests() = default; + ~ControllerRepcaTests() = default; + + static constexpr RealT kTol = + static_cast(100.0) * std::numeric_limits::epsilon(); + + /// Validate construction, defaults, parameters, signals, and time floors. + TestOutcome validation() + { + TestStatus success = true; + + noteExpectedLogs("Testing REPCA defaults, parameter floors, and invalid " + "configurations. Logged errors and warnings are expected."); + + PhasorDynamics::Bus bus(1.0, 0.0); + + PhasorDynamics::Controller::Repca empty(&bus); + success *= (empty.size() == static_cast(index(Vars::MAXIMUM))); + success *= (empty.getMonitor() == nullptr); + success *= (empty.verify() > 0); + + Fixture configured(makeData()); + configured.attachAllInputs(); + success *= (configured.repca.size() == static_cast(index(Vars::MAXIMUM))); + success *= (configured.repca.getMonitor() != nullptr); + success *= (configured.repca.verify() == 0); + + Fixture documented_defaults(makeMinimalData()); + documented_defaults.attachAllInputs(); + success *= (documented_defaults.repca.verify() == 0); + success *= defaultsMatchDocumentedValues(); + + auto integer_numeric = makeData(); + integer_numeric.parameters[Params::mva] = static_cast(100); + Fixture integer_parameter(integer_numeric); + integer_parameter.attachAllInputs(); + success *= (integer_parameter.repca.verify() == 0); + + PhasorDynamics::Controller::Repca missing_signals(&bus, makeData()); + success *= (missing_signals.verify() > 0); + + success *= invalidParameterCase(Params::mva, 0.0); + success *= invalidParameterCase(Params::Tfv, -0.1); + success *= invalidParameterCase(Params::dbdlow, 0.1); + success *= invalidParameterCase(Params::dbdupper, -0.1); + success *= invalidParameterCase(Params::emin, 0.1); + success *= invalidParameterCase(Params::emax, -0.1); + success *= invalidParameterCase(Params::Qmin, 1.1); + success *= invalidParameterCase(Params::fdbd1, 0.1); + success *= invalidParameterCase(Params::fdbd2, -0.1); + success *= invalidParameterCase(Params::Ddn, -0.1); + success *= invalidParameterCase(Params::Dup, -0.1); + success *= invalidParameterCase(Params::femin, 0.1); + success *= invalidParameterCase(Params::femax, -0.1); + success *= invalidParameterCase(Params::Pmin, 2.1); + success *= invalidParameterCase(Params::mva, true); + + success *= invalidParameterCase(Params::Tfltr, -0.2); + success *= invalidParameterCase(Params::Tft, -0.1); + success *= invalidParameterCase(Params::Tp, -0.3); + success *= invalidParameterCase(Params::Tlag, -0.4); + + const RealT nan = std::numeric_limits::quiet_NaN(); + const RealT infinity = std::numeric_limits::infinity(); + const std::array real_parameters{{ + Params::mva, + Params::Tfltr, + Params::Vfrz, + Params::Rc, + Params::Xc, + Params::Kc, + Params::dbdlow, + Params::dbdupper, + Params::emax, + Params::emin, + Params::Kp, + Params::Ki, + Params::Qmax, + Params::Qmin, + Params::Tft, + Params::Tfv, + Params::Tp, + Params::fdbd1, + Params::fdbd2, + Params::Ddn, + Params::Dup, + Params::femax, + Params::femin, + Params::Kpg, + Params::Kig, + Params::Pmax, + Params::Pmin, + Params::Tlag, + }}; + for (const Params parameter : real_parameters) + { + success *= invalidParameterCase(parameter, nan); + success *= invalidParameterCase(parameter, infinity); + success *= invalidParameterCase(parameter, -infinity); + } + success *= invalidParameterCase(Params::mva, std::numeric_limits::max()); + + { + Fixture nonfinite_system_base(makeData(), 1.0, 0.0, infinity); + nonfinite_system_base.attachAllInputs(); + success *= (nonfinite_system_base.repca.verify() > 0); + } + { + auto tiny_base_data = makeData(); + tiny_base_data.parameters[Params::mva] = std::numeric_limits::min(); + Fixture overflowing_base_ratio(tiny_base_data, + 1.0, + 0.0, + std::numeric_limits::max()); + overflowing_base_ratio.attachAllInputs(); + success *= (overflowing_base_ratio.repca.verify() > 0); + } + + const std::array flag_parameters{{ + Params::VcompFlag, + Params::RefFlag, + Params::Freqflag, + }}; + 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); + model.attachAllInputs(); + success *= (model.repca.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::Repca busless(nullptr, makeData()); + success *= (busless.verify() > 0); + + success *= unlinkedSignalRejected(); + success *= unlinkedSignalRejected(); + success *= unlinkedSignalRejected(); + success *= unlinkedSignalRejected(); + success *= unlinkedSignalRejected(); + success *= unlinkedSignalRejected(); + success *= unlinkedSignalRejected(); + success *= unlinkedSignalRejected(); + success *= unlinkedSignalRejected(); + + auto floor_data = makeInitializationData(); + floor_data.parameters[Params::Tfltr] = 0.0; + floor_data.parameters[Params::Tfv] = 0.0; + floor_data.parameters[Params::Tp] = 0.0; + floor_data.parameters[Params::Tlag] = 0.0; + + Fixture floored(floor_data); + floored.attachAllInputs(); + setInitializationInputs(floored); + success *= floored.initialize(0.25, 0.45); + success *= (floored.repca.evaluateResidual() == 0); + success *= allResidualsWithinInitTolerance(floored.repca); + + auto* y = floored.repca.y().getData(); + y[index(Vars::VMEAS)] -= 0.001; + y[index(Vars::QMEAS)] -= 0.002; + y[index(Vars::PMEAS)] -= 0.003; + y[index(Vars::PREF)] -= 0.004; + floored.repca.y().setDataUpdated(); + success *= (floored.repca.evaluateResidual() == 0); + const std::array floored_residuals{{ + {Vars::VMEAS, 1.0}, + {Vars::QMEAS, 2.0}, + {Vars::PMEAS, 3.0}, + {Vars::PREF, 4.0}, + }}; + success *= residualsMatch(floored.repca, + floored_residuals, + "floored time constants"); + + return success.report(__func__); + } + + /// Check initialization state, signal publication, monitors, and flag modes. + TestOutcome initializationAndSignals() + { + TestStatus success = true; + + noteExpectedLogs("Testing REPCA initialization without an attached frequency " + "signal. A warning is expected."); + + Fixture fixture(makeInitializationData(), 0.8, 0.6); + fixture.attachAllInputs(99.0); + setInitializationInputs(fixture); + success *= fixture.initialize(0.25, 0.45); + success *= (fixture.repca.tagDifferentiable() == 0); + success *= (fixture.repca.evaluateResidual() == 0); + + const std::array initial_state{{ + {Vars::VMEAS, 0.984002032518226}, + {Vars::QMEAS, 0.2}, + {Vars::XQPI, 0.5}, + {Vars::XQLAG, 0.5}, + {Vars::PMEAS, 0.8}, + {Vars::XPPI, 0.9}, + {Vars::PREF, 0.9}, + {Vars::V, 1.0}, + {Vars::VLDC, 0.984002032518226}, + {Vars::VDROOP, 1.08}, + {Vars::VCTRL, 0.984002032518226}, + {Vars::SFRZ, 1.0}, + {Vars::ERQ, 0.0}, + {Vars::ERQDB, 0.0}, + {Vars::ERQLIM, 0.0}, + {Vars::QPI, 0.5}, + {Vars::QEXT, 0.25}, + {Vars::EF, 0.0}, + {Vars::EP, 0.0}, + {Vars::EPLIM, 0.0}, + {Vars::PPI, 0.9}, + {Vars::PEXT, 0.45}, + }}; + success *= stateMatches(fixture.repca, initial_state, "initialization"); + + success *= scalarPreserved(fixture.input(Ext::IR), 0.2, "preserved ir"); + success *= scalarPreserved(fixture.input(Ext::II), -0.1, "preserved ii"); + success *= scalarPreserved(fixture.input(Ext::P), 0.4, "preserved p"); + success *= scalarPreserved(fixture.input(Ext::Q), 0.1, "preserved q"); + success *= scalarPreserved(fixture.input(Ext::FREQ), 0.99, "preserved freq"); + success *= scalarPreserved(fixture.qext(), 0.25, "preserved qext"); + success *= scalarPreserved(fixture.pext(), 0.45, "preserved pext"); + success *= scalarMatches(fixture.input(Ext::VREF), 0.984002032518226, "published vref"); + success *= scalarMatches(fixture.input(Ext::PREF), 0.4, "published pref"); + success *= scalarMatches(fixture.input(Ext::QREF), 0.1, "published qref"); + success *= scalarMatches(fixture.input(Ext::FREQREF), 0.99, "published freqref"); + + for (size_t row = 0; row < index(Vars::MAXIMUM); ++row) + { + const bool expected = row <= index(Vars::PREF); + if (fixture.repca.tag()[row] != expected) + { + std::cout << "REPCA differentiability tag " << row << " mismatch\n"; + success = false; + } + } + + constexpr RealT absolute_tolerance = 2.5e-7; + + success *= (fixture.repca.setAbsoluteTolerance(absolute_tolerance) == 0); + + const auto* tolerances = fixture.repca.absoluteTolerance().getData(); + for (size_t row = 0; row < index(Vars::MAXIMUM); ++row) + { + success *= valueUnchanged(tolerances[row], + absolute_tolerance, + "absolute tolerance", + row); + } + + const std::array initial_monitors{{ + 0.25, + 0.45, + 0.984002032518226, + 0.2, + 0.8, + }}; + success *= monitorMatches(fixture.repca, initial_monitors, "initialization"); + + success *= allResidualsWithinInitTolerance(fixture.repca); + + { + auto exact_data = makeInitializationData(); + exact_data.parameters[Params::mva] = 73.0; + + // This non-binary base ratio exposes any component-base round trip. + Fixture exact_commands(exact_data, 0.8, 0.6); + exact_commands.attachAllInputs(); + setInitializationInputs(exact_commands); + success *= exact_commands.initialize(0.25, 0.45); + success *= scalarPreserved(exact_commands.qext(), 0.25, "qext signal"); + success *= scalarPreserved(exact_commands.pext(), 0.45, "pext signal"); + success *= scalarPreserved(exact_commands.input(Ext::QREF), 0.1, "qref signal"); + success *= (exact_commands.repca.evaluateResidual() == 0); + success *= allResidualsWithinInitTolerance(exact_commands.repca); + } + + Fixture fallback(makeInitializationData(), 0.8, 0.6); + fallback.attachAllInputs(0.0, false); + setInitializationInputs(fallback); + success *= fallback.initialize(0.25, 0.45); + success *= scalarMatches(fallback.input(Ext::FREQREF), 1.0, "default frequency"); + success *= (fallback.repca.evaluateResidual() == 0); + success *= allResidualsWithinInitTolerance(fallback.repca); + + Fixture outputless(makeInitializationData(), + 0.8, + 0.6, + 100.0e6, + false); + outputless.attachAllInputs(); + setInitializationInputs(outputless); + success *= outputless.initialize(0.25, 0.45); + success *= (outputless.repca.evaluateResidual() == 0); + const std::array outputless_state{{ + {Vars::QEXT, 0.25}, + {Vars::PEXT, 0.45}, + }}; + success *= stateMatches(outputless.repca, + outputless_state, + "unassigned command outputs"); + success *= monitorMatches(outputless.repca, + initial_monitors, + "unassigned command outputs"); + success *= allResidualsWithinInitTolerance(outputless.repca); + + struct FlagCase + { + const char* label; + bool voltage_compensation; + bool voltage_reference; + bool frequency_control; + RealT voltage; + RealT pref; + RealT pext; + }; + + const std::array flag_cases{{ + {"droop/reactive-reference/disabled-frequency", false, false, false, 1.08, 0.8, 0.0}, + {"droop/reactive-reference/enabled-frequency", false, false, true, 1.08, 0.9, 0.45}, + {"droop/voltage-reference/disabled-frequency", false, true, false, 1.08, 0.8, 0.0}, + {"droop/voltage-reference/enabled-frequency", false, true, true, 1.08, 0.9, 0.45}, + {"line-drop/reactive-reference/disabled-frequency", true, false, false, 0.984002032518226, 0.8, 0.0}, + {"line-drop/reactive-reference/enabled-frequency", true, false, true, 0.984002032518226, 0.9, 0.45}, + {"line-drop/voltage-reference/disabled-frequency", true, true, false, 0.984002032518226, 0.8, 0.0}, + {"line-drop/voltage-reference/enabled-frequency", true, true, true, 0.984002032518226, 0.9, 0.45}, + }}; + for (const auto& test_case : flag_cases) + { + auto data = makeInitializationData(); + data.parameters[Params::VcompFlag] = test_case.voltage_compensation; + data.parameters[Params::RefFlag] = test_case.voltage_reference; + data.parameters[Params::Freqflag] = test_case.frequency_control; + + Fixture scenario(data, 0.8, 0.6); + scenario.attachAllInputs(99.0); + setInitializationInputs(scenario); + success *= scenario.initialize(0.25, 0.45); + success *= (scenario.repca.evaluateResidual() == 0); + + const std::array expected_state{{ + {Vars::VMEAS, test_case.voltage}, + {Vars::VCTRL, test_case.voltage}, + {Vars::PREF, test_case.pref}, + {Vars::PPI, test_case.pref}, + {Vars::PEXT, test_case.pext}, + }}; + success *= stateMatches(scenario.repca, + expected_state, + test_case.label); + success *= scalarMatches(scenario.qext(), 0.25, test_case.label); + success *= scalarMatches(scenario.pext(), test_case.pext, test_case.label); + success *= allResidualsWithinInitTolerance(scenario.repca); + } + + return success.report(__func__); + } + + /// Check initialization domains, adjusted limits, and atomicity. + TestOutcome initializationDomain() + { + TestStatus success = true; + + noteExpectedLogs("Testing REPCA adjusted limits and inadmissible " + "initialization points. Logged warnings and errors are expected."); + + const auto data = makeInitializationData(); + + struct RejectionCase + { + const char* label; + RealT qext; + RealT pext; + }; + + const RealT nan = std::numeric_limits::quiet_NaN(); + const RealT infinity = std::numeric_limits::infinity(); + const std::array rejection_cases{{ + {"nonfinite qext", infinity, 0.45}, + {"nonfinite pext", 0.25, infinity}, + {"nan qext", nan, 0.45}, + {"nan pext", 0.25, nan}, + }}; + for (const auto& test_case : rejection_cases) + { + success *= initializationRejectedAtomically(data, + test_case.qext, + test_case.pext, + test_case.label); + } + const std::array required_ports{{ + Ext::IR, + Ext::II, + Ext::P, + Ext::Q, + Ext::FREQ, + }}; + const std::array nonfinite_values{{infinity, -infinity, nan}}; + for (const Ext port : required_ports) + { + for (const RealT value : nonfinite_values) + { + success *= initializationRejectedAtomically(data, + 0.25, + 0.45, + "nonfinite required signal", + NonfiniteTarget::INPUT, + 0.8, + 0.6, + port, + value); + } + } + success *= initializationRejectedAtomically(data, + 0.25, + 0.45, + "nonfinite bus voltage", + NonfiniteTarget::BUS_VOLTAGE); + + auto collapsed_data = data; + + collapsed_data.parameters[Params::Qmin] = 0.5; + collapsed_data.parameters[Params::Qmax] = 0.5; + collapsed_data.parameters[Params::Pmin] = 0.9; + collapsed_data.parameters[Params::Pmax] = 0.9; + + auto reactive_aw_data = makeInitializationData(); + reactive_aw_data.parameters[Params::dbdupper] = 0.03; + reactive_aw_data.parameters[Params::emin] = 0.0; + + const std::array voltage_reference_values{{false, true}}; + const std::array, 2> voltage_cases{{ + {0.8, 0.6}, + {0.2, 0.0}, + }}; + const std::array asymmetric_reactive_state{{ + {Vars::ERQDB, -0.1}, + {Vars::ERQLIM, 0.0}, + }}; + for (const bool voltage_reference : voltage_reference_values) + { + auto data = reactive_aw_data; + data.parameters[Params::RefFlag] = voltage_reference; + for (const auto& voltage : voltage_cases) + { + Fixture asymmetric_reactive(data, + voltage.first, + voltage.second); + asymmetric_reactive.attachAllInputs(); + setInitializationInputs(asymmetric_reactive); + success *= asymmetric_reactive.initialize(0.25, 0.45); + success *= (asymmetric_reactive.repca.evaluateResidual() == 0); + success *= stateMatches(asymmetric_reactive.repca, + asymmetric_reactive_state, + "asymmetric reactive initialization"); + success *= allResidualsWithinInitTolerance(asymmetric_reactive.repca); + } + } + + auto frozen_data = reactive_aw_data; + frozen_data.parameters[Params::Vfrz] = 0.9; + + { + Fixture frozen_reactive_rate(frozen_data, 0.05, 0.0); + frozen_reactive_rate.attachAllInputs(); + setInitializationInputs(frozen_reactive_rate); + success *= frozen_reactive_rate.initialize(0.25, 0.45); + success *= (frozen_reactive_rate.repca.evaluateResidual() == 0); + success *= allResidualsWithinInitTolerance(frozen_reactive_rate.repca); + } + + { + auto active_aw_data = makeInitializationData(); + active_aw_data.parameters[Params::fdbd2] = 0.025; + active_aw_data.parameters[Params::femin] = 0.0; + Fixture asymmetric_active(active_aw_data, 0.8, 0.6); + asymmetric_active.attachAllInputs(); + setInitializationInputs(asymmetric_active); + success *= asymmetric_active.initialize(0.25, 0.45); + success *= (asymmetric_active.repca.evaluateResidual() == 0); + const std::array lower_active_state{{ + {Vars::EF, 0.0}, + {Vars::EP, -0.1}, + {Vars::EPLIM, 0.0}, + }}; + success *= stateMatches(asymmetric_active.repca, + lower_active_state, + "lower active-error boundary"); + success *= scalarMatches(asymmetric_active.input(Ext::FREQREF), + 0.995, + "asymmetric frequency reference"); + success *= scalarMatches(asymmetric_active.input(Ext::PREF), + 0.35, + "lower-bound plant reference"); + success *= allResidualsWithinInitTolerance(asymmetric_active.repca); + } + + { + auto active_aw_data = makeInitializationData(); + active_aw_data.parameters[Params::fdbd1] = -0.025; + active_aw_data.parameters[Params::femax] = 0.0; + Fixture asymmetric_active(active_aw_data, 0.8, 0.6); + asymmetric_active.attachAllInputs(); + setInitializationInputs(asymmetric_active); + success *= asymmetric_active.initialize(0.25, 0.45); + success *= (asymmetric_active.repca.evaluateResidual() == 0); + const std::array upper_active_state{{ + {Vars::EF, 0.0}, + {Vars::EP, 0.1}, + {Vars::EPLIM, 0.0}, + }}; + success *= stateMatches(asymmetric_active.repca, + upper_active_state, + "upper active-error boundary"); + success *= scalarMatches(asymmetric_active.input(Ext::FREQREF), + 0.985, + "asymmetric frequency reference"); + success *= scalarMatches(asymmetric_active.input(Ext::PREF), + 0.45, + "upper-bound plant reference"); + success *= allResidualsWithinInitTolerance(asymmetric_active.repca); + } + + auto overflow_data = makeInitializationData(); + overflow_data.parameters[Params::Rc] = std::numeric_limits::max(); + success *= initializationRejectedAtomically(overflow_data, + 0.25, + 0.45, + "nonfinite derived initialization candidate"); + + // An invalid configuration is rejected before any state is written. + { + auto invalid_data = data; + invalid_data.parameters[Params::Tfv] = -0.1; + Fixture invalid_fixture(invalid_data); + invalid_fixture.attachAllInputs(); + setInitializationInputs(invalid_fixture); + success *= (invalid_fixture.repca.allocate() == 0); + poisonState(invalid_fixture, 0.25, 0.45); + const auto invalid_y = copyVector(invalid_fixture.repca.y()); + const auto invalid_yp = copyVector(invalid_fixture.repca.yp()); + if (invalid_fixture.repca.initialize() == 0) + { + std::cout << "Expected REPCA initialization rejection: invalid configuration\n"; + success = false; + } + success *= vectorUnchanged(invalid_fixture.repca.y(), invalid_y, "state"); + success *= vectorUnchanged(invalid_fixture.repca.yp(), invalid_yp, "derivative"); + } + + // A command exactly on a limit is reconstructed through the offset + // branch of the limiter inverse, which leaves a smoothing-scaled + // residual, so only the reconstructed state is checked. + { + Fixture qmax_pmin_boundary(data, 0.8, 0.6); + qmax_pmin_boundary.attachAllInputs(); + setInitializationInputs(qmax_pmin_boundary); + success *= qmax_pmin_boundary.initialize(0.75, 0.0); + success *= (qmax_pmin_boundary.repca.evaluateResidual() == 0); + const std::array qmax_pmin_state{{ + {Vars::QPI, 1.5}, + {Vars::XQLAG, 1.5}, + {Vars::XQPI, 1.6}, + {Vars::QEXT, 0.75}, + {Vars::PREF, 0.0}, + {Vars::PPI, 0.0}, + {Vars::XPPI, -0.1}, + {Vars::PEXT, 0.0}, + }}; + success *= stateMatches(qmax_pmin_boundary.repca, + qmax_pmin_state, + "Qmax/Pmin command boundary"); + } + + { + Fixture qmin_pmax_boundary(data, 0.8, 0.6); + qmin_pmax_boundary.attachAllInputs(); + setInitializationInputs(qmin_pmax_boundary); + success *= qmin_pmax_boundary.initialize(-0.4, 1.0); + success *= (qmin_pmax_boundary.repca.evaluateResidual() == 0); + const std::array qmin_pmax_state{{ + {Vars::QPI, -0.8}, + {Vars::XQLAG, -0.8}, + {Vars::XQPI, -0.9}, + {Vars::QEXT, -0.4}, + {Vars::PREF, 2.0}, + {Vars::PPI, 2.0}, + {Vars::XPPI, 2.1}, + {Vars::PEXT, 1.0}, + }}; + success *= stateMatches(qmin_pmax_boundary.repca, + qmin_pmax_state, + "Qmin/Pmax command boundary"); + } + + { + Fixture collapsed_limits(collapsed_data, 0.8, 0.6); + collapsed_limits.attachAllInputs(); + setInitializationInputs(collapsed_limits); + success *= collapsed_limits.initialize(0.25, 0.45); + success *= (collapsed_limits.repca.evaluateResidual() == 0); + const std::array collapsed_state{{ + {Vars::XQPI, 0.5}, + {Vars::XQLAG, 0.5}, + {Vars::QPI, 0.5}, + {Vars::QEXT, 0.25}, + {Vars::XPPI, 0.9}, + {Vars::PREF, 0.9}, + {Vars::PPI, 0.9}, + {Vars::PEXT, 0.45}, + }}; + success *= stateMatches(collapsed_limits.repca, + collapsed_state, + "collapsed Q/P limits"); + success *= allResidualsWithinInitTolerance(collapsed_limits.repca); + } + + struct LimitCase + { + const char* label; + RealT qext; + RealT pext; + Vars output; + RealT expected; + }; + + const std::array limit_cases{{ + {"adjusted Qmin", -0.5, 0.45, Vars::QPI, -1.0}, + {"adjusted Qmax", 0.9, 0.45, Vars::QPI, 1.8}, + {"adjusted Pmin", 0.25, -0.1, Vars::PPI, -0.2}, + {"adjusted Pmax", 0.25, 1.1, Vars::PPI, 2.2}, + }}; + + for (const auto& test_case : limit_cases) + { + Fixture adjusted(data, 0.8, 0.6); + adjusted.attachAllInputs(); + setInitializationInputs(adjusted); + success *= adjusted.initialize(test_case.qext, test_case.pext); + success *= stateMatches(adjusted.repca, + {{test_case.output, test_case.expected}}, + test_case.label); + success *= (adjusted.repca.evaluateResidual() == 0); + success *= allResidualsWithinInitTolerance(adjusted.repca); + } + + { + auto disabled_data = data; + disabled_data.parameters[Params::Freqflag] = false; + disabled_data.parameters[Params::Pmax] = 0.5; + + Fixture disabled_frequency(disabled_data, 0.8, 0.6); + disabled_frequency.attachAllInputs(); + setInitializationInputs(disabled_frequency); + success *= disabled_frequency.initialize(0.25, 0.45); + success *= stateMatches(disabled_frequency.repca, + {{Vars::PREF, 0.8}, + {Vars::PPI, 0.8}, + {Vars::PEXT, 0.0}}, + "measured-power limit"); + success *= (disabled_frequency.repca.evaluateResidual() == 0); + success *= allResidualsWithinInitTolerance(disabled_frequency.repca); + + setState(disabled_frequency.repca, + {{Vars::PPI, 0.65}, {Vars::EPLIM, 0.1}}); + setDerivative(disabled_frequency.repca, {{Vars::XPPI, 0.0}}); + success *= (disabled_frequency.repca.evaluateResidual() == 0); + success *= residualsMatch(disabled_frequency.repca, + {{Vars::XPPI, 0.18}}, + "measured-power limit"); + } + + { + Fixture adjusted(data, 0.8, 0.6); + adjusted.attachAllInputs(); + setInitializationInputs(adjusted); + success *= adjusted.initialize(1.0, 1.25); + + setState(adjusted.repca, + {{Vars::QPI, 1.75}, + {Vars::ERQLIM, 0.1}, + {Vars::SFRZ, 1.0}, + {Vars::PPI, 2.25}, + {Vars::EPLIM, 0.1}}); + setDerivative(adjusted.repca, {{Vars::XQPI, 0.0}, {Vars::XPPI, 0.0}}); + success *= (adjusted.repca.evaluateResidual() == 0); + success *= residualsMatch(adjusted.repca, + {{Vars::XQPI, 0.3}, {Vars::XPPI, 0.18}}, + "adjusted antiwindup limits"); + + setState(adjusted.repca, + {{Vars::QPI, 0.0}, + {Vars::ERQLIM, 0.0}, + {Vars::XQPI, 1.75}, + {Vars::PPI, 0.0}, + {Vars::EPLIM, 0.0}, + {Vars::XPPI, 2.25}}); + success *= (adjusted.repca.evaluateResidual() == 0); + success *= residualsMatch(adjusted.repca, + {{Vars::QPI, 1.75}, {Vars::PPI, 2.25}}, + "adjusted command limits"); + } + + return success.report(__func__); + } + + /// Check every residual row against an independent numerical answer key. + /// The expected values are literals, not a second implementation of REPCA. + TestOutcome residualEquations() + { + TestStatus success = true; + + Fixture fixture(makeResidualData(), kStateVr, kStateVi); + fixture.attachAllInputs(); + setAnswerKeyInputs(fixture); + success *= fixture.prepare(0.0, 0.0); + setAnswerKeyState(fixture.repca); + success *= (fixture.repca.evaluateResidual() == 0); + + const std::array expected_residuals{{ + {Vars::VMEAS, 0.4}, + {Vars::QMEAS, 0.45}, + {Vars::XQPI, 0.3}, + {Vars::XQLAG, 0.46}, + {Vars::PMEAS, 0.5}, + {Vars::XPPI, -0.3}, + {Vars::PREF, 0.4}, + {Vars::V, -1.28}, + {Vars::VLDC, 0.028}, + {Vars::VDROOP, 0.1}, + {Vars::VCTRL, 0.05}, + {Vars::SFRZ, 0.5}, + {Vars::ERQ, -0.63}, + {Vars::ERQDB, 0.75}, + {Vars::ERQLIM, -0.35}, + {Vars::QPI, -0.05}, + {Vars::QEXT, -1.345}, + {Vars::EF, -0.015}, + {Vars::EP, -0.4}, + {Vars::EPLIM, 1.1}, + {Vars::PPI, 0.1}, + {Vars::PEXT, 0.05}, + }}; + + success *= (static_cast(fixture.repca.getResidual().getSize()) + == expected_residuals.size()); + const auto* residual = fixture.repca.getResidual().getData(); + for (size_t row = 0; row < expected_residuals.size(); ++row) + { + const auto variable = expected_residuals[row].variable; + if (index(variable) != row) + { + std::cout << "REPCA residual key position " << row << " names row " + << variableName(variable) << '\n'; + success = false; + } + success *= scalarMatches(residual[index(variable)], + expected_residuals[row].value, + variableName(variable)); + } + + return success.report(__func__); + } + + /// Check reactive modes, smooth limits, antiwindup, and lead-lag behavior. + TestOutcome reactiveControl() + { + TestStatus success = true; + + struct FlagCase + { + const char* label; + bool voltage_compensation; + bool voltage_reference; + RealT vctrl; + RealT erq; + }; + + const std::array flag_cases{{ + {"droop/reactive-reference", false, false, 0.08, 0.30}, + {"droop/voltage-reference", false, true, 0.08, 0.10}, + {"line-drop/reactive-reference", true, false, -0.08, 0.30}, + {"line-drop/voltage-reference", true, true, -0.08, 0.10}, + }}; + for (const auto& test_case : flag_cases) + { + auto data = makeResidualData(); + data.parameters[Params::VcompFlag] = test_case.voltage_compensation; + data.parameters[Params::RefFlag] = test_case.voltage_reference; + + Fixture fixture(data); + fixture.attachAllInputs(); + setAnswerKeyInputs(fixture); + fixture.input(Ext::VREF) = 1.05; + fixture.input(Ext::QREF) = 0.20; + + success *= fixture.prepare(0.0, 0.0); + setState(fixture.repca, + {{Vars::VMEAS, 0.95}, + {Vars::QMEAS, 0.10}, + {Vars::VLDC, 0.92}, + {Vars::VDROOP, 1.08}, + {Vars::VCTRL, 1.0}, + {Vars::ERQ, 0.0}}); + success *= (fixture.repca.evaluateResidual() == 0); + + const std::array expected_residuals{{ + {Vars::VCTRL, test_case.vctrl}, + {Vars::ERQ, test_case.erq}, + }}; + success *= residualsMatch(fixture.repca, + expected_residuals, + test_case.label); + } + + Fixture fixture(makeResidualData()); + fixture.attachAllInputs(); + setAnswerKeyInputs(fixture); + success *= fixture.prepare(0.0, 0.0); + + // A voltage of zero must clear the freeze threshold by the same + // margin the enabled probe clears it, so the threshold is raised. + { + auto freeze_data = makeResidualData(); + freeze_data.parameters[Params::Vfrz] = 0.8; + + Fixture freeze(freeze_data); + freeze.attachAllInputs(); + setAnswerKeyInputs(freeze); + success *= freeze.prepare(0.0, 0.0); + + const std::array freeze_cases{{ + {0.0, 0.0}, + {1.6, 1.0}, + }}; + for (const auto& test_case : freeze_cases) + { + setState(freeze.repca, {{Vars::V, test_case.input}, {Vars::SFRZ, 0.0}}); + success *= (freeze.repca.evaluateResidual() == 0); + success *= residualsMatch(freeze.repca, + {{Vars::SFRZ, test_case.expected}}, + "freeze gate"); + } + } + + // The interior probe sits at the midpoint of the band, where the + // smooth deadband cancels exactly. + const std::array deadband_cases{{ + {-0.82, -0.8}, + {0.005, 0.0}, + {0.83, 0.8}, + }}; + for (const auto& test_case : deadband_cases) + { + setState(fixture.repca, {{Vars::ERQ, test_case.input}, {Vars::ERQDB, 0.0}}); + success *= (fixture.repca.evaluateResidual() == 0); + success *= residualsMatch(fixture.repca, + {{Vars::ERQDB, test_case.expected}}, + "reactive-power deadband"); + } + + const std::array error_limit_cases{{ + {-1.5, -0.7}, + {0.05, 0.05}, + {1.6, 0.8}, + }}; + for (const auto& test_case : error_limit_cases) + { + setState(fixture.repca, {{Vars::ERQDB, test_case.input}, {Vars::ERQLIM, 0.0}}); + success *= (fixture.repca.evaluateResidual() == 0); + success *= residualsMatch(fixture.repca, + {{Vars::ERQLIM, test_case.expected}}, + "reactive-power error limit"); + } + + const std::array command_limit_cases{{ + {-1.6, -0.8}, + {0.05, 0.05}, + {1.7, 0.9}, + }}; + for (const auto& test_case : command_limit_cases) + { + setState(fixture.repca, + {{Vars::XQPI, test_case.input}, + {Vars::ERQLIM, 0.0}, + {Vars::QPI, 0.0}}); + success *= (fixture.repca.evaluateResidual() == 0); + success *= residualsMatch(fixture.repca, + {{Vars::QPI, test_case.expected}}, + "reactive-power command limit"); + } + + // 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{{ + {-1.6, -0.4, 0.0}, + {-1.6, 0.4, 1.2}, + {0.05, -0.4, -1.2}, + {0.05, 0.4, 1.2}, + {1.7, -0.4, -1.2}, + {1.7, 0.4, 0.0}, + }}; + for (const auto& test_case : antiwindup_cases) + { + setState(fixture.repca, + {{Vars::QPI, test_case.output}, + {Vars::ERQLIM, test_case.error}, + {Vars::SFRZ, 1.0}}); + setDerivative(fixture.repca, {{Vars::XQPI, 0.0}}); + success *= (fixture.repca.evaluateResidual() == 0); + success *= residualsMatch(fixture.repca, + {{Vars::XQPI, test_case.expected}}, + "reactive-power antiwindup"); + } + + setState(fixture.repca, + {{Vars::XQLAG, 0.14}, + {Vars::QPI, 0.27}, + {Vars::QEXT, 0.20}}); + setDerivative(fixture.repca, {{Vars::XQLAG, -0.04}}); + success *= (fixture.repca.evaluateResidual() == 0); + const std::array lead_lag_residuals{{ + {Vars::XQLAG, 0.092}, + {Vars::QEXT, -0.624}, + }}; + success *= residualsMatch(fixture.repca, + lead_lag_residuals, + "reactive-command lead-lag"); + + // The command sits beyond Qmax with the error driving further out, + // so the blocked gate leaves the PI state with no sensitivity to the + // gate, the error, or the command. + { + Fixture blocked(makeResidualData()); + blocked.attachAllInputs(); + setAnswerKeyInputs(blocked); + success *= blocked.prepare(0.0, 0.0); + setState(blocked.repca, + {{Vars::QPI, 1.7}, {Vars::ERQLIM, 0.4}, {Vars::SFRZ, 1.0}}); + setDerivative(blocked.repca, {{Vars::XQPI, 0.0}}); + numberVariables(blocked, 1.0); + success *= (blocked.repca.evaluateResidual() == 0); + + const DependencyTracking::Variable::DependencyMap expected{ + {index(Vars::XQPI), -1.0}, + {index(Vars::SFRZ), 0.0}, + {index(Vars::ERQLIM), 0.0}, + {index(Vars::QPI), 0.0}, + }; + success *= jacobianRowMatches( + blocked.repca.getResidual().getData()[index(Vars::XQPI)].getDependencies(), + expected, + index(Vars::XQPI), + "blocked reactive-power antiwindup", + kTol); + } + + return success.report(__func__); + } + + /// Check active-power modes, smooth limits, antiwindup, and command lag. + TestOutcome activePowerControl() + { + TestStatus success = true; + + struct FlagCase + { + const char* label; + bool frequency_control; + RealT pext; + }; + + const std::array flag_cases{{ + {"disabled frequency control", false, -0.6}, + {"enabled frequency control", true, 0.2}, + }}; + for (const auto& test_case : flag_cases) + { + auto data = makeResidualData(); + data.parameters[Params::Freqflag] = test_case.frequency_control; + Fixture fixture(data); + fixture.attachAllInputs(); + setAnswerKeyInputs(fixture); + success *= fixture.prepare(0.0, 0.0); + setState(fixture.repca, {{Vars::PREF, 0.8}, {Vars::PEXT, 0.3}}); + success *= (fixture.repca.evaluateResidual() == 0); + success *= residualsMatch(fixture.repca, + {{Vars::PEXT, test_case.pext}}, + test_case.label); + } + + Fixture fixture(makeResidualData()); + fixture.attachAllInputs(); + setAnswerKeyInputs(fixture); + success *= fixture.prepare(0.0, 0.0); + + // The interior probe sits at the midpoint of the band, where the + // smooth deadband cancels exactly. + const std::array frequency_deadband_cases{{ + {-0.81, -0.8}, + {0.0025, 0.0}, + {0.815, 0.8}, + }}; + for (const auto& test_case : frequency_deadband_cases) + { + fixture.input(Ext::FREQ) = 1.0; + fixture.input(Ext::FREQREF) = 1.0 + test_case.input; + setState(fixture.repca, {{Vars::EF, 0.0}}); + success *= (fixture.repca.evaluateResidual() == 0); + success *= residualsMatch(fixture.repca, + {{Vars::EF, test_case.expected}}, + "frequency deadband"); + } + + const std::array droop_cases{{ + {-0.9, -1.8}, + {0.0, 0.0}, + {0.9, 0.9}, + }}; + fixture.input(Ext::PREF) = 0.2; + for (const auto& test_case : droop_cases) + { + setState(fixture.repca, + {{Vars::EF, test_case.input}, + {Vars::EP, 0.0}, + {Vars::PMEAS, 0.4}}); + success *= (fixture.repca.evaluateResidual() == 0); + success *= residualsMatch(fixture.repca, + {{Vars::EP, test_case.expected}}, + "frequency droop"); + } + + const std::array error_limit_cases{{ + {-1.3, -0.5}, + {0.05, 0.05}, + {1.4, 0.6}, + }}; + for (const auto& test_case : error_limit_cases) + { + setState(fixture.repca, {{Vars::EP, test_case.input}, {Vars::EPLIM, 0.0}}); + success *= (fixture.repca.evaluateResidual() == 0); + success *= residualsMatch(fixture.repca, + {{Vars::EPLIM, test_case.expected}}, + "active-power error limit"); + } + + const std::array command_limit_cases{{ + {-0.8, 0.0}, + {1.0, 1.0}, + {2.8, 2.0}, + }}; + for (const auto& test_case : command_limit_cases) + { + setState(fixture.repca, + {{Vars::XPPI, test_case.input}, + {Vars::EPLIM, 0.0}, + {Vars::PPI, 0.0}}); + success *= (fixture.repca.evaluateResidual() == 0); + success *= residualsMatch(fixture.repca, + {{Vars::PPI, test_case.expected}}, + "active-power command limit"); + } + + // 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{{ + {-0.8, -0.5, 0.0}, + {-0.8, 0.5, 0.9}, + {1.0, -0.5, -0.9}, + {1.0, 0.5, 0.9}, + {2.8, -0.5, -0.9}, + {2.8, 0.5, 0.0}, + }}; + for (const auto& test_case : antiwindup_cases) + { + setState(fixture.repca, + {{Vars::PPI, test_case.output}, + {Vars::EPLIM, test_case.error}}); + setDerivative(fixture.repca, {{Vars::XPPI, 0.0}}); + success *= (fixture.repca.evaluateResidual() == 0); + success *= residualsMatch(fixture.repca, + {{Vars::XPPI, test_case.expected}}, + "active-power antiwindup"); + } + + setState(fixture.repca, {{Vars::PPI, 0.66}, {Vars::PREF, 0.60}}); + setDerivative(fixture.repca, {{Vars::PREF, 0.05}}); + success *= (fixture.repca.evaluateResidual() == 0); + success *= residualsMatch(fixture.repca, + {{Vars::PREF, 0.07}}, + "active-power command lag"); + + return success.report(__func__); + } + + /// Check every differential residual with nonzero explicit derivatives. + TestOutcome derivatives() + { + TestStatus success = true; + + Fixture fixture(makeResidualData(), kStateVr, kStateVi); + fixture.attachAllInputs(); + setAnswerKeyInputs(fixture); + success *= fixture.prepare(0.0, 0.0); + setAnswerKeyState(fixture.repca); + setDerivative(fixture.repca, + {{Vars::VMEAS, 0.2}, + {Vars::QMEAS, -0.1}, + {Vars::XQPI, 0.4}, + {Vars::XQLAG, -0.3}, + {Vars::PMEAS, 0.6}, + {Vars::XPPI, -0.5}, + {Vars::PREF, 0.8}}); + success *= (fixture.repca.evaluateResidual() == 0); + const std::array expected_residuals{{ + {Vars::VMEAS, 0.3}, + {Vars::QMEAS, 0.35}, + {Vars::XQPI, 0.2}, + {Vars::XQLAG, 0.36}, + {Vars::PMEAS, 0.4}, + {Vars::XPPI, -0.4}, + {Vars::PREF, 0.3}, + }}; + success *= residualsMatch(fixture.repca, + expected_residuals, + "explicit derivatives"); + + return success.report(__func__); + } + + /// Check every dependency-tracking Jacobian row against an independent + /// numerical and structural answer key, with both selector settings and + /// a non-unit alpha. + TestOutcome dependencyTracking() + { + TestStatus success = true; + + const auto data = makeResidualData(); + const auto dependency = dependencyTrackingJacobian(data, success); + + success *= jacobianMatches(dependency, + expectedJacobian(), + "dependency tracking", + kTol); + + auto all_flags_off_data = data; + all_flags_off_data.parameters[Params::VcompFlag] = false; + all_flags_off_data.parameters[Params::RefFlag] = false; + all_flags_off_data.parameters[Params::Freqflag] = false; + const auto all_flags_off = + dependencyTrackingJacobian(all_flags_off_data, success); + success *= jacobianMatches(all_flags_off, + expectedJacobianAllFlagsOff(), + "all-flags-off dependency tracking", + kTol); + + const auto nonunit_alpha_dependency = + dependencyTrackingJacobian(data, success, kNonunitAlpha); + success *= jacobianMatches(nonunit_alpha_dependency, + expectedJacobianNonunitAlpha(), + "non-unit-alpha dependency tracking", + kTol); + + return success.report(__func__); + } + +#ifdef GRIDKIT_ENABLE_ENZYME + /// One rich state, both selector settings, and a non-unit alpha drive both + /// sensitivity paths; every Enzyme CSR row must match dependency tracking. + TestOutcome jacobian() + { + TestStatus success = true; + + const auto data = makeResidualData(); + + success *= jacobianMatches(enzymeJacobian(data, success), + dependencyTrackingJacobian(data, success), + "Enzyme versus dependency tracking", + kTol); + + auto all_flags_off_data = data; + all_flags_off_data.parameters[Params::VcompFlag] = false; + all_flags_off_data.parameters[Params::RefFlag] = false; + all_flags_off_data.parameters[Params::Freqflag] = false; + success *= jacobianMatches( + enzymeJacobian(all_flags_off_data, success), + dependencyTrackingJacobian(all_flags_off_data, success), + "all-flags-off Enzyme versus dependency tracking", + kTol); + + success *= jacobianMatches( + enzymeJacobian(data, success, kNonunitAlpha), + dependencyTrackingJacobian(data, success, kNonunitAlpha), + "non-unit-alpha Enzyme versus dependency tracking", + kTol); + + return success.report(__func__); + } +#endif + + private: + using Params = PhasorDynamics::Controller::RepcaParameters; + using Vars = PhasorDynamics::Controller::RepcaInternalVariables; + using Ext = PhasorDynamics::Controller::RepcaExternalVariables; + using Mon = PhasorDynamics::Controller::RepcaMonitorableVariables; + using Data = PhasorDynamics::Controller::RepcaData; + using RepcaT = PhasorDynamics::Controller::Repca; + + 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 output; + RealT error; + RealT expected; + }; + + enum class NonfiniteTarget + { + NONE, + INPUT, + BUS_VOLTAGE + }; + + /// Owns the regulated bus, REPCA, assigned command nodes, and attached + /// input nodes. Signal storage precedes the model so every referenced node + /// outlives REPCA; 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 qext_node_; + PhasorDynamics::SignalNode pext_node_; + + public: + explicit Fixture(const Data& data, + RealT vr = 1.0, + RealT vi = 0.0, + RealT system_va_base = 100.0e6, + bool assign_command_outputs = true) + : bus(static_cast(vr), static_cast(vi)), + repca(&bus, data) + { + repca.setSystemBase(60.0, system_va_base); + if (assign_command_outputs) + { + repca.getSignals().template assignSignalNode(&qext_node_); + repca.getSignals().template assignSignalNode(&pext_node_); + } + } + + Fixture(const Fixture&) = delete; + Fixture& operator=(const Fixture&) = delete; + + void attachRequiredInputs(RealT initial_value = 0.0) + { + const IdxT external_index_base = repca.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 = repca.getSignals(); + signals.template attachSignalNode(&input_nodes_[index(Ext::IR)]); + signals.template attachSignalNode(&input_nodes_[index(Ext::II)]); + signals.template attachSignalNode(&input_nodes_[index(Ext::P)]); + signals.template attachSignalNode(&input_nodes_[index(Ext::Q)]); + } + + void attachAllInputs(RealT initial_value = 0.0, + bool attach_frequency = true) + { + attachRequiredInputs(initial_value); + + auto& signals = repca.getSignals(); + if (attach_frequency) + { + signals.template attachSignalNode(&input_nodes_[index(Ext::FREQ)]); + } + signals.template attachSignalNode(&input_nodes_[index(Ext::VREF)]); + signals.template attachSignalNode(&input_nodes_[index(Ext::PREF)]); + signals.template attachSignalNode(&input_nodes_[index(Ext::QREF)]); + signals.template attachSignalNode(&input_nodes_[index(Ext::FREQREF)]); + } + + void setCommands(RealT qext, RealT pext) + { + auto* y = repca.y().getData(); + y[index(Vars::QEXT)] = static_cast(qext); + y[index(Vars::PEXT)] = static_cast(pext); + repca.y().setDataUpdated(); + } + + /// Arrange the allocation, verification, bus, and command prerequisites. + bool prepare(RealT qext, RealT pext) + { + const bool success = (bus.allocate() == 0) && (repca.allocate() == 0) + && (repca.verify() == 0) && (bus.initialize() == 0); + if (!success) + { + std::cout << "REPCA fixture preparation failed\n"; + return false; + } + setCommands(qext, pext); + return true; + } + + bool initialize(RealT qext, RealT pext) + { + if (!prepare(qext, pext)) + { + return false; + } + if (repca.initialize() != 0) + { + std::cout << "REPCA initialization failed\n"; + return false; + } + return true; + } + + T qext() const + { + return repca.y().getData()[index(Vars::QEXT)]; + } + + T pext() const + { + return repca.y().getData()[index(Vars::PEXT)]; + } + + T& input(Ext port) + { + return input_values_[index(port)]; + } + + IdxT inputIndex(Ext port) const + { + return input_indices_[index(port)]; + } + + PhasorDynamics::Bus bus; + PhasorDynamics::Controller::Repca repca; + }; + + static constexpr RealT kStateVr = 0.9; + static constexpr RealT kStateVi = 0.4; + static constexpr RealT kNonunitAlpha = 2.5; + + 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 externalColumn(size_t port) + { + return kExternalColumnBase + port; + } + + Data makeMinimalData() const + { + Data data; + data.device_class = "Repca"; + data.disambiguation_string = "repca_test"; + data.monitored_variables.insert(Mon::qext); + data.monitored_variables.insert(Mon::pext); + data.monitored_variables.insert(Mon::vmeas); + data.monitored_variables.insert(Mon::qmeas); + data.monitored_variables.insert(Mon::pmeas); + return data; + } + + Data makeExplicitDefaultData() const + { + auto data = makeMinimalData(); + data.parameters[Params::mva] = 100.0; + data.parameters[Params::VcompFlag] = true; + data.parameters[Params::RefFlag] = true; + data.parameters[Params::Freqflag] = false; + data.parameters[Params::Tfltr] = 0.05; + data.parameters[Params::Vfrz] = 0.7; + data.parameters[Params::Rc] = 0.0; + data.parameters[Params::Xc] = 0.0; + data.parameters[Params::Kc] = 1.0; + data.parameters[Params::dbdlow] = 0.0; + data.parameters[Params::dbdupper] = 0.0; + data.parameters[Params::emax] = 1.0; + data.parameters[Params::emin] = -1.0; + data.parameters[Params::Kp] = 10.0; + data.parameters[Params::Ki] = 10.0; + data.parameters[Params::Qmax] = 1.0; + data.parameters[Params::Qmin] = -1.0; + data.parameters[Params::Tft] = 0.0; + data.parameters[Params::Tfv] = 3.0; + data.parameters[Params::Tp] = 0.0; + data.parameters[Params::fdbd1] = 0.0; + data.parameters[Params::fdbd2] = 0.0; + data.parameters[Params::Ddn] = 20.0; + data.parameters[Params::Dup] = 0.0; + data.parameters[Params::femax] = 1.0; + data.parameters[Params::femin] = -1.0; + data.parameters[Params::Kpg] = 10.0; + data.parameters[Params::Kig] = 10.0; + data.parameters[Params::Pmax] = 2.0; + data.parameters[Params::Pmin] = 0.0; + data.parameters[Params::Tlag] = 3.0; + return data; + } + + Data makeData() const + { + auto data = makeExplicitDefaultData(); + data.parameters[Params::Freqflag] = true; + data.parameters[Params::Tp] = 0.05; + return data; + } + + /// Distinct nonzero values for every parameter. The limiter 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::mva] = 50.0; + data.parameters[Params::Tfltr] = 0.2; + data.parameters[Params::Rc] = 0.02; + data.parameters[Params::Xc] = 0.03; + data.parameters[Params::Kc] = 0.4; + data.parameters[Params::dbdlow] = -0.02; + data.parameters[Params::dbdupper] = 0.03; + data.parameters[Params::emax] = 0.8; + data.parameters[Params::emin] = -0.7; + data.parameters[Params::Kp] = 2.0; + data.parameters[Params::Ki] = 3.0; + data.parameters[Params::Qmax] = 0.9; + data.parameters[Params::Qmin] = -0.8; + data.parameters[Params::Tft] = 0.2; + data.parameters[Params::Tfv] = 2.5; + data.parameters[Params::Tp] = 0.4; + data.parameters[Params::fdbd1] = -0.01; + data.parameters[Params::fdbd2] = 0.015; + data.parameters[Params::Ddn] = 2.0; + data.parameters[Params::Dup] = 1.0; + data.parameters[Params::femax] = 0.6; + data.parameters[Params::femin] = -0.5; + data.parameters[Params::Kpg] = 1.7; + data.parameters[Params::Kig] = 1.8; + data.parameters[Params::Pmax] = 2.0; + data.parameters[Params::Tlag] = 0.5; + return data; + } + + /// Both deadbands and both error limits are symmetric and the droop + /// gains are equal, so an operating point with no error reconstructs + /// exactly; the commands and the freeze threshold clear their limits. + Data makeInitializationData() const + { + auto data = makeResidualData(); + data.parameters[Params::Vfrz] = 0.2; + data.parameters[Params::dbdupper] = 0.02; + data.parameters[Params::emin] = -0.8; + data.parameters[Params::Qmax] = 1.5; + data.parameters[Params::fdbd1] = -0.015; + data.parameters[Params::Dup] = 2.0; + data.parameters[Params::femin] = -0.6; + return data; + } + + template + void setInitializationInputs(Fixture& fixture) const + { + fixture.input(Ext::IR) = static_cast(0.2); + fixture.input(Ext::II) = static_cast(-0.1); + fixture.input(Ext::P) = static_cast(0.4); + fixture.input(Ext::Q) = static_cast(0.1); + fixture.input(Ext::FREQ) = static_cast(0.99); + } + + template + void setAnswerKeyInputs(Fixture& fixture) const + { + fixture.input(Ext::IR) = static_cast(1.0); + fixture.input(Ext::II) = static_cast(2.0); + fixture.input(Ext::P) = static_cast(0.35); + fixture.input(Ext::Q) = static_cast(0.25); + fixture.input(Ext::FREQ) = static_cast(0.2); + fixture.input(Ext::VREF) = static_cast(1.05); + fixture.input(Ext::PREF) = static_cast(0.55); + fixture.input(Ext::QREF) = static_cast(0.3); + fixture.input(Ext::FREQREF) = static_cast(1.0); + } + + template + void setAnswerKeyState(PhasorDynamics::Controller::Repca& repca) const + { + // Every smooth-transition argument keeps a saturation margin, and + // every clamp that must pass its input through sits at the midpoint + // of its limits, so each row carries its ideal value. + setState(repca, + {{Vars::VMEAS, 0.85}, + {Vars::QMEAS, 0.45}, + {Vars::XQPI, -0.75}, + {Vars::XQLAG, -0.05}, + {Vars::PMEAS, 0.3}, + {Vars::XPPI, 1.85}, + {Vars::PREF, 0.35}, + {Vars::V, 1.5}, + {Vars::VLDC, 1.0}, + {Vars::VDROOP, 1.6}, + {Vars::VCTRL, 0.95}, + {Vars::SFRZ, 0.5}, + {Vars::ERQ, 0.83}, + {Vars::ERQDB, 0.05}, + {Vars::ERQLIM, 0.4}, + {Vars::QPI, 0.1}, + {Vars::QEXT, 0.25}, + {Vars::EF, 0.8}, + {Vars::EP, 2.0}, + {Vars::EPLIM, -0.5}, + {Vars::PPI, 0.9}, + {Vars::PEXT, 0.15}}); + setDerivative(repca, + {{Vars::VMEAS, 0.1}, + {Vars::QMEAS, -0.2}, + {Vars::XQPI, 0.3}, + {Vars::XQLAG, -0.4}, + {Vars::PMEAS, 0.5}, + {Vars::XPPI, -0.6}, + {Vars::PREF, 0.7}}); + } + + bool defaultsMatchDocumentedValues() const + { + Fixture implicit_defaults(makeMinimalData(), 0.9, 0.4); + Fixture explicit_defaults(makeExplicitDefaultData(), 0.9, 0.4); + implicit_defaults.attachAllInputs(); + explicit_defaults.attachAllInputs(); + + implicit_defaults.input(Ext::P) = 0.2; + implicit_defaults.input(Ext::Q) = 0.1; + implicit_defaults.input(Ext::FREQ) = 1.0; + explicit_defaults.input(Ext::P) = 0.2; + explicit_defaults.input(Ext::Q) = 0.1; + explicit_defaults.input(Ext::FREQ) = 1.0; + + bool success = implicit_defaults.initialize(0.1, 0.2) + && explicit_defaults.initialize(0.1, 0.2); + if (!success) + { + std::cout << "REPCA documented-default comparison failed to initialize\n"; + return false; + } + + if (implicit_defaults.repca.evaluateResidual() != 0) + { + success = false; + } + if (explicit_defaults.repca.evaluateResidual() != 0) + { + success = false; + } + if (!vectorsMatch(implicit_defaults.repca.y(), + explicit_defaults.repca.y(), + "documented-default state")) + { + success = false; + } + if (!vectorsMatch(implicit_defaults.repca.yp(), + explicit_defaults.repca.yp(), + "documented-default derivative")) + { + success = false; + } + if (!vectorsMatch(implicit_defaults.repca.getResidual(), + explicit_defaults.repca.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.repca); + setAnswerKeyState(explicit_defaults.repca); + if (implicit_defaults.repca.evaluateResidual() != 0) + { + success = false; + } + if (explicit_defaults.repca.evaluateResidual() != 0) + { + success = false; + } + if (!vectorsMatch(implicit_defaults.repca.getResidual(), + explicit_defaults.repca.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); + fixture.attachAllInputs(); + return fixture.repca.verify() > 0; + } + + template + bool unlinkedSignalRejected() const + { + Fixture fixture(makeData()); + fixture.attachAllInputs(); + PhasorDynamics::SignalNode unlinked_node; + fixture.repca.getSignals().template attachSignalNode(&unlinked_node); + return fixture.repca.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 qext, RealT pext) const + { + auto* y = fixture.repca.y().getData(); + auto* yp = fixture.repca.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(qext, pext); + fixture.repca.yp().setDataUpdated(); + } + + bool initializationRejectedAtomically(const Data& data, + RealT qext, + RealT pext, + const char* label, + NonfiniteTarget target = NonfiniteTarget::NONE, + RealT initial_vr = 0.8, + RealT initial_vi = 0.6, + Ext poisoned_port = Ext::FREQ, + RealT poison_value = + std::numeric_limits::infinity()) const + { + Fixture fixture(data, initial_vr, initial_vi); + fixture.attachAllInputs(77.0); + setInitializationInputs(fixture); + if (!fixture.prepare(qext, pext)) + { + return false; + } + + if (target == NonfiniteTarget::INPUT) + { + fixture.input(poisoned_port) = poison_value; + } + if (target == NonfiniteTarget::BUS_VOLTAGE) + { + fixture.bus.Vr() = poison_value; + fixture.bus.y().setDataUpdated(); + } + + poisonState(fixture, qext, pext); + + const auto y_before = copyVector(fixture.repca.y()); + const auto yp_before = copyVector(fixture.repca.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.repca.initialize() == 0) + { + std::cout << "Expected REPCA initialization rejection: " << label << '\n'; + success = false; + } + + if (!scalarPreserved(fixture.qext(), qext, "rejected qext preservation")) + { + success = false; + } + if (!scalarPreserved(fixture.pext(), pext, "rejected pext preservation")) + { + success = false; + } + if (!vectorUnchanged(fixture.repca.y(), y_before, "state")) + { + success = false; + } + if (!vectorUnchanged(fixture.repca.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::Repca& repca, + std::initializer_list values) const + { + auto* y = repca.y().getData(); + for (const auto& [variable, value] : values) + { + y[index(variable)] = static_cast(value); + } + repca.y().setDataUpdated(); + } + + template + void setDerivative(PhasorDynamics::Controller::Repca& repca, + std::initializer_list values) const + { + auto* yp = repca.yp().getData(); + for (const auto& [variable, value] : values) + { + yp[index(variable)] = static_cast(value); + } + repca.yp().setDataUpdated(); + } + + static const char* variableName(Vars variable) + { + static constexpr std::array names{{ + "VMEAS", + "QMEAS", + "XQPI", + "XQLAG", + "PMEAS", + "XPPI", + "PREF", + "V", + "VLDC", + "VDROOP", + "VCTRL", + "SFRZ", + "ERQ", + "ERQDB", + "ERQLIM", + "QPI", + "QEXT", + "EF", + "EP", + "EPLIM", + "PPI", + "PEXT", + }}; + 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 << "REPCA " << 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 << "REPCA " << 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 << label << " mismatch: " << std::setprecision(std::numeric_limits::max_digits10) << actual + << " != " << expected << '\n'; + return false; + } + + bool monitorMatches(const RepcaT& repca, + const std::array& expected, + const char* context) const + { + RealT time = 0.0; + Model::VariableMonitorController monitor(time); + monitor.addMonitor(repca.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,Repca_repca_test_qext,Repca_repca_test_pext," + "Repca_repca_test_vmeas,Repca_repca_test_qmeas," + "Repca_repca_test_pmeas"; + + const auto values = Tokenizer(values_line, ',')(); + if (values.size() != expected.size() + 1) + { + std::cout << "REPCA 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; + } + + /// 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 << 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 index) + { + if (preserved(actual, expected)) + { + return true; + } + std::cout << "REPCA " << what << ' ' << index + << " 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) + { + const size_t row = index(variable); + + if (!variableMatches(static_cast(vector_values[row]), + expected, + what, + variable, + context)) + { + success = false; + } + } + return success; + } + + bool residualsMatch(const RepcaT& repca, + std::initializer_list values, + const char* context = "") const + { + return rowsMatch(repca.getResidual(), values, "residual", context); + } + + template + bool residualsMatch(const RepcaT& repca, + const std::array& values, + const char* context = "") const + { + return rowsMatch(repca.getResidual(), values, "residual", context); + } + + bool stateMatches(const RepcaT& repca, + std::initializer_list values, + const char* context = "") const + { + return rowsMatch(repca.y(), values, "state", context); + } + + template + bool stateMatches(const RepcaT& repca, + const std::array& values, + const char* context = "") const + { + return rowsMatch(repca.y(), values, "state", context); + } + + bool allResidualsWithinInitTolerance(const RepcaT& repca) const + { + bool success = true; + const auto* f = repca.getResidual().getData(); + const auto* yp = repca.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", + RepcaT::INITIALIZATION_TOLERANCE)) + { + success = false; + } + if (!valueUnchanged(yp[row], 0.0, "derivative", row)) + { + 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 << "REPCA " << 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 noteExpectedLogs(const char* message) const + { + const auto previous_verbosity = Log::verbosity(); + Log::setVerbosity(Log::Verbosity::EVERYTHING); + Log::misc() << message << '\n'; + Log::setVerbosity(previous_verbosity); + } + + std::vector expectedJacobian() const + { + return { + {{index(Vars::VMEAS), -6.0}, {index(Vars::VCTRL), 5.0}}, + {{index(Vars::QMEAS), -6.0}, {externalColumn(index(Ext::Q)), 10.0}}, + {{index(Vars::XQPI), -1.0}, + {index(Vars::SFRZ), 1.2}, + {index(Vars::ERQLIM), 1.5}, + {index(Vars::QPI), 0.0}}, + {{index(Vars::XQLAG), -1.4}, {index(Vars::QPI), 0.4}}, + {{index(Vars::PMEAS), -3.5}, {externalColumn(index(Ext::P)), 5.0}}, + {{index(Vars::XPPI), -1.0}, + {index(Vars::EPLIM), 1.8}, + {index(Vars::PPI), 0.0}}, + {{index(Vars::PREF), -3.0}, {index(Vars::PPI), 2.0}}, + {{index(Vars::V), -3.0}, {kBusVrColumn, 1.8}, {kBusViColumn, 0.8}}, + {{index(Vars::VLDC), -2.0}, + {kBusVrColumn, 1.96}, + {kBusViColumn, 0.52}, + {externalColumn(index(Ext::IR)), -0.1096}, + {externalColumn(index(Ext::II)), 0.0968}}, + {{index(Vars::V), 1.0}, {index(Vars::VDROOP), -1.0}, {externalColumn(index(Ext::Q)), 0.8}}, + {{index(Vars::VLDC), 1.0}, + {index(Vars::VDROOP), 0.0}, + {index(Vars::VCTRL), -1.0}}, + {{index(Vars::V), 0.0}, {index(Vars::SFRZ), -1.0}}, + {{index(Vars::VMEAS), -1.0}, + {index(Vars::QMEAS), 0.0}, + {index(Vars::ERQ), -1.0}, + {externalColumn(index(Ext::VREF)), 1.0}, + {externalColumn(index(Ext::QREF)), 0.0}}, + {{index(Vars::ERQ), 1.0}, {index(Vars::ERQDB), -1.0}}, + {{index(Vars::ERQDB), 1.0}, {index(Vars::ERQLIM), -1.0}}, + {{index(Vars::XQPI), 1.0}, {index(Vars::ERQLIM), 2.0}, {index(Vars::QPI), -1.0}}, + {{index(Vars::XQLAG), 2.3}, {index(Vars::QPI), 0.2}, {index(Vars::QEXT), -5.0}}, + {{index(Vars::EF), -1.0}, + {externalColumn(index(Ext::FREQ)), -1.0}, + {externalColumn(index(Ext::FREQREF)), 1.0}}, + {{index(Vars::PMEAS), -1.0}, + {index(Vars::EF), 1.0}, + {index(Vars::EP), -1.0}, + {externalColumn(index(Ext::PREF)), 2.0}}, + {{index(Vars::EP), 0.0}, {index(Vars::EPLIM), -1.0}}, + {{index(Vars::XPPI), 1.0}, {index(Vars::EPLIM), 1.7}, {index(Vars::PPI), -1.0}}, + {{index(Vars::PREF), 1.0}, {index(Vars::PEXT), -2.0}}, + }; + } + + std::vector expectedJacobianAllFlagsOff() const + { + auto expected = expectedJacobian(); + expected[index(Vars::VCTRL)] = { + {index(Vars::VLDC), 0.0}, + {index(Vars::VDROOP), 1.0}, + {index(Vars::VCTRL), -1.0}, + }; + expected[index(Vars::ERQ)] = { + {index(Vars::VMEAS), 0.0}, + {index(Vars::QMEAS), -1.0}, + {index(Vars::ERQ), -1.0}, + {externalColumn(index(Ext::VREF)), 0.0}, + {externalColumn(index(Ext::QREF)), 2.0}, + }; + expected[index(Vars::PEXT)] = { + {index(Vars::PREF), 0.0}, + {index(Vars::PEXT), -2.0}, + }; + return expected; + } + + std::vector expectedJacobianNonunitAlpha() const + { + auto expected = expectedJacobian(); + expected[index(Vars::VMEAS)][index(Vars::VMEAS)] = -7.5; + expected[index(Vars::QMEAS)][index(Vars::QMEAS)] = -7.5; + expected[index(Vars::XQPI)][index(Vars::XQPI)] = -2.5; + expected[index(Vars::XQLAG)][index(Vars::XQLAG)] = -2.9; + expected[index(Vars::PMEAS)][index(Vars::PMEAS)] = -5.0; + expected[index(Vars::XPPI)][index(Vars::XPPI)] = -2.5; + expected[index(Vars::PREF)][index(Vars::PREF)] = -4.5; + return expected; + } + + bool jacobianRowMatches( + const DependencyTracking::Variable::DependencyMap& actual, + const DependencyTracking::Variable::DependencyMap& expected, + size_t row, + const char* source, + RealT tolerance) const + { + if (isEqual(actual, expected, tolerance)) + { + return true; + } + + std::cout << "REPCA " << source << " Jacobian row " << row + << " mismatch\n"; + return false; + } + + bool jacobianMatches( + const std::vector& actual, + const std::vector& expected, + const char* source, + RealT tolerance) const + { + if (actual.size() != expected.size()) + { + std::cout << "REPCA " << source << " Jacobian row-count mismatch\n"; + return false; + } + + bool success = true; + for (size_t row = 0; row < index(Vars::MAXIMUM); ++row) + { + if (!jacobianRowMatches(actual[row], expected[row], row, source, tolerance)) + { + success = false; + } + } + return success; + } + + void numberVariables(Fixture& fixture, + RealT alpha) const + { + auto* y = fixture.repca.y().getData(); + auto* yp = fixture.repca.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.repca.y().setDataUpdated(); + fixture.repca.yp().setDataUpdated(); + fixture.bus.y().setDataUpdated(); + } + + std::vector dependencyTrackingJacobian( + const Data& data, + TestStatus& success, + RealT alpha = 1.0) const + { + using DepVar = DependencyTracking::Variable; + + Fixture fixture(data, kStateVr, kStateVi); + fixture.attachAllInputs(); + setAnswerKeyInputs(fixture); + success *= fixture.prepare(0.0, 0.0); + setAnswerKeyState(fixture.repca); + numberVariables(fixture, alpha); + success *= (fixture.repca.evaluateResidual() == 0); + + std::vector rows(index(Vars::MAXIMUM)); + const auto* f = fixture.repca.getResidual().getData(); + for (size_t row = 0; row < index(Vars::MAXIMUM); ++row) + { + rows[row] = f[row].getDependencies(); + } + return rows; + } + +#ifdef GRIDKIT_ENABLE_ENZYME + std::vector enzymeJacobian( + const Data& data, + TestStatus& success, + RealT alpha = 1.0) const + { + Fixture fixture(data, kStateVr, kStateVi); + fixture.attachAllInputs(); + setAnswerKeyInputs(fixture); + success *= fixture.prepare(0.0, 0.0); + + for (IdxT row = 0; row < fixture.bus.size(); ++row) + { + fixture.bus.setVariableIndex(row, fixture.repca.size() + row); + } + + setAnswerKeyState(fixture.repca); + fixture.repca.updateTime(0.0, alpha); + success *= (fixture.repca.evaluateResidual() == 0); + success *= (fixture.repca.evaluateJacobian() == 0); + success *= (fixture.repca.constructCsr() == 0); + return MapFromCsr(fixture.repca.getCsrJacobian()); + } +#endif + }; + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp b/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp index de98a440f..0feb2a8c7 100644 --- a/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp +++ b/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp @@ -1,4 +1,3 @@ -#include #include #include @@ -252,6 +251,51 @@ namespace GridKit return success.report(__func__); } + TestOutcome repca() + { + using Buses = PhasorDynamics::Controller::RepcaBuses; + using Inputs = PhasorDynamics::Controller::RepcaSignalInputs; + using Vars = PhasorDynamics::Controller::RepcaInternalVariables; + + constexpr IdxT bus_id = static_cast(1); + constexpr IdxT input_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.signal.resize(1); + data.signal[0].signal_id = input_id; + + typename PhasorDynamics::SystemModelData::RepcaDataT repca_data; + repca_data.buses[Buses::bus] = bus_id; + repca_data.signal_inputs[Inputs::ir] = input_id; + repca_data.signal_inputs[Inputs::ii] = input_id; + repca_data.signal_inputs[Inputs::p] = input_id; + repca_data.signal_inputs[Inputs::q] = input_id; + data.repca.push_back(repca_data); + + ScalarT input_value{}; + IdxT input_index = INVALID_INDEX; + + PhasorDynamics::SystemModel system(data); + system.getSignal(input_id)->set(&input_value, &input_index); + + 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/runComponentConnectionTests.cpp b/tests/UnitTests/PhasorDynamics/runComponentConnectionTests.cpp index 127b5107b..b8c7f36f5 100644 --- a/tests/UnitTests/PhasorDynamics/runComponentConnectionTests.cpp +++ b/tests/UnitTests/PhasorDynamics/runComponentConnectionTests.cpp @@ -9,6 +9,7 @@ int main() result += test.genrouEsdc1a(); result += test.genrouHygov(); + result += test.regcaRepca(); return result.summary(); } diff --git a/tests/UnitTests/PhasorDynamics/runControllerRepcaTests.cpp b/tests/UnitTests/PhasorDynamics/runControllerRepcaTests.cpp new file mode 100644 index 000000000..ab10fca94 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/runControllerRepcaTests.cpp @@ -0,0 +1,22 @@ +#include "ControllerRepcaTests.hpp" + +int main() +{ + GridKit::Testing::TestingResults result; + + GridKit::Testing::ControllerRepcaTests test; + + result += test.validation(); + result += test.initializationAndSignals(); + result += test.initializationDomain(); + result += test.residualEquations(); + result += test.reactiveControl(); + result += test.activePowerControl(); + result += test.derivatives(); + result += test.dependencyTracking(); +#ifdef GRIDKIT_ENABLE_ENZYME + result += test.jacobian(); +#endif + + return result.summary(); +} diff --git a/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp b/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp index abb274a0f..65ec715a8 100644 --- a/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp +++ b/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp @@ -16,6 +16,7 @@ int main() result += test.load(); result += test.loadZIP(); result += test.regca(); + result += test.repca(); result += test.genrou(); result += test.genClassical(); result += test.tgov1(); diff --git a/tests/UnitTests/Utilities/CaseFormatTests.hpp b/tests/UnitTests/Utilities/CaseFormatTests.hpp index ce342a321..e4f46bcc4 100644 --- a/tests/UnitTests/Utilities/CaseFormatTests.hpp +++ b/tests/UnitTests/Utilities/CaseFormatTests.hpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include #include @@ -193,6 +194,7 @@ namespace GridKit using BusType = typename BusData::BusType; using Esdc1aData = Exciter::Esdc1aData; using HygovData = Governor::HygovData; + using RepcaData = Controller::RepcaData; const char data[] = R"({ @@ -220,7 +222,18 @@ namespace GridKit { "signal_id": 6, "name": "Under-excitation Limiter"}, { "signal_id": 7, "name": "Hydro Mechanical Power"}, { "signal_id": 8, "name": "Governor Load Reference"}, - { "signal_id": 9, "name": "Governor Auxiliary Power"} + { "signal_id": 9, "name": "Governor Auxiliary Power"}, + { "signal_id": 11, "name": "Branch Current Real"}, + { "signal_id": 12, "name": "Branch Current Imaginary"}, + { "signal_id": 13, "name": "Branch Active Power"}, + { "signal_id": 14, "name": "Branch Reactive Power"}, + { "signal_id": 15, "name": "Frequency"}, + { "signal_id": 16, "name": "Plant Voltage Reference"}, + { "signal_id": 17, "name": "Plant Active Power Reference"}, + { "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"} ], "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} }, @@ -231,6 +244,7 @@ namespace GridKit "Velm": 0.2, "Gmax": 0.98, "Gmin": 0.02, "Tw": 1.2, "At": 1.1, "Dturb": 0.4, "Qnl": 0.08, "Tn": 0.7, "Tnp": 1.4, "db1": 0.01, "db2": 0.02, "Hdam": 1.05, "Gv0": 0.0, "Gv1": 0.2, "Gv2": 0.4, "Gv3": 0.6, "Gv4": 0.8, "Gv5": 1.0, "Pgv0": 0.0, "Pgv1": 0.15, "Pgv2": 0.42, "Pgv3": 0.66, "Pgv4": 0.85, "Pgv5": 1.0}, "mon": ["pmech", "filter", "desiredgate", "gate", "flow", "head"]}, + { "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": "BusFault", "ports": {"bus":1}, "id": "1", "params": {"state0": false, "R":0.0, "X":1e-3} } @@ -256,9 +270,11 @@ namespace GridKit success *= result.gov.size() == 1; success *= result.esdc1a.size() == 1; success *= result.hygov.size() == 1; + success *= result.repca.size() == 1; success *= result.loadz.size() == 0; success *= result.exciter.size() == 1; success *= result.sexspti.size() == 1; + success *= result.signal.size() == 20; success *= result.bus[0].bus_id == 1; success *= result.bus[0].bus_type == BusType::DEFAULT; @@ -418,6 +434,56 @@ namespace GridKit success *= result.hygov[0].monitored_variables.contains(HygovData::MonitorableVariables::flow); success *= result.hygov[0].monitored_variables.contains(HygovData::MonitorableVariables::head); + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::mva]) == 50; + success *= !std::get(result.repca[0].parameters[RepcaData::Parameters::VcompFlag]); + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::RefFlag]); + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::Freqflag]); + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::Tfltr]) == 0.2; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::Vfrz]) == 0.65; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::Rc]) == 0.02; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::Xc]) == 0.03; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::Kc]) == 0.4; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::dbdlow]) == -0.02; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::dbdupper]) == 0.03; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::emax]) == 0.8; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::emin]) == -0.7; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::Kp]) == 2.0; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::Ki]) == 3.0; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::Qmax]) == 0.9; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::Qmin]) == -0.8; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::Tft]) == 0.2; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::Tfv]) == 1.5; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::Tp]) == 0.4; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::fdbd1]) == -0.01; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::fdbd2]) == 0.015; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::Ddn]) == 2.0; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::Dup]) == 1.0; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::femax]) == 0.6; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::femin]) == -0.5; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::Kpg]) == 1.7; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::Kig]) == 1.8; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::Pmax]) == 1.2; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::Pmin]) == 0.1; + success *= std::get(result.repca[0].parameters[RepcaData::Parameters::Tlag]) == 0.5; + success *= result.repca[0].buses[RepcaData::Buses::bus] == 1; + success *= result.repca[0].signal_inputs[RepcaData::SignalInputs::ir] == 11; + success *= result.repca[0].signal_inputs[RepcaData::SignalInputs::ii] == 12; + success *= result.repca[0].signal_inputs[RepcaData::SignalInputs::p] == 13; + success *= result.repca[0].signal_inputs[RepcaData::SignalInputs::q] == 14; + success *= result.repca[0].signal_inputs[RepcaData::SignalInputs::freq] == 15; + success *= result.repca[0].signal_inputs[RepcaData::SignalInputs::vref] == 16; + success *= result.repca[0].signal_inputs[RepcaData::SignalInputs::pref] == 17; + success *= result.repca[0].signal_inputs[RepcaData::SignalInputs::qref] == 18; + success *= result.repca[0].signal_inputs[RepcaData::SignalInputs::freqref] == 19; + success *= result.repca[0].signal_outputs[RepcaData::SignalOutputs::qext] == 20; + success *= result.repca[0].signal_outputs[RepcaData::SignalOutputs::pext] == 21; + success *= result.repca[0].disambiguation_string == "PC1"; + success *= result.repca[0].monitored_variables.contains(RepcaData::MonitorableVariables::qext); + success *= result.repca[0].monitored_variables.contains(RepcaData::MonitorableVariables::pext); + success *= result.repca[0].monitored_variables.contains(RepcaData::MonitorableVariables::vmeas); + success *= result.repca[0].monitored_variables.contains(RepcaData::MonitorableVariables::qmeas); + success *= result.repca[0].monitored_variables.contains(RepcaData::MonitorableVariables::pmeas); + success *= std::get(result.exciter[0].parameters[Exciter::Ieeet1Parameters::Tr]) == 0.0; success *= std::get(result.exciter[0].parameters[Exciter::Ieeet1Parameters::Ka]) == 50.0; success *= std::get(result.exciter[0].parameters[Exciter::Ieeet1Parameters::Ta]) == 0.04;