diff --git a/CHANGELOG.md b/CHANGELOG.md index 066786e09..5f4dffff2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -62,6 +62,7 @@ - Added cmake-format hooks, including in pre-commit. - Added off-nominal tap ratio and phase shift support to the PhasorDynamics `Branch` model. - Added portable Vector class to GridKit +- Added `REECB` converter model implementation for PhasorDynamics. ## v0.1 diff --git a/GridKit/Model/PhasorDynamics/CMakeLists.txt b/GridKit/Model/PhasorDynamics/CMakeLists.txt index 9e99156e4..899a0a07c 100644 --- a/GridKit/Model/PhasorDynamics/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/CMakeLists.txt @@ -57,6 +57,7 @@ add_subdirectory(Branch) add_subdirectory(Bus) add_subdirectory(BusFault) add_subdirectory(BusToSignalAdapter) +add_subdirectory(Converter) add_subdirectory(Exciter) add_subdirectory(Governor) add_subdirectory(Load) diff --git a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp index e2ae075d2..869f8c74e 100644 --- a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp +++ b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include diff --git a/GridKit/Model/PhasorDynamics/Converter/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Converter/CMakeLists.txt new file mode 100644 index 000000000..2002173ee --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/CMakeLists.txt @@ -0,0 +1,6 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +add_subdirectory(REECB) diff --git a/GridKit/Model/PhasorDynamics/Converter/README.md b/GridKit/Model/PhasorDynamics/Converter/README.md index ad38ba19a..ecd71d986 100644 --- a/GridKit/Model/PhasorDynamics/Converter/README.md +++ b/GridKit/Model/PhasorDynamics/Converter/README.md @@ -9,6 +9,6 @@ models and the bus equations, typically through commanded active and reactive cu The GridKit converter documentation includes: -- Renewable Energy Generator/Converter Model REGCA (See [REGCA](REGCA/README.md)) - Renewable Energy Generator/Converter Model REGCB (See [REGCB](REGCB/README.md)) - Renewable Energy Electrical Control Model REECA (See [REECA](REECA/README.md)) +- Renewable Energy Electrical Control Model REECB (See [REECB](REECB/README.md)) diff --git a/GridKit/Model/PhasorDynamics/Converter/REECB/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Converter/REECB/CMakeLists.txt new file mode 100644 index 000000000..03492a503 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REECB/CMakeLists.txt @@ -0,0 +1,54 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +set(_install_headers Reecb.hpp ReecbData.hpp) + +if(GRIDKIT_ENABLE_ENZYME) + gridkit_add_library( + phasor_dynamics_converter_reecb + SOURCES ReecbEnzyme.cpp + HEADERS ${_install_headers} + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal + PRIVATE + ClangEnzymeFlags + COMPILE_OPTIONS + PRIVATE + -mllvm + -enzyme-auto-sparsity=1 + -fno-math-errno) +else() + gridkit_add_library( + phasor_dynamics_converter_reecb + SOURCES Reecb.cpp + HEADERS ${_install_headers} + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal) +endif() + +gridkit_add_library( + phasor_dynamics_converter_reecb_dependency_tracking + SOURCES ReecbDependencyTracking.cpp + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal_dependency_tracking) + +target_link_libraries( + phasor_dynamics_components + INTERFACE GridKit::phasor_dynamics_converter_reecb) +target_link_libraries( + phasor_dynamics_components_dependency_tracking + INTERFACE GridKit::phasor_dynamics_converter_reecb_dependency_tracking) diff --git a/GridKit/Model/PhasorDynamics/Converter/REECB/README.md b/GridKit/Model/PhasorDynamics/Converter/REECB/README.md new file mode 100644 index 000000000..fa8ce3e33 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REECB/README.md @@ -0,0 +1,376 @@ +# **Renewable Energy Electrical Control Model (REECB)** + +REECB is a WECC renewable energy electrical control model for inverter-coupled +resources. In GridKit it is represented as a signal-control model that computes +active- and reactive-current commands. + +Notes: +- Internal electrical quantities and current commands are on model base unless otherwise stated. +- `pe`/`qgen` electrical signal ports are on system base and are converted to + REECB model base through `mva`. +- The `pe` and `qgen` signal ports must be connected together. If both are + omitted, initialization recovers the constant electrical feedback from the + model-base `ipcmd` and `iqcmd` signals. +- `qext`, `pref`, `iqcmd`, and `ipcmd` are model-base command signals. +- REECB uses the REECA voltage-dip freeze structure, but does not include the REECA VDL current-limit curves, active-power speed multiplier, or post-dip hold timers. +- Source governor-response settings may modify active-power ramp-rate or order limits before the model equations are evaluated; the equations below use those effective limits. + +## Block Diagram + +Standard REECB block diagram. + +
+ + + Figure 1: REECB block diagram. Figure courtesy of [PowerWorld](https://www.powerworld.com/WebHelp/) +
+ +## Model Parameters + +Symbol | Units | JSON | Description | Typical Value | Note +------------------------------------|----------|----------|---------------------------------------------------------|---------------|------ +$S^{\mathrm{base}}$ | [MVA] | `mva` | REECB model power base | 0.0 | Block name: `MVABase` +$s_{\mathrm{pf}}$ | [binary] | `PfFlag` | Power-factor control flag | 0 | Block name: `PfFlag`; 1 = power-factor control, 0 = Q control +$s_V$ | [binary] | `VFlag` | Voltage-control mode flag | 0 | Block name: `VFlag`; 1 = Q control, 0 = voltage control +$s_Q$ | [binary] | `QFlag` | Reactive-power control flag | 0 | Block name: `QFlag`; 1 = voltage/Q control, 0 = constant pf or Q control +$s_{PQ}$ | [binary] | `Pqflag` | P/Q priority flag for converter current limit | 0 | Block name: `Pqflag`; 0 = Q priority, 1 = P priority +$T_{\mathrm{rv}}$ | [sec] | `Trv` | Voltage-measurement filter time constant | 0.02 | Block name: `Trv`; if zero, $V_{\mathrm{meas}}$ is algebraic +$T_{\mathrm{p}}$ | [sec] | `Tp` | Electrical-power measurement filter time constant | 0.0 | Block name: `Tp`; if zero, $P_{\mathrm{meas}}$ is algebraic +$V_{\mathrm{ref0}}$ | [p.u.] | `Vref0` | Outer-loop voltage reference | 0.0 | Block name: `Vref0`; initialized to terminal voltage if omitted +$V_{\mathrm{dip}}$ | [p.u.] | `Vdip` | Low-voltage threshold for reactive-current injection logic | 0.85 | Block name: `Vdip` +$V_{\mathrm{up}}$ | [p.u.] | `Vup` | High-voltage threshold for reactive-current injection logic | 1.15 | Block name: `Vup` +$D_{\mathrm{bd1}}$ | [p.u.] | `dbd1` | Overvoltage deadband for voltage-error response | 0.0 | Block name: `dbd1` +$D_{\mathrm{bd2}}$ | [p.u.] | `dbd2` | Undervoltage deadband for voltage-error response | 0.0 | Block name: `dbd2` +$K_{\mathrm{qv}}$ | [p.u.] | `kqv` | Reactive-current injection gain during voltage dip/overvoltage logic | 5.0 | Block name: `kqv` +$I_{\mathrm{qinj}}^{\min}$ | [p.u.] | `Iql1` | Minimum reactive-current injection limit | -1.1 | Block name: `Iql1` +$I_{\mathrm{qinj}}^{\max}$ | [p.u.] | `Iqh1` | Maximum reactive-current injection limit | 1.1 | Block name: `Iqh1` +$Q^{\max}$ | [p.u.] | `Qmax` | Maximum reactive-power control limit | 0.436 | Block name: `Qmax` +$Q^{\min}$ | [p.u.] | `Qmin` | Minimum reactive-power control limit | -0.436 | Block name: `Qmin` +$K_{\mathrm{qp}}$ | [p.u.] | `Kqp` | Reactive-power control proportional gain | 0.0 | Block name: `Kqp` +$K_{\mathrm{qi}}$ | [p.u./s] | `Kqi` | Reactive-power control integral gain | 0.1 | Block name: `Kqi` +$V^{\max}$ | [p.u.] | `Vmax` | Maximum voltage-control limit | 1.1 | Block name: `Vmax` +$V^{\min}$ | [p.u.] | `Vmin` | Minimum voltage-control limit | 0.9 | Block name: `Vmin` +$K_{\mathrm{vp}}$ | [p.u.] | `Kvp` | Voltage-control proportional gain | 18.0 | Block name: `Kvp` +$K_{\mathrm{vi}}$ | [p.u./s] | `Kvi` | Voltage-control integral gain | 5.0 | Block name: `Kvi` +$T_{\mathrm{iq}}$ | [sec] | `Tiq` | Reactive-current command lag time constant | 0.02 | Block name: `Tiq` +$T_{\mathrm{pord}}$ | [sec] | `Tpord` | Active-power order filter time constant | 0.0 | Block name: `Tpord` +$R_P^{\max}$ | [p.u./s] | `dPmax` | Positive active-power order ramp-rate limit | 99.0 | Block name: `dPmax` +$R_P^{\min}$ | [p.u./s] | `dPmin` | Negative active-power order ramp-rate limit | -99.0 | Block name: `dPmin` +$P^{\max}$ | [p.u.] | `Pmax` | Maximum active-power order limit | 1.0 | Block name: `Pmax` +$P^{\min}$ | [p.u.] | `Pmin` | Minimum active-power order limit | 0.0 | Block name: `Pmin` +$I^{\max}$ | [p.u.] | `Imax` | Maximum total converter current | 1.3 | Block name: `Imax` + +### Parameter Validation + +Invalid REECB parameter sets are rejected by the following checks. If source data preprocessing adjusts active-power ramp-rate or order limits, apply these checks to the effective values used by the equations. + +The required checks are: + +```math +\begin{aligned} + &S^{\mathrm{base}} > 0 \\ + &s_{\mathrm{pf}}, s_V, s_Q, s_{PQ} \in \{0,1\} \\ + &T_{\mathrm{rv}}, T_{\mathrm{p}} \ge 0 \\ + &V_{\mathrm{dip}} < V_{\mathrm{up}} \\ + &D_{\mathrm{bd1}} \le 0 \le D_{\mathrm{bd2}} \\ + &I_{\mathrm{qinj}}^{\min} \le I_{\mathrm{qinj}}^{\max} \\ + &Q^{\min} \le Q^{\max} \\ + &V^{\min} \le V^{\max} \\ + &T_{\mathrm{iq}}, T_{\mathrm{pord}} > 0 \\ + &R_P^{\min} < 0 < R_P^{\max} \\ + &P^{\min} \le P^{\max} \\ + &I^{\max} \ge 0 +\end{aligned} +``` + +### Model Derived Parameters + +The off-mode flag complements are: + +```math +\begin{aligned} + s_{\mathrm{pf}}^{\mathrm{off}} &= 1 - s_{\mathrm{pf}} \\ + s_V^{\mathrm{off}} &= 1 - s_V \\ + s_Q^{\mathrm{off}} &= 1 - s_Q \\ + s_{PQ}^{\mathrm{off}} &= 1 - s_{PQ} +\end{aligned} +``` + +System-base electrical signals are converted to REECB model base by: + +```math +\begin{aligned} + x_{\mathrm{model}} &= x_{\mathrm{system}} + \dfrac{S^{\mathrm{sys}}}{S^{\mathrm{base}}} +\end{aligned} +``` + +The same conversion is applied to connected `pe`/`qgen` signals. + +## Model Variables + +### Internal Variables + +#### Differential + +Symbol | Units | Description | Note +------------------------|--------|-------------------------------------|------ +$V_{\mathrm{meas}}$ | [p.u.] | Filtered terminal voltage | State 1 in Fig. 1; source label: `Vmeas`; algebraic when $T_{\mathrm{rv}} = 0$ +$P_{\mathrm{meas}}$ | [p.u.] | Filtered electrical power | State 2 in Fig. 1; source label: `Pmeas`; algebraic when $T_{\mathrm{p}} = 0$ +$x_{\mathrm{PIQ}}$ | [p.u.] | Reactive-power PI controller state | State 3 in Fig. 1; source label: `PIQ` +$x_{\mathrm{PIV}}$ | [p.u.] | Voltage PI controller state | State 4 in Fig. 1; source label: `PIV` +$Q_V$ | [p.u.] | Reactive-current command lag state | State 5 in Fig. 1; source label: `Q_V` +$P_{\mathrm{ord}}$ | [p.u.] | Filtered active-power order | State 6 in Fig. 1; source label: `Pord` + +#### Algebraic + +Symbol | Units | Description | Note +--------------------------------|--------|-------------------------------------|------ +$V_T$ | [p.u.] | Terminal voltage magnitude | +$V_{\mathrm{meas}}^{\mathrm{safe}}$ | [p.u.] | Safe filtered terminal voltage for divider blocks | Lower bounded by 0.01 +$s_{\mathrm{dip}}$ | [binary] | Voltage-dip/overvoltage freeze indicator | 1 when outside voltage thresholds +$V_{\mathrm{err}}$ | [p.u.] | Deadbanded voltage error | Defined by CommonMath `deadband2` +$I_{\mathrm{qv}}$ | [p.u.] | Reactive-current injection candidate | Converter base +$Q_{\mathrm{ref}}$ | [p.u.] | Selected reactive-power reference | From power-factor or external reactive-power command +$e_Q$ | [p.u.] | Reactive-power control error | Limited $Q_{\mathrm{ref}}$ minus $Q_{\mathrm{gen}}$ +$V_{\mathrm{PIQ}}$ | [p.u.] | Reactive-power control PI output | Limited by $V^{\min}$ and $V^{\max}$ +$e_{\mathrm{PIV}}$ | [p.u.] | Voltage-control PI error | Selected voltage-control signal minus $V_{\mathrm{meas}}$ +$f_{\mathrm{pord}}$ | [p.u./s] | Active-power order derivative before ramp-rate limiting | Feeds $r_{\mathrm{pord}}$ +$r_{\mathrm{pord}}$ | [p.u./s] | Ramp-rate-limited active-power order derivative | Feeds $P_{\mathrm{ord}}$ anti-windup +$I_{\mathrm{q}}^{\mathrm{circ}}$ | [p.u.] | Reactive-current limit from converter current circle | Converter base; nonnegative algebraic branch +$I_{\mathrm{p}}^{\mathrm{circ}}$ | [p.u.] | Active-current limit from converter current circle | Converter base; nonnegative algebraic branch +$I_{\mathrm{q}}^{\max}$ | [p.u.] | Final reactive-current upper limit | Converter base; updated by current-limit logic +$I_{\mathrm{p}}^{\max}$ | [p.u.] | Final active-current upper limit | Converter base; updated by current-limit logic +$I_{\mathrm{qbase}}$ | [p.u.] | Base reactive-current command | Converter base; before $s_Q$ selection and reactive-current injection +$I_{\mathrm{q}}^{\mathrm{raw}}$ | [p.u.] | Raw reactive-current command before final limit | Converter base +$I_{\mathrm{q}}^{\mathrm{cmd}}$ | [p.u.] | Reactive-current command output | Converter base +$I_{\mathrm{p}}^{\mathrm{cmd}}$ | [p.u.] | Active-current command output | Converter base + +### External Variables + +#### Differential +None. + +#### Algebraic + +Symbol | Units | Description | Note +-------------------------------------|--------|-----------------------------------|------ +$V_{\mathrm{r}}$ | [p.u.] | Terminal voltage, real component | Owned by bus object +$V_{\mathrm{i}}$ | [p.u.] | Terminal voltage, imaginary component | Owned by bus object +$P_e$ | [p.u.] | Electrical active power | Source label: `Pe`; system-base signal converted to model base +$Q_{\mathrm{gen}}$ | [p.u.] | Reactive-power signal | Source label: `Qgen`; system-base signal converted to model base +$Q_{\mathrm{ext}}$ | [p.u.] | External reactive-power command | Model-base command; optional, defaults to initialized constant; source label: `Qext` +$\phi_{\mathrm{pf}}^{\mathrm{ref}}$ | [rad] | Power-factor angle reference | Source label: `pfaref`; used through tangent block +$P_{\mathrm{ref}}$ | [p.u.] | External active-power reference | Model-base command; optional, defaults to initialized constant + +## Model Equations + +### Differential Equations + +The state-equation residuals use compact limiter notation where applicable. The measurement filters are written in descriptor form; when $T_{\mathrm{rv}} = 0$ or $T_{\mathrm{p}} = 0$, the corresponding residual becomes algebraic. The reactive-current lag requires $T_{\mathrm{iq}} > 0$. + +```math +\begin{aligned} + 0 &= -T_{\mathrm{rv}}\dot V_{\mathrm{meas}} - V_{\mathrm{meas}} + V_T \\ + 0 &= -T_{\mathrm{p}}\dot P_{\mathrm{meas}} - P_{\mathrm{meas}} + P_e \\ + 0 &= + -\dot x_{\mathrm{PIQ}} + + (1 - s_{\mathrm{dip}}) + \text{antiwindup}\!\left( + V_{\mathrm{PIQ}}, + K_{\mathrm{qi}} e_Q, + V^{\min}, + V^{\max} + \right) \\ + 0 &= + -\dot x_{\mathrm{PIV}} + + (1 - s_{\mathrm{dip}}) + \text{antiwindup}\!\left( + I_{\mathrm{qbase}}, + K_{\mathrm{vi}} e_{\mathrm{PIV}}, + -I_{\mathrm{q}}^{\max}, + I_{\mathrm{q}}^{\max} + \right) \\ + 0 &= + -T_{\mathrm{iq}}\dot Q_V + - (1 - s_{\mathrm{dip}})Q_V + + (1 - s_{\mathrm{dip}})Q_{\mathrm{ref}}/V_{\mathrm{meas}}^{\mathrm{safe}} \\ + 0 &= + -\dot P_{\mathrm{ord}} + + (1 - s_{\mathrm{dip}}) + \text{antiwindup}\!\left( + P_{\mathrm{ord}}, + r_{\mathrm{pord}}, + P^{\min}, + P^{\max} + \right) +\end{aligned} +``` + +CommonMath defines the [Anti-Windup](../../../../CommonMath.md#anti-windup-indicator) target and smooth approximation. + +### Algebraic Equations + +The algebraic targets use CommonMath helper notation where applicable: + +```math +\begin{aligned} + 0 &= -V_T^2 + V_\mathrm r^2 + V_\mathrm i^2 \\ + 0 &= -V_\mathrm{meas}^\mathrm{safe} + \max(V_\mathrm{meas}, 0.01) \\ + 0 &= -s_\mathrm{dip} + \text{outside}(V_T, V_\mathrm{dip}, V_\mathrm{up}) \\ + 0 &= -V_\mathrm{err} + \text{deadband2}(V_\mathrm{ref0} - V_\mathrm{meas}, D_\mathrm{bd1}, D_\mathrm{bd2}) \\ + 0 &= -I_\mathrm{qv} + \text{clamp}(K_\mathrm{qv} V_\mathrm{err}, I_\mathrm{qinj}^{\min}, I_\mathrm{qinj}^{\max}) \\ + 0 &= -Q_\mathrm{ref} + + s_\mathrm{pf} P_\mathrm{meas}\tan(\phi_\mathrm{pf}^\mathrm{ref}) + + s_\mathrm{pf}^\mathrm{off} Q_\mathrm{ext} \\ + 0 &= -e_Q + \text{clamp}(Q_\mathrm{ref}, Q^{\min}, Q^{\max}) - Q_\mathrm{gen} \\ + 0 &= -V_\mathrm{PIQ} + \text{clamp}(K_\mathrm{qp} e_Q + x_\mathrm{PIQ}, V^{\min}, V^{\max}) \\ + 0 &= -e_\mathrm{PIV} + s_V V_\mathrm{PIQ} + s_V^\mathrm{off}Q_\mathrm{ref} - V_\mathrm{meas} \\ + 0 &= -T_\mathrm{pord} f_\mathrm{pord} + P_\mathrm{ref} - P_\mathrm{ord} \\ + 0 &= -r_\mathrm{pord} + \text{clamp}(f_\mathrm{pord}, R_P^{\min}, R_P^{\max}) +\end{aligned} +``` + +```math +\begin{aligned} + 0 &= -{I_\mathrm{q}^\mathrm{circ}}^2 + (I^{\max})^2 - s_{PQ}(I_\mathrm{p}^\mathrm{cmd})^2 \\ + 0 &= -{I_\mathrm{p}^\mathrm{circ}}^2 + (I^{\max})^2 - s_{PQ}^\mathrm{off}(I_\mathrm{q}^\mathrm{cmd})^2 \\ + 0 &= -I_\mathrm{q}^{\max} + s_{PQ}^{\mathrm{off}} I^{\max} + s_{PQ} I_\mathrm{q}^{\mathrm{circ}} \\ + 0 &= -I_\mathrm{p}^{\max} + s_{PQ} I^{\max} + s_{PQ}^{\mathrm{off}} I_\mathrm{p}^{\mathrm{circ}} \\ + 0 &= -I_\mathrm{qbase} + \text{clamp}(K_\mathrm{vp} e_\mathrm{PIV} + x_\mathrm{PIV}, -I_\mathrm{q}^{\max}, I_\mathrm{q}^{\max}) \\ + 0 &= -I_\mathrm{q}^\mathrm{raw} + s_Q I_\mathrm{qbase} + s_Q^\mathrm{off} Q_V + s_\mathrm{dip} I_\mathrm{qv} \\ + 0 &= -I_\mathrm{q}^\mathrm{cmd} + \text{clamp}(I_\mathrm{q}^\mathrm{raw}, -I_\mathrm{q}^{\max}, I_\mathrm{q}^{\max}) \\ + 0 &= -I_\mathrm{p}^\mathrm{cmd} + \text{clamp}(P_\mathrm{ord}/V_\mathrm{meas}^\mathrm{safe}, 0, I_\mathrm{p}^{\max}) +\end{aligned} +``` + +The $V_T$, $I_{\mathrm{q}}^{\mathrm{circ}}$, and $I_{\mathrm{p}}^{\mathrm{circ}}$ variables use nonnegative branches of squared algebraic residuals. This matches the REECB current-limit pseudo-code: for $s_{PQ}=0$, $I_{\mathrm{q}}^{\max}=I^{\max}$ and $I_{\mathrm{p}}^{\max}$ is reduced by $I_{\mathrm{q}}^{\mathrm{cmd}}$; for $s_{PQ}=1$, $I_{\mathrm{p}}^{\max}=I^{\max}$ and $I_{\mathrm{q}}^{\max}$ is reduced by $I_{\mathrm{p}}^{\mathrm{cmd}}$. A consistent solution satisfies the nonnegative branch and nonnegative radicands. + +CommonMath defines the helper targets and smooth approximations for [max, clamp, deadband2, and outside](../../../../CommonMath.md#derived-functions). + +## Initialization + +Initialization is performed by evaluating the steady-state residuals in dependency order. Let subscript $0$ denote initial values and set all internal derivatives to zero. REECB recovers the initial electrical operating point from the `ipcmd` and `iqcmd` signal values initialized upstream: + +```math +\begin{aligned} + V_{T,0} &= \sqrt{V_{\mathrm{r},0}^2 + V_{\mathrm{i},0}^2} \\ + V_{\mathrm{meas},0} &= V_{T,0} \\ + V_{\mathrm{meas},0}^{\mathrm{safe}} &= \text{max}(V_{\mathrm{meas},0}, 0.01) \\ + P_{e,0} &= I_{\mathrm{p},0}^{\mathrm{cmd}}V_{\mathrm{meas},0}^{\mathrm{safe}} \\ + Q_{\mathrm{gen},0} &= I_{\mathrm{q},0}^{\mathrm{cmd}}V_{\mathrm{meas},0}^{\mathrm{safe}} \\ + Q_{\mathrm{ext},0} &= Q_{\mathrm{gen},0} \\ + \phi_{\mathrm{pf},0}^{\mathrm{ref}} &= \begin{cases} + \tan^{-1}(Q_{\mathrm{ext},0}/P_{e,0}), & P_{e,0} \neq 0 \\ + 0, & P_{e,0} = 0 + \end{cases} \\ + P_{\mathrm{ref},0} &= \text{clamp}(P_{e,0}, P^{\min}, P^{\max}) +\end{aligned} +``` + +When optional `qext`, `pfaref`, or `pref` signals are connected, initialization writes the recovered steady-state references to those signals. + +If $V_{\mathrm{ref0}}$ is omitted, set $V_{\mathrm{ref0}} = V_{T,0}$. + +```math +\begin{aligned} + P_{\mathrm{meas},0} &= P_{e,0} +\end{aligned} +``` + +Then evaluate the upstream algebraic chain: + +```math +\begin{aligned} + s_{\mathrm{dip},0} &= \text{outside}(V_{T,0}, V_{\mathrm{dip}}, V_{\mathrm{up}}) \\ + V_{\mathrm{err},0} &= \text{deadband2}(V_{\mathrm{ref0}} - V_{\mathrm{meas},0}, D_{\mathrm{bd1}}, D_{\mathrm{bd2}}) \\ + I_{\mathrm{qv},0} &= \text{clamp}(K_{\mathrm{qv}} V_{\mathrm{err},0}, I_{\mathrm{qinj}}^{\min}, I_{\mathrm{qinj}}^{\max}) \\ + Q_{\mathrm{ref},0} &= s_{\mathrm{pf}} P_{\mathrm{meas},0}\tan(\phi_{\mathrm{pf},0}^{\mathrm{ref}}) + s_{\mathrm{pf}}^{\mathrm{off}} Q_{\mathrm{ext},0} \\ + e_{Q,0} &= \text{clamp}(Q_{\mathrm{ref},0}, Q^{\min}, Q^{\max}) - Q_{\mathrm{gen},0} \\ + Q_{V,0} &= \dfrac{Q_{\mathrm{ref},0}}{V_{\mathrm{meas},0}^{\mathrm{safe}}} \\ + P_{\mathrm{ord},0} &= P_{\mathrm{ref},0} +\end{aligned} +``` + +Initialize the reactive-power PI output so its residual and zero-derivative anti-windup condition hold: + +```math +\begin{aligned} + V_{\mathrm{PIQ},0} &= \text{clamp}(K_{\mathrm{qp}} e_{Q,0} + x_{\mathrm{PIQ},0}, V^{\min}, V^{\max}) \\ + e_{\mathrm{PIV},0} &= s_V V_{\mathrm{PIQ},0} + s_V^{\mathrm{off}}Q_{\mathrm{ref},0} - V_{\mathrm{meas},0} +\end{aligned} +``` + +For an unsaturated zero-derivative start, require $e_{Q,0}=0$ for $x_{\mathrm{PIQ}}$ and choose or verify $e_{\mathrm{PIV},0}=0$ for $x_{\mathrm{PIV}}$. Then $x_{\mathrm{PIQ},0}=V_{\mathrm{PIQ},0}-K_{\mathrm{qp}}e_{Q,0}$; when $s_V=1$, set $V_{\mathrm{PIQ},0}=V_{\mathrm{meas},0}$, and when $s_V=0$, the selected $Q_{\mathrm{ref},0}$ signal must equal $V_{\mathrm{meas},0}$. For a saturated PI start, initialize the pre-limit PI output beyond the active limit in the direction that freezes the anti-windup gate. + +Finish by evaluating the current-limit and current-command algebraic residuals in priority order. At the command steps, use the power-flow current targets before final limiting: + +```math +\begin{aligned} + I_{\mathrm{qbase},0}^{\star} &= \dfrac{Q_{\mathrm{gen},0}}{V_{\mathrm{meas},0}^{\mathrm{safe}}} \\ + I_{\mathrm{p},0}^{\star} &= \dfrac{P_{\mathrm{ord},0}}{V_{\mathrm{meas},0}^{\mathrm{safe}}} +\end{aligned} +``` + +These current targets are in the terminal-voltage reference frame. REGCA rotates +the resulting `ipcmd` and `iqcmd` signals into network-frame current injection. + +For $s_{PQ}=0$, use: + +```math +\begin{aligned} +I_{\mathrm{q},0}^{\mathrm{circ}} +\rightarrow I_{\mathrm{q},0}^{\max} +\rightarrow I_{\mathrm{qbase},0} +\rightarrow I_{\mathrm{q},0}^{\mathrm{raw}} +\rightarrow I_{\mathrm{q},0}^{\mathrm{cmd}} +\rightarrow I_{\mathrm{p},0}^{\mathrm{circ}} +\rightarrow I_{\mathrm{p},0}^{\max} +\rightarrow I_{\mathrm{p},0}^{\mathrm{cmd}} +\end{aligned} +``` + +For $s_{PQ}=1$, use: + +```math +\begin{aligned} +I_{\mathrm{p},0}^{\mathrm{circ}} +\rightarrow I_{\mathrm{p},0}^{\max} +\rightarrow I_{\mathrm{p},0}^{\mathrm{cmd}} +\rightarrow I_{\mathrm{q},0}^{\mathrm{circ}} +\rightarrow I_{\mathrm{q},0}^{\max} +\rightarrow I_{\mathrm{qbase},0} +\rightarrow I_{\mathrm{q},0}^{\mathrm{raw}} +\rightarrow I_{\mathrm{q},0}^{\mathrm{cmd}} +\end{aligned} +``` + +After $I_{\mathrm{q},0}^{\max}$ and $I_{\mathrm{qbase},0}$ are known, initialize the voltage PI state from its output residual; the unsaturated zero-derivative start also requires the $e_{\mathrm{PIV},0}=0$ condition above: + +```math +\begin{aligned} + x_{\mathrm{PIV},0} &= I_{\mathrm{qbase},0} - K_{\mathrm{vp}} e_{\mathrm{PIV},0} +\end{aligned} +``` + +The current-circle variables use the nonnegative branch of the squared algebraic residuals, so the radicands must be nonnegative. A standard steady-state initialization assumes $s_{\mathrm{dip},0}=0$. Starts inside voltage-dip or overvoltage logic are not uniquely determined by these closed-form equations. + +## Model Outputs + +Output | Units | Description | Note +----------------|--------|-------------------------------------|------ +`iqcmd` | [p.u.] | Reactive-current command output | Converter base +`ipcmd` | [p.u.] | Active-current command output | Converter base +`vmeas` | [p.u.] | Filtered terminal voltage | +`pmeas` | [p.u.] | Filtered electrical power | +`piq` | [p.u.] | Reactive-power PI controller state | +`piv` | [p.u.] | Voltage PI controller state | +`qv` | [p.u.] | Reactive-current command lag state | +`pord` | [p.u.] | Filtered active-power order | +`qref` | [p.u.] | Selected reactive-power reference | +`sdip` | [binary] | Voltage-dip/overvoltage freeze indicator | +`iqmax` | [p.u.] | Final reactive-current upper limit | Converter base +`ipmax` | [p.u.] | Final active-current upper limit | Converter base +`iqv` | [p.u.] | Reactive-current injection candidate | Converter base +`vqctrl` | [p.u.] | Reactive-power control PI output | +`iqbase` | [p.u.] | Base reactive-current command | Converter base diff --git a/GridKit/Model/PhasorDynamics/Converter/REECB/Reecb.cpp b/GridKit/Model/PhasorDynamics/Converter/REECB/Reecb.cpp new file mode 100644 index 000000000..3b6cacabf --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REECB/Reecb.cpp @@ -0,0 +1,27 @@ +/** + * @file Reecb.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Non-Enzyme instantiation for the REECB electrical-control model. + */ + +#include "ReecbImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + template + int Reecb::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Reecb..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Reecb; + template class Reecb; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REECB/Reecb.hpp b/GridKit/Model/PhasorDynamics/Converter/REECB/Reecb.hpp new file mode 100644 index 000000000..d4ed8b9a7 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REECB/Reecb.hpp @@ -0,0 +1,190 @@ +/** + * @file Reecb.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Declaration of the REECB electrical-control model. + */ + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + template + class BusBase; + + template + class SignalNode; + + namespace Converter + { + /// Internal variables of a `Reecb`. + enum class ReecbInternalVariables : size_t + { + VMEAS, ///< Filtered terminal voltage + PMEAS, ///< Filtered electrical power + XPIQ, ///< Reactive-power PI state + XPIV, ///< Voltage PI state + QV, ///< Reactive-current command lag state + PORD, ///< Filtered active-power order + VT, ///< Terminal voltage magnitude + VMEASSAFE, ///< Safe filtered terminal voltage + SDIP, ///< Voltage-dip/overvoltage freeze indicator + VERR, ///< Deadbanded voltage error + IQV, ///< Reactive-current injection candidate + QREF, ///< Selected reactive-power reference + EQ, ///< Reactive-power control error + VPIQ, ///< Reactive-power PI output + EPIV, ///< Voltage-control PI error + FPORD, ///< Active-power order derivative before ramp limiting + RPORD, ///< Ramp-rate-limited active-power derivative + IQCIRC, ///< Reactive-current limit from current circle + IPCIRC, ///< Active-current limit from current circle + IQMAX, ///< Final reactive-current upper limit + IPMAX, ///< Final active-current upper limit + IQBASE, ///< Base reactive-current command + IQRAW, ///< Raw reactive-current command + IQCMD, ///< Reactive-current command output + IPCMD, ///< Active-current command output + MAXIMUM, + }; + + /// External variables of a `Reecb`. + enum class ReecbExternalVariables : size_t + { + PE, ///< Electrical active-power signal + QGEN, ///< Reactive-power signal + QEXT, ///< External reactive-power command + PFAREF, ///< Power-factor angle reference + PREF, ///< External active-power reference + MAXIMUM, + }; + + template + class Reecb : public Component + { + using Component::alpha_; + using Component::f_; + using Component::gridkit_component_id_; + using Component::J_; + using Component::J_cols_buffer_; + using Component::J_rows_buffer_; + using Component::J_vals_buffer_; + 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 = ReecbData; + using MonitorT = Model::VariableMonitor; + + Reecb(BusT* bus); + Reecb(BusT* bus, const ModelDataT& data); + ~Reecb(); + + int setGridKitComponentID(IdxT) override final; + int allocate() override final; + int verify() const override final; + int initialize() override final; + int tagDifferentiable() override final; + int evaluateResidual() override final; + int evaluateJacobian() override final; + + auto getSignals() + -> ComponentSignals& + { + return signals_; + } + + const Model::VariableMonitorBase* getMonitor() const override; + + __attribute__((always_inline)) inline int evaluateInternalResidual( + ScalarT*, ScalarT*, ScalarT*, ScalarT*, ScalarT*); + + private: + void initModelParams(const ModelDataT& data); + void initializeMonitor(); + void setDerivedParameters(); + + ScalarT toComponentBase(ScalarT value) const; + + ScalarT& Vr(); + ScalarT& Vi(); + + BusT* bus_{nullptr}; + + RealT mva_base_{0}; + RealT PfFlag_{0}; + RealT VFlag_{0}; + RealT QFlag_{0}; + RealT Pqflag_{0}; + RealT Trv_{0}; + RealT Tp_{0}; + RealT Vref0_{0}; + RealT Vdip_{0}; + RealT Vup_{0}; + RealT dbd1_{0}; + RealT dbd2_{0}; + RealT kqv_{0}; + RealT Iql1_{0}; + RealT Iqh1_{0}; + RealT Qmax_{0}; + RealT Qmin_{0}; + RealT Kqp_{0}; + RealT Kqi_{0}; + RealT Vmax_{0}; + RealT Vmin_{0}; + RealT Kvp_{0}; + RealT Kvi_{0}; + RealT Tiq_{0}; + RealT Tpord_{0}; + RealT dPmax_{0}; + RealT dPmin_{0}; + RealT Pmax_{0}; + RealT Pmin_{0}; + RealT Imax_{0}; + RealT va_converter_base_{0}; + RealT pf_off_{1}; + RealT v_off_{1}; + RealT q_off_{1}; + + bool Vref0_given_{false}; + IdxT parameter_error_count_{0}; + + ScalarT qext_set_{0}; + ScalarT pe_set_{0}; + ScalarT qgen_set_{0}; + ScalarT pfaref_set_{0}; + ScalarT pref_set_{0}; + + ComponentSignals signals_; + std::unique_ptr monitor_; + + std::vector ws_; + std::vector ws_indices_; + }; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbData.hpp b/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbData.hpp new file mode 100644 index 000000000..2a0601373 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbData.hpp @@ -0,0 +1,100 @@ +/** + * @file ReecbData.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Modeling data for the REECB electrical-control model. + */ + +#pragma once + +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + /// Parameter keys for the REECB electrical-control model. + enum class ReecbParameters + { + mva, ///< Model MVA base + PfFlag, ///< Power-factor control flag + VFlag, ///< Voltage-control mode flag + QFlag, ///< Reactive-power control flag + Pqflag, ///< P/Q current-priority flag + Trv, ///< Voltage-measurement filter time constant + Tp, ///< Electrical-power measurement filter time constant + Vref0, ///< Outer-loop voltage reference + Vdip, ///< Low-voltage dip threshold + Vup, ///< High-voltage threshold + dbd1, ///< Overvoltage deadband threshold + dbd2, ///< Undervoltage deadband threshold + kqv, ///< Reactive-current injection gain + Iql1, ///< Minimum reactive-current injection limit + Iqh1, ///< Maximum reactive-current injection limit + Qmax, ///< Maximum reactive-power control limit + Qmin, ///< Minimum reactive-power control limit + Kqp, ///< Reactive-power proportional gain + Kqi, ///< Reactive-power integral gain + Vmax, ///< Maximum voltage-control limit + Vmin, ///< Minimum voltage-control limit + Kvp, ///< Voltage-control proportional gain + Kvi, ///< Voltage-control integral gain + Tiq, ///< Reactive-current command lag time constant + Tpord, ///< Active-power order filter time constant + dPmax, ///< Positive active-power ramp-rate limit + dPmin, ///< Negative active-power ramp-rate limit + Pmax, ///< Maximum active-power order limit + Pmin, ///< Minimum active-power order limit + Imax ///< Maximum converter current + }; + + /// Ports for the REECB electrical-control model. + enum class ReecbPorts + { + bus, ///< Terminal bus ID + pe, ///< Electrical active-power signal ID + qgen, ///< Reactive-power signal ID + qext, ///< Optional reactive-power command signal ID + pfaref, ///< Optional power-factor angle reference signal ID + pref, ///< Optional active-power reference signal ID + iqcmd, ///< Reactive-current command output signal ID + ipcmd ///< Active-current command output signal ID + }; + + /// Variables available through the monitor interface. + enum class ReecbMonitorableVariables + { + iqcmd, ///< Reactive-current command output + ipcmd, ///< Active-current command output + vmeas, ///< Filtered terminal voltage + pmeas, ///< Filtered electrical power + piq, ///< Reactive-power PI controller state + piv, ///< Voltage PI controller state + qv, ///< Reactive-current command lag state + pord, ///< Filtered active-power order + qref, ///< Selected reactive-power reference + sdip, ///< Voltage-dip/overvoltage freeze indicator + iqmax, ///< Final reactive-current upper limit + ipmax, ///< Final active-current upper limit + iqv, ///< Reactive-current injection candidate + vqctrl, ///< Reactive-power control PI output + iqbase ///< Base reactive-current command + }; + + template + struct ReecbData : public ComponentData + { + ReecbData() = default; + + using Parameters = ReecbParameters; + using Ports = ReecbPorts; + using MonitorableVariables = ReecbMonitorableVariables; + }; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbDependencyTracking.cpp new file mode 100644 index 000000000..2f81e8274 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbDependencyTracking.cpp @@ -0,0 +1,27 @@ +/** + * @file ReecbDependencyTracking.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Dependency-tracking instantiations for the REECB electrical-control model. + */ + +#include "ReecbImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + template + int Reecb::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Reecb..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Reecb; + template class Reecb; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbEnzyme.cpp b/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbEnzyme.cpp new file mode 100644 index 000000000..0863ce010 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbEnzyme.cpp @@ -0,0 +1,92 @@ +/** + * @file ReecbEnzyme.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Enzyme sparse Jacobian for the REECB electrical-control model. + */ + +#include + +#include "ReecbImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + template + int Reecb::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Reecb..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + + J_.zeroMatrix(); + if (J_rows_buffer_ == nullptr) + { + J_rows_buffer_ = new IdxT[static_cast(size_) * static_cast(size_)]; + J_cols_buffer_ = new IdxT[static_cast(size_) * static_cast(size_)]; + J_vals_buffer_ = new RealT[static_cast(size_) * static_cast(size_)]; + } + + using ModelT = GridKit::PhasorDynamics::Converter::Reecb; + using Fn = GridKit::Enzyme::Sparse::MemberFunctions; + + GridKit::Enzyme::Sparse::DfDy::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + J_); + + GridKit::Enzyme::Sparse::DfDwb::eval(this, + f_.size(), + static_cast(bus_->size()), + (this->getResidualIndices()).data(), + (bus_->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + J_); + + GridKit::Enzyme::Sparse::DfDws::eval(this, + f_.size(), + ws_.size(), + (this->getResidualIndices()).data(), + ws_indices_.data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + J_); + return 0; + } + + template class Reecb; + template class Reecb; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbImpl.hpp b/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbImpl.hpp new file mode 100644 index 000000000..79759da4d --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbImpl.hpp @@ -0,0 +1,724 @@ +/** + * @file ReecbImpl.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Definition of the REECB electrical-control model. + */ + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + using Log = ::GridKit::Utilities::Logger; + + template + Reecb::Reecb(BusT* bus) + : bus_(bus) + { + size_ = static_cast(ReecbInternalVariables::MAXIMUM); + } + + template + Reecb::Reecb(BusT* bus, const ModelDataT& data) + : bus_(bus), + monitor_(std::make_unique(data)) + { + initModelParams(data); + initializeMonitor(); + size_ = static_cast(ReecbInternalVariables::MAXIMUM); + } + + template + Reecb::~Reecb() + { + } + + template + scalar_type& Reecb::Vr() + { + return bus_->Vr(); + } + + template + scalar_type& Reecb::Vi() + { + return bus_->Vi(); + } + + template + void Reecb::setDerivedParameters() + { + va_converter_base_ = mva_base_ * static_cast(1.0e6); + pf_off_ = ONE - PfFlag_; + v_off_ = ONE - VFlag_; + q_off_ = ONE - QFlag_; + } + + template + scalar_type Reecb::toComponentBase( + scalar_type value) const + { + return value * va_system_base_ / va_converter_base_; + } + + template + void Reecb::initModelParams(const ModelDataT& data) + { + using Params = typename ModelDataT::Parameters; + + parameter_error_count_ = 0; + Vref0_given_ = false; + + auto load_required_real = [&](auto key, RealT& target, const char* name) + { + if (!data.parameters.contains(key)) + { + Log::error() << "Reecb: missing required parameter '" << name << "'\n"; + ++parameter_error_count_; + return; + } + + const auto& value = data.parameters.at(key); + if (const auto* real_value = std::get_if(&value)) + { + target = *real_value; + } + else if (const auto* index_value = std::get_if(&value)) + { + target = static_cast(*index_value); + } + else + { + Log::error() << "Reecb: parameter '" << name << "' must be numeric\n"; + ++parameter_error_count_; + } + }; + + auto load_required_switch = [&](auto key, RealT& target, const char* name) + { + if (!data.parameters.contains(key)) + { + Log::error() << "Reecb: missing required parameter '" << name << "'\n"; + ++parameter_error_count_; + return; + } + + const auto& value = data.parameters.at(key); + if (const auto* bool_value = std::get_if(&value)) + { + target = ZERO; + if (*bool_value) + { + target = ONE; + } + } + else if (const auto* index_value = std::get_if(&value); + index_value && (*index_value == 0 || *index_value == 1)) + { + target = static_cast(*index_value); + } + else if (const auto* real_value = std::get_if(&value); + real_value && (*real_value == ZERO || *real_value == ONE) ) + { + target = *real_value; + } + else + { + Log::error() << "Reecb: parameter '" << name << "' must be bool or 0/1\n"; + ++parameter_error_count_; + } + }; + + load_required_real(Params::mva, mva_base_, "mva"); + load_required_switch(Params::PfFlag, PfFlag_, "PfFlag"); + load_required_switch(Params::VFlag, VFlag_, "VFlag"); + load_required_switch(Params::QFlag, QFlag_, "QFlag"); + load_required_switch(Params::Pqflag, Pqflag_, "Pqflag"); + load_required_real(Params::Trv, Trv_, "Trv"); + load_required_real(Params::Tp, Tp_, "Tp"); + if (data.parameters.contains(Params::Vref0)) + { + load_required_real(Params::Vref0, Vref0_, "Vref0"); + Vref0_given_ = true; + } + load_required_real(Params::Vdip, Vdip_, "Vdip"); + load_required_real(Params::Vup, Vup_, "Vup"); + load_required_real(Params::dbd1, dbd1_, "dbd1"); + load_required_real(Params::dbd2, dbd2_, "dbd2"); + load_required_real(Params::kqv, kqv_, "kqv"); + load_required_real(Params::Iql1, Iql1_, "Iql1"); + load_required_real(Params::Iqh1, Iqh1_, "Iqh1"); + load_required_real(Params::Qmax, Qmax_, "Qmax"); + load_required_real(Params::Qmin, Qmin_, "Qmin"); + load_required_real(Params::Kqp, Kqp_, "Kqp"); + load_required_real(Params::Kqi, Kqi_, "Kqi"); + load_required_real(Params::Vmax, Vmax_, "Vmax"); + load_required_real(Params::Vmin, Vmin_, "Vmin"); + load_required_real(Params::Kvp, Kvp_, "Kvp"); + load_required_real(Params::Kvi, Kvi_, "Kvi"); + load_required_real(Params::Tiq, Tiq_, "Tiq"); + load_required_real(Params::Tpord, Tpord_, "Tpord"); + load_required_real(Params::dPmax, dPmax_, "dPmax"); + load_required_real(Params::dPmin, dPmin_, "dPmin"); + load_required_real(Params::Pmax, Pmax_, "Pmax"); + load_required_real(Params::Pmin, Pmin_, "Pmin"); + load_required_real(Params::Imax, Imax_, "Imax"); + setDerivedParameters(); + } + + template + const Model::VariableMonitorBase* Reecb::getMonitor() const + { + return monitor_.get(); + } + + template + void Reecb::initializeMonitor() + { + using Variable = typename ModelDataT::MonitorableVariables; + auto index = [](ReecbInternalVariables variable) + { + return static_cast(variable); + }; + + monitor_->set(Variable::iqcmd, [this, index] + { return y_[index(ReecbInternalVariables::IQCMD)]; }); + monitor_->set(Variable::ipcmd, [this, index] + { return y_[index(ReecbInternalVariables::IPCMD)]; }); + monitor_->set(Variable::vmeas, [this, index] + { return y_[index(ReecbInternalVariables::VMEAS)]; }); + monitor_->set(Variable::pmeas, [this, index] + { return y_[index(ReecbInternalVariables::PMEAS)]; }); + monitor_->set(Variable::piq, [this, index] + { return y_[index(ReecbInternalVariables::XPIQ)]; }); + monitor_->set(Variable::piv, [this, index] + { return y_[index(ReecbInternalVariables::XPIV)]; }); + monitor_->set(Variable::qv, [this, index] + { return y_[index(ReecbInternalVariables::QV)]; }); + monitor_->set(Variable::pord, [this, index] + { return y_[index(ReecbInternalVariables::PORD)]; }); + monitor_->set(Variable::qref, [this, index] + { return y_[index(ReecbInternalVariables::QREF)]; }); + monitor_->set(Variable::sdip, [this, index] + { return y_[index(ReecbInternalVariables::SDIP)]; }); + monitor_->set(Variable::iqmax, [this, index] + { return y_[index(ReecbInternalVariables::IQMAX)]; }); + monitor_->set(Variable::ipmax, [this, index] + { return y_[index(ReecbInternalVariables::IPMAX)]; }); + monitor_->set(Variable::iqv, [this, index] + { return y_[index(ReecbInternalVariables::IQV)]; }); + monitor_->set(Variable::vqctrl, [this, index] + { return y_[index(ReecbInternalVariables::VPIQ)]; }); + monitor_->set(Variable::iqbase, [this, index] + { return y_[index(ReecbInternalVariables::IQBASE)]; }); + } + + template + int Reecb::setGridKitComponentID(IdxT component_id) + { + gridkit_component_id_ = component_id; + return 0; + } + + template + int Reecb::allocate() + { + size_ = static_cast(ReecbInternalVariables::MAXIMUM); + auto size = static_cast(size_); + + f_.assign(size, ScalarT{0}); + y_.assign(size, ScalarT{0}); + yp_.assign(size, ScalarT{0}); + tag_.assign(size, false); + variable_indices_.resize(size); + residual_indices_.resize(size); + + wb_.assign(2, ScalarT{0}); + + auto signal_size = static_cast(ReecbExternalVariables::MAXIMUM); + ws_.assign(signal_size, ScalarT{0}); + ws_indices_.assign(signal_size, INVALID_INDEX); + + for (IdxT j = 0; j < size_; ++j) + { + this->setVariableIndex(j, j); + this->setResidualIndex(j, j); + } + + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y_[static_cast(ReecbInternalVariables::IQCMD)], + &(this->getVariableIndex(static_cast(ReecbInternalVariables::IQCMD)))); + } + + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y_[static_cast(ReecbInternalVariables::IPCMD)], + &(this->getVariableIndex(static_cast(ReecbInternalVariables::IPCMD)))); + } + + return 0; + } + + template + int Reecb::verify() const + { + int ret = static_cast(parameter_error_count_); + + auto check = [&](bool condition, const char* message) + { + if (!condition) + { + Log::error() << "Reecb: " << message << '\n'; + ret += 1; + } + }; + + if (bus_ == nullptr) + { + Log::error() << "Reecb: bus pointer is null\n"; + ret += 1; + } + + check(mva_base_ > ZERO, "mva must be positive"); + check(va_converter_base_ > ZERO, "converter VA base must be positive"); + check(PfFlag_ == ZERO || PfFlag_ == ONE, "PfFlag must be 0 or 1"); + check(VFlag_ == ZERO || VFlag_ == ONE, "VFlag must be 0 or 1"); + check(QFlag_ == ZERO || QFlag_ == ONE, "QFlag must be 0 or 1"); + check(Pqflag_ == ZERO || Pqflag_ == ONE, "Pqflag must be 0 or 1"); + check(Trv_ >= ZERO, "Trv must be non-negative"); + check(Tp_ >= ZERO, "Tp must be non-negative"); + check(Vdip_ < Vup_, "Vdip must be less than Vup"); + check(dbd1_ <= ZERO && ZERO <= dbd2_, "dbd1 <= 0 <= dbd2 is required"); + check(Iql1_ <= Iqh1_, "Iql1 must be less than or equal to Iqh1"); + check(Qmin_ <= Qmax_, "Qmin must be less than or equal to Qmax"); + check(Vmin_ <= Vmax_, "Vmin must be less than or equal to Vmax"); + check(Tiq_ > ZERO, "Tiq must be positive"); + check(Tpord_ > ZERO, "Tpord must be positive"); + check(dPmin_ < ZERO && ZERO < dPmax_, "dPmin < 0 < dPmax is required"); + check(Pmin_ <= Pmax_, "Pmin must be less than or equal to Pmax"); + check(Imax_ >= ZERO, "Imax must be non-negative"); + + const bool has_pe = signals_.template isAttached(); + const bool has_qgen = signals_.template isAttached(); + + if (has_pe != has_qgen) + { + Log::error() << "Reecb: pe and qgen must be connected together\n"; + ret += 1; + } + + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "Reecb: pe signal attached with no linked source\n"; + ret += 1; + } + + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "Reecb: qgen signal attached with no linked source\n"; + ret += 1; + } + + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "Reecb: qext signal attached with no linked source\n"; + ret += 1; + } + + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "Reecb: pfaref signal attached with no linked source\n"; + ret += 1; + } + + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "Reecb: pref signal attached with no linked source\n"; + ret += 1; + } + + return ret; + } + + template + int Reecb::initialize() + { + if (parameter_error_count_ > 0 || verify() > 0) + { + Log::error() << "Reecb: cannot initialize with invalid configuration\n"; + return 1; + } + + const auto VMEAS = static_cast(ReecbInternalVariables::VMEAS); + const auto PMEAS = static_cast(ReecbInternalVariables::PMEAS); + const auto XPIQ = static_cast(ReecbInternalVariables::XPIQ); + const auto XPIV = static_cast(ReecbInternalVariables::XPIV); + const auto QV = static_cast(ReecbInternalVariables::QV); + const auto PORD = static_cast(ReecbInternalVariables::PORD); + const auto VT = static_cast(ReecbInternalVariables::VT); + const auto VMEASSAFE = static_cast(ReecbInternalVariables::VMEASSAFE); + const auto SDIP = static_cast(ReecbInternalVariables::SDIP); + const auto VERR = static_cast(ReecbInternalVariables::VERR); + const auto IQV = static_cast(ReecbInternalVariables::IQV); + const auto QREF = static_cast(ReecbInternalVariables::QREF); + const auto EQ = static_cast(ReecbInternalVariables::EQ); + const auto VPIQ = static_cast(ReecbInternalVariables::VPIQ); + const auto EPIV = static_cast(ReecbInternalVariables::EPIV); + const auto FPORD = static_cast(ReecbInternalVariables::FPORD); + const auto RPORD = static_cast(ReecbInternalVariables::RPORD); + const auto IQCIRC = static_cast(ReecbInternalVariables::IQCIRC); + const auto IPCIRC = static_cast(ReecbInternalVariables::IPCIRC); + const auto IQMAX = static_cast(ReecbInternalVariables::IQMAX); + const auto IPMAX = static_cast(ReecbInternalVariables::IPMAX); + const auto IQBASE = static_cast(ReecbInternalVariables::IQBASE); + const auto IQRAW = static_cast(ReecbInternalVariables::IQRAW); + const auto IQCMD = static_cast(ReecbInternalVariables::IQCMD); + const auto IPCMD = static_cast(ReecbInternalVariables::IPCMD); + + const ScalarT vr = Vr(); + const ScalarT vi = Vi(); + + y_[VT] = std::sqrt(vr * vr + vi * vi); + if (!Vref0_given_) + { + Vref0_ = static_cast(y_[VT]); + } + + y_[VMEAS] = y_[VT]; + y_[VMEASSAFE] = Math::max(y_[VMEAS], static_cast(0.01)); + + const ScalarT ipcmd0 = y_[IPCMD]; + const ScalarT iqcmd0 = y_[IQCMD]; + const ScalarT pe0 = ipcmd0 * y_[VMEASSAFE]; + const ScalarT qgen0 = iqcmd0 * y_[VMEASSAFE]; + + const ScalarT qref0 = qgen0; + const ScalarT pref0 = Math::clamp(pe0, Pmin_, Pmax_); + ScalarT pfaref0{ZERO}; + if (std::abs(static_cast(pe0)) > ZERO) + { + pfaref0 = static_cast(std::atan(static_cast(qref0 / pe0))); + } + + pe_set_ = pe0; + qgen_set_ = qgen0; + qext_set_ = qref0; + pfaref_set_ = pfaref0; + pref_set_ = pref0; + + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(qref0); + } + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(pfaref0); + } + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(pref0); + } + + y_[PMEAS] = pe0; + const ScalarT qref_pf = PfFlag_ * y_[PMEAS] * std::tan(pfaref_set_); + const ScalarT qref_ext = pf_off_ * qext_set_; + + y_[SDIP] = Math::outside(y_[VT], Vdip_, Vup_); + y_[VERR] = Math::deadband2(Vref0_ - y_[VMEAS], dbd1_, dbd2_); + y_[IQV] = Math::clamp(kqv_ * y_[VERR], Iql1_, Iqh1_); + y_[QREF] = qref_pf + qref_ext; + y_[EQ] = Math::clamp(y_[QREF], Qmin_, Qmax_) - qgen_set_; + y_[QV] = y_[QREF] / y_[VMEASSAFE]; + y_[PORD] = pref_set_; + y_[FPORD] = ZERO; + y_[RPORD] = ZERO; + + static constexpr RealT kSmallTol = static_cast(1.0e-9); + static constexpr RealT kSatMargin = static_cast(0.1); + + auto steadyPiOutput = [&](const ScalarT interior_target, + const ScalarT error, + const ScalarT integral_rate, + const ScalarT lower, + const ScalarT upper) -> ScalarT + { + const RealT error_value = static_cast(error); + const RealT rate_value = static_cast(integral_rate); + + if (std::abs(error_value) <= kSmallTol || std::abs(rate_value) <= kSmallTol) + { + return interior_target; + } + + if (rate_value > ZERO) + { + return upper + static_cast(kSatMargin); + } + + return lower - static_cast(kSatMargin); + }; + + const ScalarT xpiq_rate = Kqi_ * y_[EQ]; + const ScalarT vpiq_target = VFlag_ * y_[VMEAS] + v_off_ * y_[QREF]; + const ScalarT vpiq_arg = steadyPiOutput(vpiq_target, + y_[EQ], + xpiq_rate, + static_cast(Vmin_), + static_cast(Vmax_)); + y_[VPIQ] = Math::clamp(vpiq_arg, Vmin_, Vmax_); + y_[EPIV] = VFlag_ * y_[VPIQ] + v_off_ * y_[QREF] - y_[VMEAS]; + y_[XPIQ] = vpiq_arg - Kqp_ * y_[EQ]; + + const ScalarT iqbase_target = qgen0 / y_[VMEASSAFE]; + const ScalarT xpiv_rate = Kvi_ * y_[EPIV]; + const ScalarT ip_star = y_[PORD] / y_[VMEASSAFE]; + + auto initializeReactiveBase = [&]() + { + const ScalarT piv_lower = -y_[IQMAX]; + const ScalarT piv_upper = y_[IQMAX]; + const ScalarT piv_arg = steadyPiOutput(iqbase_target, + y_[EPIV], + xpiv_rate, + piv_lower, + piv_upper); + y_[IQBASE] = Math::clamp(piv_arg, -y_[IQMAX], y_[IQMAX]); + y_[XPIV] = piv_arg - Kvp_ * y_[EPIV]; + }; + + auto initializeReactiveCommand = [&]() + { + y_[IQRAW] = QFlag_ * y_[IQBASE] + q_off_ * y_[QV] + y_[SDIP] * y_[IQV]; + y_[IQCMD] = Math::clamp(y_[IQRAW], -y_[IQMAX], y_[IQMAX]); + }; + + if (Pqflag_ == ZERO) + { + y_[IQCIRC] = Imax_; + y_[IQMAX] = Imax_; + initializeReactiveBase(); + initializeReactiveCommand(); + + const ScalarT ip_radicand = Imax_ * Imax_ - y_[IQCMD] * y_[IQCMD]; + if (static_cast(ip_radicand) < ZERO) + { + Log::error() << "Reecb: initial active-current circle radicand is negative\n"; + return 1; + } + y_[IPCIRC] = std::sqrt(ip_radicand); + y_[IPMAX] = y_[IPCIRC]; + y_[IPCMD] = Math::clamp(ip_star, ZERO, y_[IPMAX]); + } + else + { + y_[IPCIRC] = Imax_; + y_[IPMAX] = Imax_; + y_[IPCMD] = Math::clamp(ip_star, ZERO, y_[IPMAX]); + + const ScalarT iq_radicand = Imax_ * Imax_ - y_[IPCMD] * y_[IPCMD]; + if (static_cast(iq_radicand) < ZERO) + { + Log::error() << "Reecb: initial reactive-current circle radicand is negative\n"; + return 1; + } + y_[IQCIRC] = std::sqrt(iq_radicand); + y_[IQMAX] = y_[IQCIRC]; + initializeReactiveBase(); + initializeReactiveCommand(); + } + + std::fill(yp_.begin(), yp_.end(), ZERO); + + return 0; + } + + template + int Reecb::tagDifferentiable() + { + std::fill(tag_.begin(), tag_.end(), false); + tag_[static_cast(ReecbInternalVariables::VMEAS)] = (Trv_ > ZERO); + tag_[static_cast(ReecbInternalVariables::PMEAS)] = (Tp_ > ZERO); + tag_[static_cast(ReecbInternalVariables::XPIQ)] = true; + tag_[static_cast(ReecbInternalVariables::XPIV)] = true; + tag_[static_cast(ReecbInternalVariables::QV)] = true; + tag_[static_cast(ReecbInternalVariables::PORD)] = true; + return 0; + } + + template + __attribute__((always_inline)) inline int + Reecb::evaluateInternalResidual( + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + ScalarT* ws, + ScalarT* f) + { + const auto VMEAS = static_cast(ReecbInternalVariables::VMEAS); + const auto PMEAS = static_cast(ReecbInternalVariables::PMEAS); + const auto XPIQ = static_cast(ReecbInternalVariables::XPIQ); + const auto XPIV = static_cast(ReecbInternalVariables::XPIV); + const auto QV = static_cast(ReecbInternalVariables::QV); + const auto PORD = static_cast(ReecbInternalVariables::PORD); + const auto VT = static_cast(ReecbInternalVariables::VT); + const auto VMEASSAFE = static_cast(ReecbInternalVariables::VMEASSAFE); + const auto SDIP = static_cast(ReecbInternalVariables::SDIP); + const auto VERR = static_cast(ReecbInternalVariables::VERR); + const auto IQV = static_cast(ReecbInternalVariables::IQV); + const auto QREF = static_cast(ReecbInternalVariables::QREF); + const auto EQ = static_cast(ReecbInternalVariables::EQ); + const auto VPIQ = static_cast(ReecbInternalVariables::VPIQ); + const auto EPIV = static_cast(ReecbInternalVariables::EPIV); + const auto FPORD = static_cast(ReecbInternalVariables::FPORD); + const auto RPORD = static_cast(ReecbInternalVariables::RPORD); + const auto IQCIRC = static_cast(ReecbInternalVariables::IQCIRC); + const auto IPCIRC = static_cast(ReecbInternalVariables::IPCIRC); + const auto IQMAX = static_cast(ReecbInternalVariables::IQMAX); + const auto IPMAX = static_cast(ReecbInternalVariables::IPMAX); + const auto IQBASE = static_cast(ReecbInternalVariables::IQBASE); + const auto IQRAW = static_cast(ReecbInternalVariables::IQRAW); + const auto IQCMD = static_cast(ReecbInternalVariables::IQCMD); + const auto IPCMD = static_cast(ReecbInternalVariables::IPCMD); + + const auto PE = static_cast(ReecbExternalVariables::PE); + const auto QGEN = static_cast(ReecbExternalVariables::QGEN); + const auto QEXT = static_cast(ReecbExternalVariables::QEXT); + const auto PFAREF = static_cast(ReecbExternalVariables::PFAREF); + const auto PREF = static_cast(ReecbExternalVariables::PREF); + + const ScalarT vr = wb[0]; + const ScalarT vi = wb[1]; + + const ScalarT pe = ws[PE]; + const ScalarT qgen = ws[QGEN]; + const ScalarT qext = ws[QEXT]; + const ScalarT pfaref = ws[PFAREF]; + const ScalarT pref = ws[PREF]; + + const ScalarT not_sdip = ONE - y[SDIP]; + const ScalarT xpiq_arg = Kqp_ * y[EQ] + y[XPIQ]; + const ScalarT xpiv_arg = Kvp_ * y[EPIV] + y[XPIV]; + const ScalarT xpiq_limited_rate = Math::antiwindup(xpiq_arg, Kqi_ * y[EQ], Vmin_, Vmax_); + const ScalarT xpiv_limited_rate = + Math::antiwindup(xpiv_arg, Kvi_ * y[EPIV], -y[IQMAX], y[IQMAX]); + const ScalarT pord_limited_rate = Math::antiwindup(y[PORD], y[RPORD], Pmin_, Pmax_); + const ScalarT qref_pf = PfFlag_ * y[PMEAS] * std::tan(pfaref); + const ScalarT qref_ext = pf_off_ * qext; + const ScalarT epiv_target = VFlag_ * y[VPIQ] + v_off_ * y[QREF]; + const ScalarT iqcirc_current_sq = Pqflag_ * y[IPCMD] * y[IPCMD]; + const ScalarT ipcirc_current_sq = (ONE - Pqflag_) * y[IQCMD] * y[IQCMD]; + const ScalarT iqmax_target = (ONE - Pqflag_) * Imax_ + Pqflag_ * y[IQCIRC]; + const ScalarT ipmax_target = Pqflag_ * Imax_ + (ONE - Pqflag_) * y[IPCIRC]; + const ScalarT iqraw_target = QFlag_ * y[IQBASE] + q_off_ * y[QV] + y[SDIP] * y[IQV]; + + f[VMEAS] = -Trv_ * yp[VMEAS] - y[VMEAS] + y[VT]; + f[PMEAS] = -Tp_ * yp[PMEAS] - y[PMEAS] + pe; + f[XPIQ] = -yp[XPIQ] + not_sdip * xpiq_limited_rate; + f[XPIV] = -yp[XPIV] + not_sdip * xpiv_limited_rate; + f[QV] = -Tiq_ * yp[QV] - not_sdip * y[QV] + not_sdip * y[QREF] / y[VMEASSAFE]; + f[PORD] = -yp[PORD] + not_sdip * pord_limited_rate; + + f[VT] = -y[VT] * y[VT] + vr * vr + vi * vi; + f[VMEASSAFE] = -y[VMEASSAFE] + Math::max(y[VMEAS], static_cast(0.01)); + f[SDIP] = -y[SDIP] + Math::outside(y[VT], Vdip_, Vup_); + f[VERR] = -y[VERR] + Math::deadband2(Vref0_ - y[VMEAS], dbd1_, dbd2_); + f[IQV] = -y[IQV] + Math::clamp(kqv_ * y[VERR], Iql1_, Iqh1_); + f[QREF] = -y[QREF] + qref_pf + qref_ext; + f[EQ] = -y[EQ] + Math::clamp(y[QREF], Qmin_, Qmax_) - qgen; + f[VPIQ] = -y[VPIQ] + Math::clamp(xpiq_arg, Vmin_, Vmax_); + f[EPIV] = -y[EPIV] + epiv_target - y[VMEAS]; + f[FPORD] = -Tpord_ * y[FPORD] + pref - y[PORD]; + f[RPORD] = -y[RPORD] + Math::clamp(y[FPORD], dPmin_, dPmax_); + + f[IQCIRC] = -y[IQCIRC] * y[IQCIRC] + Imax_ * Imax_ - iqcirc_current_sq; + f[IPCIRC] = -y[IPCIRC] * y[IPCIRC] + Imax_ * Imax_ - ipcirc_current_sq; + f[IQMAX] = -y[IQMAX] + iqmax_target; + f[IPMAX] = -y[IPMAX] + ipmax_target; + f[IQBASE] = -y[IQBASE] + Math::clamp(xpiv_arg, -y[IQMAX], y[IQMAX]); + f[IQRAW] = -y[IQRAW] + iqraw_target; + f[IQCMD] = -y[IQCMD] + Math::clamp(y[IQRAW], -y[IQMAX], y[IQMAX]); + f[IPCMD] = -y[IPCMD] + Math::clamp(y[PORD] / y[VMEASSAFE], ZERO, y[IPMAX]); + + return 0; + } + + template + int Reecb::evaluateResidual() + { + const auto PE = static_cast(ReecbExternalVariables::PE); + const auto QGEN = static_cast(ReecbExternalVariables::QGEN); + const auto QEXT = static_cast(ReecbExternalVariables::QEXT); + const auto PFAREF = static_cast(ReecbExternalVariables::PFAREF); + const auto PREF = static_cast(ReecbExternalVariables::PREF); + + ws_[PE] = pe_set_; + ws_[QGEN] = qgen_set_; + ws_[QEXT] = qext_set_; + ws_[PFAREF] = pfaref_set_; + ws_[PREF] = pref_set_; + std::fill(ws_indices_.begin(), ws_indices_.end(), INVALID_INDEX); + + if (signals_.template isAttached()) + { + ws_[PE] = toComponentBase( + signals_.template readExternalVariable()); + ws_indices_[PE] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[QGEN] = toComponentBase( + signals_.template readExternalVariable()); + ws_indices_[QGEN] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[QEXT] = signals_.template readExternalVariable(); + ws_indices_[QEXT] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[PFAREF] = signals_.template readExternalVariable(); + ws_indices_[PFAREF] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[PREF] = signals_.template readExternalVariable(); + ws_indices_[PREF] = + signals_.template readExternalVariableIndex(); + } + + wb_[0] = Vr(); + wb_[1] = Vi(); + + evaluateInternalResidual(y_.data(), yp_.data(), wb_.data(), ws_.data(), f_.data()); + return 0; + } + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md index 32df36cdc..a897ec589 100644 --- a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md +++ b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md @@ -144,16 +144,17 @@ are specified: `Genrou` | 6th order machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Tdop`, `Tdopp`, `Tqop`, `Tqopp`, `Xd`, `Xdp`, `Xdpp`, `Xq`, `Xqp`, `Xqpp`, `Xl`, `S10`, `S12`, `mva` | `ir`, `ii`, `p`, `q`, `delta`, `omega`, `speed` `Gensal` | 5th order salient-pole machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Tdop`, `Tdopp`, `Tqopp`, `Xd`, `Xdp`, `Xdpp`, `Xq`, `Xl`, `S10`, `S12`, `mva` | `ir`, `ii`, `p`, `q`, `delta`, `omega`, `speed`, `Eqp`, `psidp`, `psiqpp`, `psidpp`, `vd`, `vq`, `te`, `id`, `iq` `GenClassical` | the classical machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Xdp`, `mva` | `ir`, `ii`, `p`, `q`, `delta`, `omega` + `Reecb` | the REECB renewable electrical-control model | `bus`, `pe`\*, `qgen`\*, `qext`\*, `pfaref`\*, `pref`\*, `iqcmd`, `ipcmd` | `mva`, `PfFlag`, `VFlag`, `QFlag`, `Pqflag`, `Trv`, `Tp`, `Vref0`, `Vdip`, `Vup`, `dbd1`, `dbd2`, `kqv`, `Iql1`, `Iqh1`, `Qmax`, `Qmin`, `Kqp`, `Kqi`, `Vmax`, `Vmin`, `Kvp`, `Kvi`, `Tiq`, `Tpord`, `dPmax`, `dPmin`, `Pmax`, `Pmin`, `Imax` | `iqcmd`, `ipcmd`, `vmeas`, `pmeas`, `piq`, `piv`, `qv`, `pord`, `qref`, `sdip`, `iqmax`, `ipmax`, `iqv`, `vqctrl`, `iqbase` `Tgov1 ` | the TGOV1 governor model | `pmech`, `speed` | `R`, `T1`, `T2`, `T3`, `Pvmax`, `Pvmin`, `Dt` | `none` `Ieeet1` | the IEEET1 exciter model | `bus`, `speed`, `efd`, `vs`\* | `Tr`, `Ka`, `Ta`, `Ke`, `Te`, `Kf`, `Tf`, `Vrmin`, `Vrmax`, `E1`, `E2`, `Se1`, `Se2`, `Ispdlim` | `efd`, `ksat` `SexsPti` | the SEXS-PTI simplified exciter model | `bus`, `efd`, `vs`\* | `Ta`, `Tb`, `Te`, `K`, `Efdmax`, `Efdmin` | `efd` - `Ieeest` | the IEEEST stabilizer model | `input`, `output` | `A1`, `A2`, `A3`, `A4`, `A5`, `A6`, `T1`, `T2`, `T3`, `T4`, `T5`, `T6`, `Ks`, `Lsmin`, `Lsmax`, `Vcl`, `Vcu`, `Tdelay` | `vss` + `Ieeest` | the IEEEST stabilizer model | `input`, `output` | `A1`, `A2`, `A3`, `A4`, `A5`, `A6`, `T1`, `T2`, `T3`, `T4`, `T5`, `T6`, `Ks`, `Lsmin`, `Lsmax`, `Vcl`, `Vcu`, `Tdelay` | `vss` `BusFault` | simple impedance-based fault at a bus | `bus`, `status`\* | `state0`, `R`, `X` | `state`, `ir`, `ii` `BusToSignalAdapter` | signal adapter component for a bus | `bus`, `vr`, `vi`, `ir`, `ii` | | Ports marked with \* are optional and, if missing, will be assumed to be -connected to a constant value. This list is subject to change. - +connected to a constant value. +This list is subject to change. For `Branch`, `tap` and `phase` are optional parameters. If omitted, `tap` defaults to `1.0` and `phase` defaults to `0.0` radians. Bus `bus1` is the tap diff --git a/GridKit/Model/PhasorDynamics/SystemModelData.hpp b/GridKit/Model/PhasorDynamics/SystemModelData.hpp index d4deef66e..06c345916 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelData.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelData.hpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -40,6 +41,7 @@ namespace GridKit using BusDataT = BusData; using BusToSignalAdapterDataT = BusToSignalAdapterData; using BusFaultDataT = BusFaultData; + using ReecbDataT = Converter::ReecbData; using Tgov1DataT = Governor::Tgov1Data; using Ieeet1DataT = Exciter::Ieeet1Data; using SexsPtiDataT = Exciter::SexsPtiData; @@ -93,6 +95,7 @@ namespace GridKit std::vector adapter; ///< bus-to-signal adapters within the model std::vector branch; ///< Branches within the model std::vector bus_fault; ///< Bus faults within the model + std::vector reecb; ///< REECB electrical 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 00a6278c2..102375a5c 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp @@ -136,6 +136,12 @@ namespace GridKit raw_component.get_to(loadzip); sm.loadzip.push_back(loadzip); } + else if (kind == "Reecb") + { + typename SystemModelData::ReecbDataT reecb; + raw_component.get_to(reecb); + sm.reecb.push_back(reecb); + } else if (kind == "Tgov1") { typename SystemModelData::Tgov1DataT gov; diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index 891c2c57d..3f8e57393 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -46,6 +46,7 @@ namespace GridKit using namespace Governor; using namespace Exciter; using namespace Stabilizer; + using namespace Converter; // Set system model tolerances rel_tol_ = 1e-7; @@ -115,6 +116,63 @@ namespace GridKit addComponent(adapter); } + // Add REECB electrical controllers + for (const auto& reecbdata : data.reecb) + { + IdxT bus_index = 0; + if (reecbdata.ports.contains(ReecbPorts::bus)) + { + bus_index = reecbdata.ports.at(ReecbPorts::bus); + } + + auto* reecb = new Reecb(getBus(bus_index), reecbdata); + + if (reecbdata.ports.contains(ReecbPorts::pe)) + { + const IdxT signal = reecbdata.ports.at(ReecbPorts::pe); + constexpr auto PE = ReecbExternalVariables::PE; + reecb->getSignals().template attachSignalNode(getSignal(signal)); + } + if (reecbdata.ports.contains(ReecbPorts::qgen)) + { + const IdxT signal = reecbdata.ports.at(ReecbPorts::qgen); + constexpr auto QGEN = ReecbExternalVariables::QGEN; + reecb->getSignals().template attachSignalNode(getSignal(signal)); + } + if (reecbdata.ports.contains(ReecbPorts::qext)) + { + const IdxT signal = reecbdata.ports.at(ReecbPorts::qext); + constexpr auto QEXT = ReecbExternalVariables::QEXT; + reecb->getSignals().template attachSignalNode(getSignal(signal)); + } + if (reecbdata.ports.contains(ReecbPorts::pfaref)) + { + const IdxT signal = reecbdata.ports.at(ReecbPorts::pfaref); + constexpr auto PFAREF = ReecbExternalVariables::PFAREF; + reecb->getSignals().template attachSignalNode(getSignal(signal)); + } + if (reecbdata.ports.contains(ReecbPorts::pref)) + { + const IdxT signal = reecbdata.ports.at(ReecbPorts::pref); + constexpr auto PREF = ReecbExternalVariables::PREF; + reecb->getSignals().template attachSignalNode(getSignal(signal)); + } + if (reecbdata.ports.contains(ReecbPorts::iqcmd)) + { + const IdxT signal = reecbdata.ports.at(ReecbPorts::iqcmd); + constexpr auto IQCMD = ReecbInternalVariables::IQCMD; + reecb->getSignals().template assignSignalNode(getSignal(signal)); + } + if (reecbdata.ports.contains(ReecbPorts::ipcmd)) + { + const IdxT signal = reecbdata.ports.at(ReecbPorts::ipcmd); + constexpr auto IPCMD = ReecbInternalVariables::IPCMD; + reecb->getSignals().template assignSignalNode(getSignal(signal)); + } + + addComponent(reecb); + } + // Add branches for (const auto& branchdata : data.branch) { diff --git a/docs/Figures/PhasorDynamics_REECB_Diagram.png b/docs/Figures/PhasorDynamics_REECB_Diagram.png new file mode 100644 index 000000000..c3aaa1478 Binary files /dev/null and b/docs/Figures/PhasorDynamics_REECB_Diagram.png differ diff --git a/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp b/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp index 1ce2d8b1f..9c85c6388 100644 --- a/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp +++ b/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp @@ -55,7 +55,31 @@ namespace GridKit return success.report(__func__); } - TestOutcome deadband() + TestOutcome deadband1() + { + TestStatus success = true; + + const ScalarT lower = scalar(-0.05); + const ScalarT upper = scalar(0.10); + const ScalarT tol = scalar(kSmoothTolerance); + + // Type 1 (no-offset) deadband returns ~0 inside the band ... + const ScalarT inside_input = scalar(0.02); + success *= (std::abs(Math::deadband1(inside_input, lower, upper)) < tol * tol); + + // ... and passes the input through unchanged (no offset) outside the band. + const ScalarT far_above_input = scalar(4.0); + const ScalarT far_below_input = scalar(-4.0); + + success *= std::isfinite(Math::deadband1(far_above_input, lower, upper)); + success *= within(Math::deadband1(far_above_input, lower, upper), far_above_input, tol); + success *= std::isfinite(Math::deadband1(far_below_input, lower, upper)); + success *= within(Math::deadband1(far_below_input, lower, upper), far_below_input, tol); + + return success.report(__func__); + } + + TestOutcome deadband2() { TestStatus success = true; @@ -69,12 +93,12 @@ namespace GridKit const ScalarT expected_below = below_input - lower; const ScalarT expected_above = above_input - upper; - success *= within(Math::deadband(below_input, lower, upper), expected_below, tol); - success *= (std::abs(Math::deadband(inside_input, lower, upper)) < tol * tol); - success *= within(Math::deadband(above_input, lower, upper), expected_above, tol); + success *= within(Math::deadband2(below_input, lower, upper), expected_below, tol); + success *= (std::abs(Math::deadband2(inside_input, lower, upper)) < tol * tol); + success *= within(Math::deadband2(above_input, lower, upper), expected_above, tol); - const ScalarT lower_breakpoint = Math::deadband(lower, lower, upper); - const ScalarT upper_breakpoint = Math::deadband(upper, lower, upper); + const ScalarT lower_breakpoint = Math::deadband2(lower, lower, upper); + const ScalarT upper_breakpoint = Math::deadband2(upper, lower, upper); success *= (lower_breakpoint < scalar(0.0)); success *= (upper_breakpoint > scalar(0.0)); @@ -82,7 +106,7 @@ namespace GridKit success *= (std::abs(upper_breakpoint) < tol); const ScalarT x = scalar(-0.4); - success *= (std::abs(Math::deadband(x, lower, upper) + success *= (std::abs(Math::deadband2(x, lower, upper) - (x - Math::clamp(x, lower, upper))) < scalar(kRoundoffTolerance)); @@ -91,17 +115,17 @@ namespace GridKit const ScalarT expected_far_above = far_above_input - upper; const ScalarT expected_far_below = far_below_input - lower; - success *= std::isfinite(Math::deadband(far_above_input, lower, upper)); - success *= within(Math::deadband(far_above_input, lower, upper), expected_far_above, tol); - success *= std::isfinite(Math::deadband(far_below_input, lower, upper)); - success *= within(Math::deadband(far_below_input, lower, upper), expected_far_below, tol); + success *= std::isfinite(Math::deadband2(far_above_input, lower, upper)); + success *= within(Math::deadband2(far_above_input, lower, upper), expected_far_above, tol); + success *= std::isfinite(Math::deadband2(far_below_input, lower, upper)); + success *= within(Math::deadband2(far_below_input, lower, upper), expected_far_below, tol); const ScalarT point = scalar(0.25); const ScalarT above_point = scalar(0.75); const ScalarT below_point = scalar(-0.25); - success *= (std::abs(Math::deadband(above_point, point, point) - (above_point - point)) + success *= (std::abs(Math::deadband2(above_point, point, point) - (above_point - point)) < scalar(kRoundoffTolerance)); - success *= (std::abs(Math::deadband(below_point, point, point) - (below_point - point)) + success *= (std::abs(Math::deadband2(below_point, point, point) - (below_point - point)) < scalar(kRoundoffTolerance)); return success.report(__func__); diff --git a/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp b/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp index 233c36ddd..aecdc7ae8 100644 --- a/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp +++ b/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp @@ -7,7 +7,8 @@ int main() GridKit::Testing::SmoothnessIndicatorTests test; result += test.clamp(); - result += test.deadband(); + result += test.deadband1(); + result += test.deadband2(); result += test.limitIndicators(); result += test.slew(); result += test.linseg(); diff --git a/tests/UnitTests/PhasorDynamics/CMakeLists.txt b/tests/UnitTests/PhasorDynamics/CMakeLists.txt index 932cb13b8..bae05a53f 100644 --- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt +++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt @@ -94,6 +94,14 @@ target_link_libraries( GridKit::phasor_dynamics_components_dependency_tracking GridKit::testing) +add_executable(test_phasor_converter_reecb runConverterReecbTests.cpp) +target_link_libraries( + test_phasor_converter_reecb + GridKit::definitions + GridKit::phasor_dynamics_components + GridKit::phasor_dynamics_components_dependency_tracking + GridKit::testing) + add_executable(test_phasor_stabilizer_ieeest runStabilizerIeeestTests.cpp) target_link_libraries( test_phasor_stabilizer_ieeest @@ -138,6 +146,7 @@ add_test(NAME PhasorDynamicsGovernorTgov1Test COMMAND test_phasor_governor_tgov1 add_test(NAME PhasorDynamicsExciterIeeet1Test COMMAND test_phasor_exciter_ieeet1) add_test(NAME PhasorDynamicsGensalTest COMMAND test_phasor_gensal) add_test(NAME PhasorDynamicsExciterSexsPtiTest COMMAND test_phasor_exciter_sexspti) +add_test(NAME PhasorDynamicsConverterReecbTest COMMAND test_phasor_converter_reecb) add_test(NAME PhasorDynamicsStabilizerIeeestTest COMMAND test_phasor_stabilizer_ieeest) add_test(NAME PhasorDynamicsGenClassicalTest COMMAND test_phasor_gen_classical) add_test(NAME PhasorDynamicsLoadTest COMMAND test_phasor_load) @@ -159,6 +168,7 @@ install( test_phasor_exciter_ieeet1 test_phasor_gensal test_phasor_exciter_sexspti + test_phasor_converter_reecb test_phasor_stabilizer_ieeest test_phasor_gen_classical test_phasor_system diff --git a/tests/UnitTests/PhasorDynamics/ConverterReecbTests.hpp b/tests/UnitTests/PhasorDynamics/ConverterReecbTests.hpp new file mode 100644 index 000000000..2a7a5533a --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/ConverterReecbTests.hpp @@ -0,0 +1,547 @@ +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + template + class ConverterReecbTests + { + public: + using ScalarT = scalar_type; + using IdxT = index_type; + using RealT = typename PhasorDynamics::Component::RealT; + + static constexpr ScalarT kTol = static_cast(1.0e-8); + + TestOutcome constructionAndValidation() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + + PhasorDynamics::Converter::Reecb reecb(&bus, makeReecbData()); + success *= + (reecb.size() + == static_cast(PhasorDynamics::Converter::ReecbInternalVariables::MAXIMUM)); + success *= (reecb.getMonitor() != nullptr); + success *= (reecb.verify() == 0); + + auto wide_band_reecb = makeReecbData(); + wide_band_reecb.parameters[PhasorDynamics::Converter::ReecbParameters::Vdip] = + static_cast(-99.0); + wide_band_reecb.parameters[PhasorDynamics::Converter::ReecbParameters::Vup] = + static_cast(99.0); + PhasorDynamics::Converter::Reecb wide_band_reecb_model( + &bus, + wide_band_reecb); + success *= (wide_band_reecb_model.verify() == 0); + + auto bad_reecb_band = makeReecbData(); + bad_reecb_band.parameters[PhasorDynamics::Converter::ReecbParameters::Vdip] = + static_cast(1.2); + bad_reecb_band.parameters[PhasorDynamics::Converter::ReecbParameters::Vup] = + static_cast(1.2); + PhasorDynamics::Converter::Reecb bad_reecb_band_model( + &bus, + bad_reecb_band); + success *= (bad_reecb_band_model.verify() > 0); + + auto bad_reecb = makeReecbData(); + bad_reecb.parameters[PhasorDynamics::Converter::ReecbParameters::Imax] = + static_cast(-1.0); + PhasorDynamics::Converter::Reecb bad_reecb_model(&bus, bad_reecb); + success *= (bad_reecb_model.verify() > 0); + + return success.report(__func__); + } + + TestOutcome reecbSignalsInitializationAndResidual() + { + using Var = PhasorDynamics::Converter::ReecbInternalVariables; + using Ext = PhasorDynamics::Converter::ReecbExternalVariables; + + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + ScalarT pe_value{0.6}; + ScalarT qgen_value{0.2}; + ScalarT iqcmd_value{0.2}; + ScalarT ipcmd_value{0.6}; + IdxT pe_index = 20; + IdxT qgen_index = 21; + IdxT iqcmd_index = 22; + IdxT ipcmd_index = 23; + + PhasorDynamics::SignalNode pe_node; + PhasorDynamics::SignalNode qgen_node; + PhasorDynamics::SignalNode iqcmd_node; + PhasorDynamics::SignalNode ipcmd_node; + pe_node.set(&pe_value, &pe_index); + qgen_node.set(&qgen_value, &qgen_index); + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + + auto data = makeReecbData(); + data.parameters[PhasorDynamics::Converter::ReecbParameters::QFlag] = static_cast(1); + PhasorDynamics::Converter::Reecb reecb(&bus, data); + reecb.getSignals().template attachSignalNode(&pe_node); + reecb.getSignals().template attachSignalNode(&qgen_node); + reecb.getSignals().template assignSignalNode(&iqcmd_node); + reecb.getSignals().template assignSignalNode(&ipcmd_node); + + success *= (reecb.allocate() == 0); + success *= (reecb.verify() == 0); + iqcmd_node.init(qgen_value); + ipcmd_node.init(pe_value); + success *= (reecb.initialize() == 0); + success *= (reecb.tagDifferentiable() == 0); + success *= (reecb.evaluateResidual() == 0); + + success *= isEqual(reecb.y()[index(Var::VMEAS)], static_cast(1.0), kTol); + success *= isEqual(reecb.y()[index(Var::PMEAS)], pe_value, kTol); + success *= isEqual(reecb.y()[index(Var::QREF)], qgen_value, kTol); + success *= isEqual(iqcmd_node.read(), reecb.y()[index(Var::IQCMD)], kTol); + success *= isEqual(ipcmd_node.read(), reecb.y()[index(Var::IPCMD)], kTol); + success *= (reecb.tag()[index(Var::VMEAS)] == false); + success *= (reecb.tag()[index(Var::PMEAS)] == true); + + for (size_t i = 0; i < reecb.getResidual().size(); ++i) + { + success *= isEqual(reecb.getResidual()[i], static_cast(0.0), kTol); + success *= isEqual(reecb.yp()[i], static_cast(0.0), kTol); + } + + return success.report(__func__); + } + + TestOutcome rejectsHalfConnectedElectricalFeedback() + { + using Ext = PhasorDynamics::Converter::ReecbExternalVariables; + + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + + ScalarT signal_value{0.6}; + IdxT signal_index = 24; + + PhasorDynamics::SignalNode signal_node; + signal_node.set(&signal_value, &signal_index); + + { + PhasorDynamics::Converter::Reecb reecb(&bus, makeReecbData()); + reecb.getSignals().template attachSignalNode(&signal_node); + success *= (reecb.verify() > 0); + } + + { + PhasorDynamics::Converter::Reecb reecb(&bus, makeReecbData()); + reecb.getSignals().template attachSignalNode(&signal_node); + success *= (reecb.verify() > 0); + } + + return success.report(__func__); + } + + TestOutcome reecbCommandSignalInitialization() + { + using Var = PhasorDynamics::Converter::ReecbInternalVariables; + + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + ScalarT iqcmd_value{0.2}; + ScalarT ipcmd_value{0.6}; + IdxT iqcmd_index = 22; + IdxT ipcmd_index = 23; + + PhasorDynamics::SignalNode iqcmd_node; + PhasorDynamics::SignalNode ipcmd_node; + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + + PhasorDynamics::Converter::Reecb reecb(&bus, makeReecbData()); + reecb.getSignals().template assignSignalNode(&iqcmd_node); + reecb.getSignals().template assignSignalNode(&ipcmd_node); + + success *= (reecb.allocate() == 0); + success *= (reecb.verify() == 0); + iqcmd_node.init(static_cast(0.2)); + ipcmd_node.init(static_cast(0.6)); + success *= (reecb.initialize() == 0); + success *= (reecb.tagDifferentiable() == 0); + success *= (reecb.evaluateResidual() == 0); + + success *= isEqual(reecb.y()[index(Var::PMEAS)], static_cast(0.6), kTol); + success *= isEqual(reecb.y()[index(Var::QREF)], static_cast(0.2), kTol); + success *= isEqual(reecb.y()[index(Var::PORD)], static_cast(0.6), kTol); + success *= isEqual(reecb.y()[index(Var::IPCMD)], static_cast(0.6), kTol); + success *= isEqual(reecb.y()[index(Var::IQCMD)], static_cast(0.2), kTol); + + for (size_t i = 0; i < reecb.getResidual().size(); ++i) + { + success *= isEqual(reecb.getResidual()[i], static_cast(0.0), kTol); + success *= isEqual(reecb.yp()[i], static_cast(0.0), kTol); + } + + return success.report(__func__); + } + + TestOutcome reecbElectricalFeedbackUsesMvaBase() + { + using Var = PhasorDynamics::Converter::ReecbInternalVariables; + using Params = PhasorDynamics::Converter::ReecbParameters; + using Ext = PhasorDynamics::Converter::ReecbExternalVariables; + + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + auto data = makeReecbData(); + data.parameters[Params::mva] = static_cast(50.0); + + ScalarT pe_value{0.25}; + ScalarT qgen_value{0.05}; + ScalarT iqcmd_value{0.1}; + ScalarT ipcmd_value{0.5}; + IdxT pe_index = 25; + IdxT qgen_index = 26; + IdxT iqcmd_index = 27; + IdxT ipcmd_index = 28; + + PhasorDynamics::SignalNode pe_node; + PhasorDynamics::SignalNode qgen_node; + PhasorDynamics::SignalNode iqcmd_node; + PhasorDynamics::SignalNode ipcmd_node; + pe_node.set(&pe_value, &pe_index); + qgen_node.set(&qgen_value, &qgen_index); + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + + PhasorDynamics::Converter::Reecb reecb(&bus, data); + reecb.getSignals().template attachSignalNode(&pe_node); + reecb.getSignals().template attachSignalNode(&qgen_node); + reecb.getSignals().template assignSignalNode(&iqcmd_node); + reecb.getSignals().template assignSignalNode(&ipcmd_node); + + success *= (reecb.allocate() == 0); + success *= (reecb.verify() == 0); + iqcmd_node.init(static_cast(0.1)); + ipcmd_node.init(static_cast(0.5)); + success *= (reecb.initialize() == 0); + success *= (reecb.evaluateResidual() == 0); + + success *= isEqual(reecb.y()[index(Var::PMEAS)], static_cast(0.5), kTol); + success *= isEqual(reecb.y()[index(Var::QREF)], static_cast(0.1), kTol); + success *= isEqual(reecb.y()[index(Var::PORD)], static_cast(0.5), kTol); + success *= isEqual(reecb.y()[index(Var::IPCMD)], static_cast(0.5), kTol); + success *= isEqual(reecb.y()[index(Var::IQCMD)], static_cast(0.1), kTol); + + return success.report(__func__); + } + + TestOutcome reecbReferenceFallbackAtAngle() + { + using Var = PhasorDynamics::Converter::ReecbInternalVariables; + using Params = PhasorDynamics::Converter::ReecbParameters; + + TestStatus success = true; + + PhasorDynamics::Bus bus(0.8, 0.6); + bus.allocate(); + bus.initialize(); + + auto data = makeReecbData(); + data.parameters[Params::PfFlag] = static_cast(1); + data.parameters[Params::Qmin] = static_cast(-2.0); + data.parameters[Params::Qmax] = static_cast(2.0); + data.parameters[Params::Pmax] = static_cast(2.0); + data.parameters[Params::Imax] = static_cast(2.0); + + ScalarT iqcmd_value{-0.2}; + ScalarT ipcmd_value{1.6}; + IdxT iqcmd_index = 22; + IdxT ipcmd_index = 23; + + PhasorDynamics::SignalNode iqcmd_node; + PhasorDynamics::SignalNode ipcmd_node; + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + + PhasorDynamics::Converter::Reecb reecb(&bus, data); + reecb.getSignals().template assignSignalNode(&iqcmd_node); + reecb.getSignals().template assignSignalNode(&ipcmd_node); + + success *= (reecb.allocate() == 0); + success *= (reecb.verify() == 0); + iqcmd_node.init(static_cast(-0.2)); + ipcmd_node.init(static_cast(1.6)); + success *= (reecb.initialize() == 0); + success *= (reecb.tagDifferentiable() == 0); + success *= (reecb.evaluateResidual() == 0); + + success *= isEqual(reecb.y()[index(Var::PORD)], static_cast(1.6), kTol); + success *= isEqual(reecb.y()[index(Var::QREF)], static_cast(-0.2), kTol); + success *= isEqual(reecb.y()[index(Var::IPCMD)], static_cast(1.6), kTol); + success *= isEqual(reecb.y()[index(Var::IQCMD)], static_cast(-0.2), kTol); + + for (size_t i = 0; i < reecb.getResidual().size(); ++i) + { + success *= isEqual(reecb.getResidual()[i], static_cast(0.0), kTol); + success *= isEqual(reecb.yp()[i], static_cast(0.0), kTol); + } + + return success.report(__func__); + } + + TestOutcome initializesSaturatedStartConsistently() + { + using Var = PhasorDynamics::Converter::ReecbInternalVariables; + using Params = PhasorDynamics::Converter::ReecbParameters; + + TestStatus success = true; + + PhasorDynamics::Bus bus(1.13, 0.0); + bus.allocate(); + bus.initialize(); + + auto data = makeReecbData(); + data.parameters[Params::QFlag] = static_cast(0); + data.parameters[Params::Pqflag] = static_cast(0); + data.parameters[Params::Vmin] = static_cast(0.9); + data.parameters[Params::Vmax] = static_cast(1.05); + data.parameters[Params::Kvp] = static_cast(10.0); + data.parameters[Params::Kvi] = static_cast(60.0); + data.parameters[Params::Vup] = static_cast(99.0); + data.parameters[Params::Vdip] = static_cast(-99.0); + data.parameters[Params::Imax] = static_cast(1.1); + + ScalarT iqcmd_value{0.15 / 1.13}; + ScalarT ipcmd_value{0.5 / 1.13}; + IdxT iqcmd_index = 22; + IdxT ipcmd_index = 23; + + PhasorDynamics::SignalNode iqcmd_node; + PhasorDynamics::SignalNode ipcmd_node; + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + + PhasorDynamics::Converter::Reecb reecb(&bus, data); + reecb.getSignals().template assignSignalNode(&iqcmd_node); + reecb.getSignals().template assignSignalNode(&ipcmd_node); + + success *= (reecb.allocate() == 0); + success *= (reecb.verify() == 0); + iqcmd_node.init(static_cast(0.15 / 1.13)); + ipcmd_node.init(static_cast(0.5 / 1.13)); + success *= (reecb.initialize() == 0); + + const ScalarT piv_arg = static_cast(10.0) * reecb.y()[index(Var::EPIV)] + + reecb.y()[index(Var::XPIV)]; + success *= (piv_arg < -reecb.y()[index(Var::IQMAX)]); + + success *= (reecb.evaluateResidual() == 0); + for (size_t i = 0; i < reecb.getResidual().size(); ++i) + { + success *= isEqual(reecb.getResidual()[i], static_cast(0.0), kTol); + } + + return success.report(__func__); + } + +#ifdef GRIDKIT_ENABLE_ENZYME + TestOutcome reecbSaturatedVoltagePiJacobianFinite() + { + using Var = PhasorDynamics::Converter::ReecbInternalVariables; + using Params = PhasorDynamics::Converter::ReecbParameters; + + TestStatus success = true; + + PhasorDynamics::Bus bus(1.13, 0.0); + bus.allocate(); + bus.initialize(); + + auto data = makeReecbData(); + data.parameters[Params::QFlag] = static_cast(0); + data.parameters[Params::Pqflag] = static_cast(0); + data.parameters[Params::Vmin] = static_cast(0.9); + data.parameters[Params::Vmax] = static_cast(1.05); + data.parameters[Params::Kvp] = static_cast(10.0); + data.parameters[Params::Kvi] = static_cast(60.0); + data.parameters[Params::Vup] = static_cast(99.0); + data.parameters[Params::Vdip] = static_cast(-99.0); + data.parameters[Params::Imax] = static_cast(1.1); + + ScalarT iqcmd_value{0.15 / 1.13}; + ScalarT ipcmd_value{0.5 / 1.13}; + IdxT iqcmd_index = 22; + IdxT ipcmd_index = 23; + + PhasorDynamics::SignalNode iqcmd_node; + PhasorDynamics::SignalNode ipcmd_node; + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + + PhasorDynamics::Converter::Reecb reecb(&bus, data); + reecb.getSignals().template assignSignalNode(&iqcmd_node); + reecb.getSignals().template assignSignalNode(&ipcmd_node); + + success *= (reecb.allocate() == 0); + success *= (reecb.verify() == 0); + iqcmd_node.init(static_cast(0.15 / 1.13)); + ipcmd_node.init(static_cast(0.5 / 1.13)); + success *= (reecb.initialize() == 0); + success *= (reecb.tagDifferentiable() == 0); + success *= (reecb.evaluateResidual() == 0); + success *= (reecb.evaluateJacobian() == 0); + + const ScalarT piv_arg = static_cast(10.0) * reecb.y()[index(Var::EPIV)] + + reecb.y()[index(Var::XPIV)]; + success *= (piv_arg < -reecb.y()[index(Var::IQMAX)]); + + auto jacobian_entries = reecb.getJacobian().getEntries(false); + const auto& values = std::get<2>(jacobian_entries); + success *= (!values.empty()); + for (const auto value : values) + { + success *= std::isfinite(value); + } + + return success.report(__func__); + } +#endif + + TestOutcome jsonParseAndSystemAssembly() + { + TestStatus success = true; + + std::istringstream input(R"json( +{ + "header": { + "format_version": 0, + "format_revision": 1, + "case_name": "renewable electrical control", + "case_description": "REECB parser test", + "case_comments": "", + "freq_base": 60.0, + "va_base": 100000000.0 + }, + "buses": [ + { "number": 1, "class": "bus", "name": "Bus 1", "init": { "Vr": 1.0, "Vi": 0.0 }, "v_base": 1.0 } + ], + "signals": [ + { "signal_id": 12, "name": "Iqcmd" }, + { "signal_id": 13, "name": "Ipcmd" } + ], + "devices": [ + { + "class": "Reecb", + "ports": { "bus": 1, "iqcmd": 12, "ipcmd": 13 }, + "id": "REE1", + "params": { + "mva": 100.0, "PfFlag": 0, "VFlag": 1, "QFlag": 0, "Pqflag": 1, + "Trv": 0.0, "Tp": 0.02, "Vdip": 0.7, "Vup": 1.2, + "dbd1": -0.01, "dbd2": 0.01, "kqv": 0.0, "Iql1": -1.0, "Iqh1": 1.0, + "Qmax": 1.0, "Qmin": -1.0, "Kqp": 1.0, "Kqi": 0.0, + "Vmax": 1.2, "Vmin": 0.8, "Kvp": 1.0, "Kvi": 0.0, + "Tiq": 0.02, "Tpord": 0.02, "dPmax": 1.0, "dPmin": -1.0, + "Pmax": 1.0, "Pmin": 0.0, "Imax": 2.0 + } + } + ] +} +)json"); + + auto data = PhasorDynamics::parseSystemModelData(input); + success *= (data.reecb.size() == 1); + success *= (std::get( + data.reecb[0].parameters.at(PhasorDynamics::Converter::ReecbParameters::Pqflag)) + == static_cast(1)); + success *= (std::get( + data.reecb[0].parameters.at(PhasorDynamics::Converter::ReecbParameters::mva)) + == static_cast(100.0)); + + PhasorDynamics::SystemModel system(data); + success *= (system.allocate() == 0); + success *= (system.initialize() == 0); + success *= (system.evaluateResidual() == 0); + success *= (system.size() == 27); + + return success.report(__func__); + } + + private: + static size_t index(PhasorDynamics::Converter::ReecbInternalVariables variable) + { + return static_cast(variable); + } + + auto makeReecbData() -> PhasorDynamics::Converter::ReecbData + { + using Params = PhasorDynamics::Converter::ReecbParameters; + using Mon = PhasorDynamics::Converter::ReecbMonitorableVariables; + + PhasorDynamics::Converter::ReecbData data; + data.device_class = "Reecb"; + data.disambiguation_string = "reecb_test"; + data.monitored_variables.insert(Mon::iqcmd); + data.monitored_variables.insert(Mon::ipcmd); + + data.parameters[Params::mva] = static_cast(100.0); + data.parameters[Params::PfFlag] = static_cast(0); + data.parameters[Params::VFlag] = static_cast(1); + data.parameters[Params::QFlag] = static_cast(0); + data.parameters[Params::Pqflag] = static_cast(1); + data.parameters[Params::Trv] = static_cast(0.0); + data.parameters[Params::Tp] = static_cast(0.02); + data.parameters[Params::Vdip] = static_cast(0.7); + data.parameters[Params::Vup] = static_cast(1.2); + data.parameters[Params::dbd1] = static_cast(-0.01); + data.parameters[Params::dbd2] = static_cast(0.01); + data.parameters[Params::kqv] = static_cast(0.0); + data.parameters[Params::Iql1] = static_cast(-1.0); + data.parameters[Params::Iqh1] = static_cast(1.0); + data.parameters[Params::Qmax] = static_cast(1.0); + data.parameters[Params::Qmin] = static_cast(-1.0); + data.parameters[Params::Kqp] = static_cast(1.0); + data.parameters[Params::Kqi] = static_cast(0.0); + data.parameters[Params::Vmax] = static_cast(1.2); + data.parameters[Params::Vmin] = static_cast(0.8); + data.parameters[Params::Kvp] = static_cast(1.0); + data.parameters[Params::Kvi] = static_cast(0.0); + data.parameters[Params::Tiq] = static_cast(0.02); + data.parameters[Params::Tpord] = static_cast(0.02); + data.parameters[Params::dPmax] = static_cast(1.0); + data.parameters[Params::dPmin] = static_cast(-1.0); + data.parameters[Params::Pmax] = static_cast(1.0); + data.parameters[Params::Pmin] = static_cast(0.0); + data.parameters[Params::Imax] = static_cast(2.0); + + return data; + } + }; + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/runConverterReecbTests.cpp b/tests/UnitTests/PhasorDynamics/runConverterReecbTests.cpp new file mode 100644 index 000000000..0645f6133 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/runConverterReecbTests.cpp @@ -0,0 +1,22 @@ +#include "ConverterReecbTests.hpp" + +int main() +{ + GridKit::Testing::TestingResults result; + + GridKit::Testing::ConverterReecbTests test; + + result += test.constructionAndValidation(); + result += test.reecbSignalsInitializationAndResidual(); + result += test.rejectsHalfConnectedElectricalFeedback(); + result += test.reecbCommandSignalInitialization(); + result += test.reecbElectricalFeedbackUsesMvaBase(); + result += test.reecbReferenceFallbackAtAngle(); + result += test.initializesSaturatedStartConsistently(); +#ifdef GRIDKIT_ENABLE_ENZYME + result += test.reecbSaturatedVoltagePiJacobianFinite(); +#endif + result += test.jsonParseAndSystemAssembly(); + + return result.summary(); +}