diff --git a/CHANGELOG.md b/CHANGELOG.md index 0c31a6e41..652845e95 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -61,6 +61,8 @@ - Added `BusToSignalAdapter` component for communicating bus voltages and injection currents. - Added cmake-format hooks, including in pre-commit. - Added off-nominal tap ratio and phase shift support to the PhasorDynamics `Branch` model. +- Added `REGCA` converter model implementation for PhasorDynamics. +- Added `HYGOV` governor model implementation for PhasorDynamics. ## v0.1 diff --git a/GridKit/CommonMath.hpp b/GridKit/CommonMath.hpp index 9bb3a5c3c..daaa11e94 100644 --- a/GridKit/CommonMath.hpp +++ b/GridKit/CommonMath.hpp @@ -153,7 +153,31 @@ namespace GridKit } /** - * @brief Smooth two-sided deadband function + * @brief Smooth Type 1 no-offset two-sided deadband function + * + * Smooth approximation to a deadband that returns zero inside the band and + * passes the input through unchanged outside the band. + * + * @tparam ScalarT - scalar data type + * @tparam RealT - Real data type (see GridKit::ScalarTraits::RealT) + * + * @param[in] x - Input signal + * @param[in] lower - Lower breakpoint + * @param[in] upper - Upper breakpoint + * @return Smooth no-offset deadbanded value + */ + template + __attribute__((always_inline)) inline ScalarT deadband1( + const ScalarT x, + const RealT lower, + const RealT upper) + { + assert(lower <= upper); + return x * (sigmoid(lower - x) + sigmoid(x - upper)); + } + + /** + * @brief Smooth Type 2 offset two-sided deadband function * * Smooth approximation to x - min(max(x, lower), upper), composed from the * smooth ramp function. @@ -164,10 +188,10 @@ namespace GridKit * @param[in] x - Input signal * @param[in] lower - Lower breakpoint * @param[in] upper - Upper breakpoint - * @return Smooth deadbanded value + * @return Smooth offset deadbanded value */ template - __attribute__((always_inline)) inline ScalarT deadband( + __attribute__((always_inline)) inline ScalarT deadband2( const ScalarT x, const RealT lower, const RealT upper) diff --git a/GridKit/CommonMath.md b/GridKit/CommonMath.md index a654487a2..29a9a7901 100644 --- a/GridKit/CommonMath.md +++ b/GridKit/CommonMath.md @@ -49,7 +49,8 @@ The scale $\mu=4\cdot f_{\text{sync}}=240$ is chosen so $\sigma$ behaves like a | `max` | Smooth binary maximum | `REECA` | | `min` | Smooth binary minimum | `REECA` | | `clamp` | Bounded saturation | `IEEEST` | -| `deadband` | Signed two-sided deadband | `REECA` | +| `deadband1` | Type 1 no-offset signed two-sided deadband | - | +| `deadband2` | Type 2 offset signed two-sided deadband | `REECA` | | `slew` | Symmetric slew-rate limiter | - | | `linseg` | Saturated linear segment contribution | `REGCA`, `REECA` | | `above` | Above-lower-limit indicator | - | @@ -84,7 +85,15 @@ x & \ell \le x \le u \\ u & x > u \end{cases} \\ -\text{deadband}(x;\ell,u) +\text{deadband1}(x;\ell,u) +&= +\begin{cases} +x & x < \ell \\ +0 & \ell \le x \le u \\ +x & x > u +\end{cases} +\\ +\text{deadband2}(x;\ell,u) &= \begin{cases} x-\ell & x < \ell \\ @@ -156,7 +165,8 @@ f & x \ge u \land f < 0 \\ \text{max}(x,y) &= y+\rho(x-y) \\ \text{min}(x,y) &= x-\rho(x-y) \\ \text{clamp}(x;\ell,u) &= \ell+\rho(x-\ell)-\rho(x-u) \\ -\text{deadband}(x;\ell,u) &= \rho(x-u)-\rho(\ell-x) \\ +\text{deadband1}(x;\ell,u) &= x\left[\sigma(\ell-x)+\sigma(x-u)\right] \\ +\text{deadband2}(x;\ell,u) &= \rho(x-u)-\rho(\ell-x) \\ \text{slew}(f;r) &= -r+\rho(f+r)-\rho(f-r) \\ \text{linseg}(x;a,b,h) &= \dfrac{h}{b-a}\left[\rho(x-a)-\rho(x-b)\right] \\ \text{above}(x;\ell) &= \sigma(x-\ell) \\ diff --git a/GridKit/Model/PhasorDynamics/CMakeLists.txt b/GridKit/Model/PhasorDynamics/CMakeLists.txt index f76e58c0f..1c7d3eb5b 100644 --- a/GridKit/Model/PhasorDynamics/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/CMakeLists.txt @@ -40,6 +40,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 e48a8a9c5..505a3f572 100644 --- a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp +++ b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp @@ -4,8 +4,10 @@ #include #include #include +#include #include #include +#include #include #include #include diff --git a/GridKit/Model/PhasorDynamics/Converter/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Converter/CMakeLists.txt new file mode 100644 index 000000000..cafb2cc36 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/CMakeLists.txt @@ -0,0 +1,6 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +add_subdirectory(REGCA) diff --git a/GridKit/Model/PhasorDynamics/Converter/REECA/README.md b/GridKit/Model/PhasorDynamics/Converter/REECA/README.md index 6e0d0b83e..834533492 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REECA/README.md +++ b/GridKit/Model/PhasorDynamics/Converter/REECA/README.md @@ -161,7 +161,7 @@ Symbol | Units | Description | $V_T$ | [p.u.] | Terminal voltage magnitude | $V_{\mathrm{meas}}^{\mathrm{safe}}$ | [p.u.] | Safe filtered terminal voltage for divider blocks | Lower bounded by 0.01 $s_{\mathrm{dip}}$ | [binary] | Voltage-dip/overvoltage freeze indicator | 1 when outside voltage thresholds -$V_{\mathrm{err}}$ | [p.u.] | Deadbanded voltage error | Defined by CommonMath `deadband` +$V_{\mathrm{err}}$ | [p.u.] | Deadbanded voltage error | Defined by CommonMath `deadband2` $I_{\mathrm{qv}}$ | [p.u.] | Reactive-current injection candidate | Converter base $Q_{\mathrm{ref}}$ | [p.u.] | Selected reactive-power reference | From power-factor or external reactive-power command $e_Q$ | [p.u.] | Reactive-power control error | Limited $Q_{\mathrm{ref}}$ minus $Q_{\mathrm{gen}}$ @@ -262,7 +262,7 @@ The algebraic targets use CommonMath helper notation where applicable: 0 &= -V_T^2 + V_\mathrm r^2 + V_\mathrm i^2 \\ 0 &= -V_\mathrm{meas}^\mathrm{safe} + \max(V_\mathrm{meas}, 0.01) \\ 0 &= -s_\mathrm{dip} + \text{outside}(V_T, V_\mathrm{dip}, V_\mathrm{up}) \\ - 0 &= -V_\mathrm{err} + \text{deadband}(V_\mathrm{ref0} - V_\mathrm{meas}, D_\mathrm{bd1}, D_\mathrm{bd2}) \\ + 0 &= -V_\mathrm{err} + \text{deadband2}(V_\mathrm{ref0} - V_\mathrm{meas}, D_\mathrm{bd1}, D_\mathrm{bd2}) \\ 0 &= -I_\mathrm{qv} + \text{clamp}(K_\mathrm{qv} V_\mathrm{err}, I_\mathrm{qinj}^{\min}, I_\mathrm{qinj}^{\max}) \\ 0 &= -Q_\mathrm{ref} + s_\mathrm{pf} P_\mathrm{meas}\tan(\phi_\mathrm{pf}^\mathrm{ref}) @@ -290,7 +290,7 @@ The algebraic targets use CommonMath helper notation where applicable: The $V_T$, $I_{\mathrm{q}}^{\mathrm{circ}}$, and $I_{\mathrm{p}}^{\mathrm{circ}}$ variables use nonnegative branches of squared algebraic residuals. This preserves the $s_{PQ}=0$ Q-priority and $s_{PQ}=1$ P-priority current-circle behavior without explicit square roots; a consistent solution should satisfy the nonnegative branch and nonnegative radicands. -CommonMath defines the helper targets and smooth approximations for [min, max, clamp, deadband, and outside](../../../../CommonMath.md#derived-functions). +CommonMath defines the helper targets and smooth approximations for [min, max, clamp, deadband2, and outside](../../../../CommonMath.md#derived-functions). ## Initialization @@ -324,7 +324,7 @@ Then evaluate the upstream algebraic chain: \begin{aligned} V_{\mathrm{meas},0}^{\mathrm{safe}} &= \text{max}(V_{\mathrm{meas},0}, 0.01) \\ s_{\mathrm{dip},0} &= \text{outside}(V_{T,0}, V_{\mathrm{dip}}, V_{\mathrm{up}}) \\ - V_{\mathrm{err},0} &= \text{deadband}(V_{\mathrm{ref0}} - V_{\mathrm{meas},0}, D_{\mathrm{bd1}}, D_{\mathrm{bd2}}) \\ + V_{\mathrm{err},0} &= \text{deadband2}(V_{\mathrm{ref0}} - V_{\mathrm{meas},0}, D_{\mathrm{bd1}}, D_{\mathrm{bd2}}) \\ I_{\mathrm{qv},0} &= \text{clamp}(K_{\mathrm{qv}} V_{\mathrm{err},0}, I_{\mathrm{qinj}}^{\min}, I_{\mathrm{qinj}}^{\max}) \\ Q_{\mathrm{ref},0} &= s_{\mathrm{pf}} P_{\mathrm{meas},0}\tan(\phi_{\mathrm{pf},0}^{\mathrm{ref}}) + s_{\mathrm{pf}}^{\mathrm{off}} Q_{\mathrm{ext},0} \\ e_{Q,0} &= \text{clamp}(Q_{\mathrm{ref},0}, Q^{\min}, Q^{\max}) - Q_{\mathrm{gen},0} \\ diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Converter/REGCA/CMakeLists.txt new file mode 100644 index 000000000..3ffe1fef7 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/CMakeLists.txt @@ -0,0 +1,54 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +set(_install_headers Regca.hpp RegcaData.hpp) + +if(GRIDKIT_ENABLE_ENZYME) + gridkit_add_library( + phasor_dynamics_converter_regca + SOURCES RegcaEnzyme.cpp + HEADERS ${_install_headers} + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal + PRIVATE + ClangEnzymeFlags + COMPILE_OPTIONS + PRIVATE + -mllvm + -enzyme-auto-sparsity=1 + -fno-math-errno) +else() + gridkit_add_library( + phasor_dynamics_converter_regca + SOURCES Regca.cpp + HEADERS ${_install_headers} + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal) +endif() + +gridkit_add_library( + phasor_dynamics_converter_regca_dependency_tracking + SOURCES RegcaDependencyTracking.cpp + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal_dependency_tracking) + +target_link_libraries( + phasor_dynamics_components + INTERFACE GridKit::phasor_dynamics_converter_regca) +target_link_libraries( + phasor_dynamics_components_dependency_tracking + INTERFACE GridKit::phasor_dynamics_converter_regca_dependency_tracking) diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md b/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md index ea352c520..66500962f 100644 --- a/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/README.md @@ -5,9 +5,10 @@ inverter-coupled resources. In GridKit it is represented as a controlled current source at the network interface. Notes: -- LVACM uses the unfiltered terminal voltage $V_T$; LVPL uses the filtered voltage $V_M$. -- Internal currents are on converter base; bus injections are converted to system base in the network interface. -- HVRCM is represented by internal algebraic current $I_{\mathrm{q}}^{\mathrm{extra}}$. +- Internal currents are on converter base. +- Bus injections and exported branch powers are on system base. +- LVACM uses $V_T$; LVPL uses $V_M$. +- HVRCM is represented by algebraic current $I_{\mathrm{q}}^{\mathrm{extra}}$. ## Block Diagram @@ -21,23 +22,26 @@ Standard REGCA converter-interface model. ## Model Parameters -Symbol | Units | Description | Typical Value | Note ----------------------------------|----------|-------------------------------------------------------|---------------|------ -$P_{\mathrm{0}}$ | [p.u.] | Initial active power injection | | On system base -$Q_{\mathrm{0}}$ | [p.u.] | Initial reactive power injection | | On system base -$S^{\mathrm{conv}}$ | [MVA] | Converter/model power base | TBD | -$T_{\mathrm{g}}$ | [sec] | Converter current-control lag time constant | TBD | -$T_M$ | [sec] | Terminal voltage sensor time constant | TBD | Block name: `Tfltr` -$R_{\mathrm{q}}^{\max}$ | [p.u./s] | Reactive-current recovery positive rate limit | TBD | Block name: `Iqrmax` -$R_{\mathrm{q}}^{\min}$ | [p.u./s] | Reactive-current recovery negative rate limit | TBD | Block name: `Iqrmin` -$R_{\mathrm{p}}^{\max}$ | [p.u./s] | Active-current magnitude recovery rate limit | TBD | Block name: `rrpwr` -$s_L$ | [binary] | LVPL switch | TBD | Block name: `LPVLSW` -$I_{L1}$ | [p.u.] | LVPL upper-current ceiling | TBD | Block name: `Lvpl1` -$V_{L0}$ | [p.u.] | LVPL zero-crossing voltage | TBD | Block name: `xerox` -$V_{L1}$ | [p.u.] | LVPL upper breakpoint voltage | TBD | Block name: `brkpt` -$V_{A0}$ | [p.u.] | LVACM lower breakpoint voltage | TBD | Block name: `Lvpnt0` -$V_{A1}$ | [p.u.] | LVACM upper breakpoint voltage | TBD | Block name: `Lvpnt1` -$V_{\mathrm{hv}}^{\max}$ | [p.u.] | Terminal-voltage ceiling for HV reactive management | TBD | Block name: `Vlim` +Symbol | Units | JSON | Description | Typical Value | Note +---------------------------------|----------|----------|-------------------------------------------------------|---------------|------ +$P_0$ | [p.u.] | `P0` | Initial active power injection on system base | 1.0 | Required initialization source +$Q_0$ | [p.u.] | `Q0` | Initial reactive power injection on system base | 0.0 | Required initialization source +$S^{\mathrm{base}}$ | [MVA] | `mva` | REGCA model power base | 100.0 | +$T_{\mathrm{g}}$ | [sec] | `Tg` | Converter current-control lag time constant | 0.02 | Block name: `Tg` +$T_M$ | [sec] | `TM` | Terminal voltage sensor time constant | 0.02 | Block name: `Tfltr` +$R_{\mathrm{q}}^{\max}$ | [p.u./s] | `Rqmax` | Reactive-current recovery positive rate limit | 999.0 | Block name: `Iqrmax`; GridKit requires $R_{\mathrm{q}}^{\max} > 0$ +$R_{\mathrm{q}}^{\min}$ | [p.u./s] | `Rqmin` | Reactive-current recovery negative rate limit | -999.0 | Block name: `Iqrmin`; GridKit requires $R_{\mathrm{q}}^{\min} < 0$ +$R_{\mathrm{p}}^{\max}$ | [p.u./s] | `Rpmax` | Active-current magnitude recovery rate limit | 999.0 | Block name: `rrpwr` +$s_L$ | [binary] | `sL` | LVPL switch | 1 | Block name: `LPVLSW` +$I_{L1}$ | [p.u.] | `IL1` | LVPL upper-current ceiling | 1.1 | Block name: `LVPL1` +$V_{L0}$ | [p.u.] | `VL0` | LVPL zero-crossing voltage | 0.4 | Block name: `zerox` +$V_{L1}$ | [p.u.] | `VL1` | LVPL upper breakpoint voltage | 0.9 | Block name: `brkpt` +$V_{A0}$ | [p.u.] | `VA0` | LVACM lower breakpoint voltage | 0.4 | Block name: `LVPnt0` +$V_{A1}$ | [p.u.] | `VA1` | LVACM upper breakpoint voltage | 0.9 | Block name: `LVPnt1` +$V_{\mathrm{hv}}^{\max}$ | [p.u.] | `Vhvmax` | Terminal-voltage ceiling for HV reactive management | 1.2 | Block name: `VLim` + +PowerWorld fields `Qmin`, `Khv`, and `Xe` are not parameters of this GridKit +REGCA implementation. ### Parameter Validation @@ -45,7 +49,7 @@ Implementations should reject or report invalid parameter sets: ```math \begin{aligned} - S^{\mathrm{conv}} &> 0 & + S^{\mathrm{base}} &> 0 & T_{\mathrm{g}} &> 0 & T_M &> 0 \\ R_{\mathrm{p}}^{\max} &> 0 & @@ -60,16 +64,14 @@ Implementations should reject or report invalid parameter sets: ### Model Derived Parameters -The smooth active-current bound equations use $M_{\mathrm{p}}$, a numerical -relaxation for inactive $\pm\infty$ rate bounds: +The active-current bounds use finite surrogate $M_{\mathrm{p}}$ for inactive +$\pm\infty$ limits: ```math M_{\mathrm{p}} = 100 R_{\mathrm{p}}^{\max} ``` -$M_{\mathrm{p}}$ is not a physical REGCA parameter; it should be large enough -that inactive bounds do not bind expected $f_{\mathrm{p}}$ values while staying -moderate enough to keep the smooth clamp well conditioned. +$M_{\mathrm{p}}$ is not a physical REGCA parameter. ## Model Variables @@ -94,6 +96,8 @@ $I_L$ | [p.u.] | LVPL upper-limit current curve $I_{\mathrm{r}}$ | [p.u.] | Injected current, real component on network reference frame | Converter base $\ell_{\mathrm{p}}$ | [p.u./s] | Smooth active-current lower rate bound | Equivalent to diagram `Rdown` $u_{\mathrm{p}}$ | [p.u./s] | Smooth active-current upper rate bound | Effective `Rup`; includes LVPL anti-windup when $s_L=1$ +$P_{\mathrm{br}}$ | [p.u.] | Active-power output | System base +$Q_{\mathrm{br}}$ | [p.u.] | Reactive-power output | System base ### External Variables @@ -106,129 +110,71 @@ Symbol | Units | Description --------------------------------|--------|------------------------------------------------------------------|------ $V_{\mathrm{r}}$ | [p.u.] | Terminal voltage, real component on network reference frame | Owned by bus object $V_{\mathrm{i}}$ | [p.u.] | Terminal voltage, imaginary component on network reference frame | Owned by bus object -$I_{\mathrm{q}}^{\mathrm{cmd}}$ | [p.u.] | Reactive-current command | Converter base; owned by REEC, constant if no REEC is connected -$I_{\mathrm{p}}^{\mathrm{cmd}}$ | [p.u.] | Active-current command | Converter base; owned by REEC, constant if no REEC is connected +$I_{\mathrm{p}}^{\mathrm{cmd}}$ | [p.u.] | Active-current command in the terminal-voltage reference frame | Converter base; owned by REEC, constant if no REEC is connected +$I_{\mathrm{q}}^{\mathrm{cmd}}$ | [p.u.] | Reactive-current command in the terminal-voltage reference frame | Converter base; owned by REEC, constant if no REEC is connected -## Model Equations - -Define the pre-limit current derivatives: +REGCA initializes its current commands from the required `P0`/`Q0` power-flow +injection and writes those resolved commands to any attached `ipcmd`/`iqcmd` +ports. If no controller is connected, the resolved initialization commands are +held constant during residual evaluation. -```math -\begin{aligned} - f_{\mathrm{q}} &= \dfrac{1}{T_{\mathrm{g}}}(I_{\mathrm{q}}^{\mathrm{cmd}} - I_{\mathrm{q}}) \\ - f_{\mathrm{p}} &= \dfrac{1}{T_{\mathrm{g}}}(I_{\mathrm{p}}^{\mathrm{cmd}} - I_{\mathrm{p}}) -\end{aligned} -``` +## Model Equations ### Differential Equations -The exact state equations are - -```math -\begin{aligned} - 0 &= -T_M \dot V_M - V_M + V_T \\ - 0 &= -\dot I_{\mathrm{q}} + - \begin{cases} - \min(f_{\mathrm{q}}, R_{\mathrm{q}}^{\max}) & Q_{\mathrm{0}} > 0 \\ - \max(f_{\mathrm{q}}, R_{\mathrm{q}}^{\min}) & Q_{\mathrm{0}} \le 0 - \end{cases} \\ - 0 &= -\dot I_{\mathrm{p}} + \text{clamp}(f_{\mathrm{p}}, \ell_{\mathrm{p}}, u_{\mathrm{p}}) -\end{aligned} -``` - -The implemented smooth state equations are +The state equations use CommonMath helper notation. The $I_{\mathrm{q}}$ +limiter branch is selected by the initialized reactive-current command. ```math \begin{aligned} 0 &= -T_M \dot V_M - V_M + V_T \\ - 0 &= -\dot I_{\mathrm{q}} + + 0 &= -T_{\mathrm{g}}\dot I_{\mathrm{q}} + + (I_{\mathrm{q}}^{\mathrm{cmd}} - I_{\mathrm{q}}) + \begin{cases} - f_{\mathrm{q}} - \rho(f_{\mathrm{q}} - R_{\mathrm{q}}^{\max}) - & Q_{\mathrm{0}} > 0 \\ - f_{\mathrm{q}} + \rho(R_{\mathrm{q}}^{\min} - f_{\mathrm{q}}) - & Q_{\mathrm{0}} \le 0 + -\rho(I_{\mathrm{q}}^{\mathrm{cmd}} - I_{\mathrm{q}} - T_{\mathrm{g}}R_{\mathrm{q}}^{\max}) + & I_{\mathrm{q},0}^{\mathrm{cmd}} > 0 \\ + \rho(T_{\mathrm{g}}R_{\mathrm{q}}^{\min} - (I_{\mathrm{q}}^{\mathrm{cmd}} - I_{\mathrm{q}})) + & I_{\mathrm{q},0}^{\mathrm{cmd}} \le 0 \end{cases} \\ - 0 &= -\dot I_{\mathrm{p}} + - \ell_{\mathrm{p}} - + \rho(f_{\mathrm{p}} - \ell_{\mathrm{p}}) - - \rho(f_{\mathrm{p}} - u_{\mathrm{p}}) + 0 &= -T_{\mathrm{g}}\dot I_{\mathrm{p}} + + T_{\mathrm{g}}\ell_{\mathrm{p}} + + \rho(I_{\mathrm{p}}^{\mathrm{cmd}} - I_{\mathrm{p}} - T_{\mathrm{g}}\ell_{\mathrm{p}}) + - \rho(I_{\mathrm{p}}^{\mathrm{cmd}} - I_{\mathrm{p}} - T_{\mathrm{g}}u_{\mathrm{p}}) \end{aligned} ``` -Here $\rho$ is GridKit's smooth ramp function. The $I_{\mathrm{q}}$ branch is -selected by initial reactive power $Q_{\mathrm{0}}$. The $I_{\mathrm{p}}$ -equation is the smooth clamp of $f_{\mathrm{p}}$ between the algebraic bounds -$\ell_{\mathrm{p}}$ and $u_{\mathrm{p}}$. +CommonMath defines the [linear segment, ramp, and sigmoid helpers](../../../../CommonMath.md#derived-functions). ### Algebraic Equations -The piecewise definitions in this section switch on continuous states, unlike -the $I_{\mathrm{q}}$ differential branch selected by initial conditions. The -exact algebraic targets are: - -```math -\begin{aligned} - 0 &= -V_T^2 + V_{\mathrm{r}}^2 + V_{\mathrm{i}}^2 \\ - I_{\mathrm{i}} &= -I_{\mathrm{q}} + I_{\mathrm{q}}^{\mathrm{extra}} \\ - 0 &= - \begin{cases} - I_{\mathrm{q}}^{\mathrm{extra}} - & V_T < V_{\mathrm{hv}}^{\max} \\ - V_T - V_{\mathrm{hv}}^{\max} - & I_{\mathrm{q}}^{\mathrm{extra}} > 0 - \end{cases} \\ - I_L &= I_{L1} - \begin{cases} - 0 & V_M \le V_{L0} \\ - \dfrac{V_M - V_{L0}}{V_{L1} - V_{L0}} & V_{L0} < V_M < V_{L1} \\ - 1 & V_M \ge V_{L1} - \end{cases} \\ - I_{\mathrm{r}} &= I_{\mathrm{p}} - \begin{cases} - 0 & V_T \le V_{A0} \\ - \dfrac{V_T - V_{A0}}{V_{A1} - V_{A0}} & V_{A0} < V_T < V_{A1} \\ - 1 & V_T \ge V_{A1} - \end{cases} \\ - \ell_{\mathrm{p}} &= - \begin{cases} - -R_{\mathrm{p}}^{\max} & I_{\mathrm{p}} \le 0 \\ - -\infty & I_{\mathrm{p}} > 0 - \end{cases} \\ - u_{\mathrm{p}} &= - \begin{cases} - R_{\mathrm{p}}^{\max} & I_{\mathrm{p}} \ge 0 \ \land\ (s_L = 0 \lor I_{\mathrm{p}} < I_L) \\ - 0 & s_L = 1 \ \land\ I_{\mathrm{p}} \ge I_L \\ - \infty & I_{\mathrm{p}} < 0 - \end{cases} -\end{aligned} -``` - -The implemented algebraic residuals use smooth $\text{linseg}$, -$\rho$, and $\sigma$ operators: +The algebraic equations are: ```math \begin{aligned} 0 &= -V_T^2 + V_{\mathrm{r}}^2 + V_{\mathrm{i}}^2 \\ - 0 &= -I_{\mathrm{i}} - I_{\mathrm{q}} + I_{\mathrm{q}}^{\mathrm{extra}} \\ + 0 &= -V_T I_{\mathrm{i}} + - (I_{\mathrm{q}} - I_{\mathrm{q}}^{\mathrm{extra}})V_{\mathrm{r}} + + V_{\mathrm{i}}I_{\mathrm{p}}\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1) \\ 0 &= -I_{\mathrm{q}}^{\mathrm{extra}} - + \rho\!\left(I_{\mathrm{q}}^{\mathrm{extra}} - (V_{\mathrm{hv}}^{\max} - V_T)\right) \\ + + \rho(V_T - V_{\mathrm{hv}}^{\max}) \\ 0 &= -I_L + \text{linseg}(V_M;\ V_{L0},\ V_{L1},\ I_{L1}) \\ - 0 &= -I_{\mathrm{r}} - + I_{\mathrm{p}}\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1) \\ + 0 &= -V_T I_{\mathrm{r}} + + (I_{\mathrm{q}} - I_{\mathrm{q}}^{\mathrm{extra}})V_{\mathrm{i}} + + V_{\mathrm{r}}I_{\mathrm{p}}\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1) \\ 0 &= -\ell_{\mathrm{p}} - R_{\mathrm{p}}^{\max} - (M_{\mathrm{p}} - R_{\mathrm{p}}^{\max})\sigma(I_{\mathrm{p}}) \\ 0 &= -u_{\mathrm{p}} - + \begin{cases} - M_{\mathrm{p}}(1-\sigma(I_{\mathrm{p}})) - + R_{\mathrm{p}}^{\max}\sigma(I_{\mathrm{p}}) - \sigma(I_L - I_{\mathrm{p}}) - & s_L = 1 \\ - M_{\mathrm{p}}(1-\sigma(I_{\mathrm{p}})) - + R_{\mathrm{p}}^{\max}\sigma(I_{\mathrm{p}}) - & s_L = 0 - \end{cases} + + M_{\mathrm{p}}(1-\sigma(I_{\mathrm{p}})) + + R_{\mathrm{p}}^{\max}\sigma(I_{\mathrm{p}}) + \left(1 - s_L + s_L\sigma(I_L - I_{\mathrm{p}})\right) \\ + 0 &= -P_{\mathrm{br}} + + \dfrac{S^{\mathrm{base}}}{S^{\mathrm{sys}}} + \left(V_{\mathrm{r}} I_{\mathrm{r}} + V_{\mathrm{i}} I_{\mathrm{i}}\right) \\ + 0 &= -Q_{\mathrm{br}} + + \dfrac{S^{\mathrm{base}}}{S^{\mathrm{sys}}} + \left(V_{\mathrm{i}} I_{\mathrm{r}} - V_{\mathrm{r}} I_{\mathrm{i}}\right) \end{aligned} ``` @@ -241,8 +187,8 @@ REGCA currents: ```math \begin{aligned} - I_{\mathrm{r}}^{\mathrm{inj}} &:= I_{\mathrm{r}}\dfrac{S^{\mathrm{conv}}}{S^{\mathrm{sys}}} \\ - I_{\mathrm{i}}^{\mathrm{inj}} &:= I_{\mathrm{i}}\dfrac{S^{\mathrm{conv}}}{S^{\mathrm{sys}}} + I_{\mathrm{r}}^{\mathrm{inj}} &:= I_{\mathrm{r}}\dfrac{S^{\mathrm{base}}}{S^{\mathrm{sys}}} \\ + I_{\mathrm{i}}^{\mathrm{inj}} &:= I_{\mathrm{i}}\dfrac{S^{\mathrm{base}}}{S^{\mathrm{sys}}} \end{aligned} ``` @@ -255,57 +201,52 @@ steady-state initial values: ```math \begin{aligned} - V_T &= \sqrt{V_\mathrm{r}^2 + V_\mathrm{i}^2} \\ - I_\mathrm{r0} &= \dfrac{P_0 V_\mathrm{r} + Q_0 V_\mathrm{i}}{V_T^2} - \dfrac{S^\mathrm{sys}}{S^\mathrm{conv}} \\ - I_\mathrm{i0} &= \dfrac{P_0 V_\mathrm{i} - Q_0 V_\mathrm{r}}{V_T^2} - \dfrac{S^\mathrm{sys}}{S^\mathrm{conv}} \\ - V_{M0} &= V_T \\ - I_{L0} &= \text{linseg}(V_T;\ V_{L0},\ V_{L1},\ I_{L1}) \\ - I_\mathrm{p0} &= \dfrac{I_\mathrm{r0}} - {\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1)} \\ - \ell_\mathrm{p0} &= -R_\mathrm{p}^{\max} - - (M_\mathrm{p} - R_\mathrm{p}^{\max})\sigma(I_\mathrm{p0}) \\ - u_\mathrm{p0} &= - \begin{cases} - M_\mathrm{p}(1-\sigma(I_\mathrm{p0})) - + R_\mathrm{p}^{\max}\sigma(I_\mathrm{p0}) - \sigma(I_{L0} - I_\mathrm{p0}) - & s_L = 1 \\ - M_\mathrm{p}(1-\sigma(I_\mathrm{p0})) - + R_\mathrm{p}^{\max}\sigma(I_\mathrm{p0}) - & s_L = 0 - \end{cases} \\ - I_\mathrm{q0}^\mathrm{extra} &= 0 \\ - I_\mathrm{q0} &= -I_\mathrm{i0} + I_\mathrm{q0}^\mathrm{extra} \\ - I_\mathrm{p0}^\mathrm{cmd} &= I_\mathrm{p0} \\ - I_\mathrm{q0}^\mathrm{cmd} &= I_\mathrm{q0} + V_T &= \sqrt{V_{\mathrm{r}}^2 + V_{\mathrm{i}}^2} \\ + V_{M,0} &= V_T \\ + I_{L,0} &= \text{linseg}(V_T;\ V_{L0},\ V_{L1},\ I_{L1}) \\ + I_{\mathrm{p},0} &= V_T^{-1}P_0 + \dfrac{S^{\mathrm{sys}}}{S^{\mathrm{base}}} + \left[\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1)\right]^{-1} \\ + I_{\mathrm{q},0} &= V_T^{-1}Q_0\dfrac{S^{\mathrm{sys}}}{S^{\mathrm{base}}} \\ + \ell_{\mathrm{p},0} &= -R_{\mathrm{p}}^{\max} + - (M_{\mathrm{p}} - R_{\mathrm{p}}^{\max})\sigma(I_{\mathrm{p},0}) \\ + u_{\mathrm{p},0} &= M_{\mathrm{p}}(1-\sigma(I_{\mathrm{p},0})) + + R_{\mathrm{p}}^{\max}\sigma(I_{\mathrm{p},0}) + \left(1 - s_L + s_L\sigma(I_{L,0} - I_{\mathrm{p},0})\right) \\ + I_{\mathrm{q},0}^{\mathrm{extra}} &= 0 \\ + I_{\mathrm{i},0} &= V_T^{-1}\left[ + -I_{\mathrm{q},0}V_{\mathrm{r}} + + V_{\mathrm{i}}I_{\mathrm{p},0}\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1) + \right] \\ + I_{\mathrm{r},0} &= V_T^{-1}\left[ + I_{\mathrm{q},0}V_{\mathrm{i}} + + V_{\mathrm{r}}I_{\mathrm{p},0}\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1) + \right] \\ + P_{\mathrm{br},0} &= \dfrac{S^{\mathrm{base}}}{S^{\mathrm{sys}}} + \left(V_{\mathrm{r}} I_{\mathrm{r},0} + V_{\mathrm{i}} I_{\mathrm{i},0}\right) \\ + Q_{\mathrm{br},0} &= \dfrac{S^{\mathrm{base}}}{S^{\mathrm{sys}}} + \left(V_{\mathrm{i}} I_{\mathrm{r},0} - V_{\mathrm{r}} I_{\mathrm{i},0}\right) \end{aligned} ``` -For normal power-flow starts, $V_T > V_{A1}$, so -$\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1) = 1$ and the -$I_{\mathrm{p0}}$ formula is well defined. - -Initialization should verify: -- $V_T \le V_{\mathrm{hv}}^{\max}$. If $V_T \ge V_{\mathrm{hv}}^{\max}$, - $I_{\mathrm{q0}}^{\mathrm{extra}} = 0$ may not satisfy the HVRCM algebraic - condition, and a nonzero value should be solved or the initialization rejected. -- $\text{linseg}(V_T;\ V_{A0},\ V_{A1},\ 1) > 0$ when - $I_{\mathrm{r0}} \ne 0$. If the LVACM gain is zero, no finite - $I_{\mathrm{p0}}$ can reproduce nonzero initial active current. - -All internal derivatives initialize to zero. +REGCA writes the resolved initial commands to attached `ipcmd` and `iqcmd` +ports. Initialization rejects nonpositive $V_T$, terminal voltage at or above +$V_{\mathrm{hv}}^{\max}$, and nonzero $P_0$ with zero LVACM gain. All internal +derivatives initialize to zero. ## Model Outputs -Real and imaginary injected currents, $I_{\mathrm{r}}$ and -$I_{\mathrm{i}}$, are converter-base algebraic variables. System-base power -outputs use the bus-facing currents: -```math -\begin{aligned} - P &= V_{\mathrm{r}} I_{\mathrm{r}}^{\mathrm{inj}} + V_{\mathrm{i}} I_{\mathrm{i}}^{\mathrm{inj}} \\ - Q &= V_{\mathrm{i}} I_{\mathrm{r}}^{\mathrm{inj}} - V_{\mathrm{r}} I_{\mathrm{i}}^{\mathrm{inj}} -\end{aligned} -``` -Power outputs are positive leaving the converter and entering the bus. +Output | Units | Description | Note +----------------|----------|--------------------------------------|------ +`ir` | [p.u.] | Real current injection | Converter base; exported through `ibranchr` when assigned +`ii` | [p.u.] | Imaginary current injection | Converter base; exported through `ibranchi` when assigned +`p` | [p.u.] | Active-power output | System base; exported through `pbranch` when assigned +`q` | [p.u.] | Reactive-power output | System base; exported through `qbranch` when assigned +`vt` | [p.u.] | Terminal voltage magnitude | +`vm` | [p.u.] | Filtered terminal voltage | +`ip` | [p.u.] | Active-current state | Converter base +`iq` | [p.u.] | Reactive-current state | Converter base +`iqextra` | [p.u.] | HVRCM extra reactive current | Converter base +`il` | [p.u.] | LVPL upper-limit current curve | +`lp` | [p.u./s] | Active-current lower rate bound | +`up` | [p.u./s] | Active-current upper rate bound | diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.cpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.cpp new file mode 100644 index 000000000..b018c4bc3 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.cpp @@ -0,0 +1,27 @@ +/** + * @file Regca.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Non-Enzyme instantiation for the REGCA converter model. + */ + +#include "RegcaImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + template + int Regca::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Regca..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Regca; + template class Regca; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp new file mode 100644 index 000000000..db10ae6c9 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp @@ -0,0 +1,195 @@ +/** + * @file Regca.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Declaration of the REGCA phasor-dynamics converter model. + */ + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + template + class BusBase; + + template + class SignalNode; + } // namespace PhasorDynamics +} // namespace GridKit + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + /// Internal variables of a `Regca` + enum class RegcaInternalVariables : size_t + { + VM, ///< Filtered terminal voltage + IQ, ///< Reactive-current state + IP, ///< Active-current state + VT, ///< Terminal voltage magnitude + II, ///< Imaginary injected current + IQEXTRA, ///< HVRCM extra reactive current + IL, ///< LVPL upper-limit current curve + IR, ///< Real injected current + LP, ///< Active-current lower rate bound + UP, ///< Active-current upper rate bound + PBR, ///< Active-power output on system base + QBR, ///< Reactive-power output on system base + MAXIMUM, + }; + + /// External variables of a `Regca` + enum class RegcaExternalVariables : size_t + { + IPCMD, ///< Active-current command signal + IQCMD, ///< Reactive-current command signal + MAXIMUM, + }; + + template + class Regca : public Component + { + using Component::gridkit_component_id_; + using Component::alpha_; + using Component::f_; + using Component::h_; + using Component::J_; + using Component::J_cols_buffer_; + using Component::J_rows_buffer_; + using Component::J_vals_buffer_; + using Component::nnz_; + using Component::residual_indices_; + using Component::size_; + using Component::tag_; + using Component::time_; + using Component::va_system_base_; + using Component::variable_indices_; + using Component::wb_; + using Component::y_; + using Component::yp_; + + public: + using RealT = typename Component::RealT; + using bus_type = BusBase; + using signal_type = SignalNode; + using model_data_type = RegcaData; + using MonitorT = Model::VariableMonitor; + + Regca(bus_type* bus); + Regca(bus_type* bus, const model_data_type& data); + ~Regca(); + + int setGridKitComponentID(IdxT) override final; + int allocate() override final; + int verify() const override final; + int initialize() override final; + int tagDifferentiable() override final; + int 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*); + __attribute__((always_inline)) inline int evaluateBusResidual( + ScalarT*, ScalarT*, ScalarT*, ScalarT*); + + private: + void initializeParameters(const model_data_type& data); + void initializeMonitor(); + void setDerivedParameters(); + + ScalarT toComponentBase(ScalarT value) const + { + return value * va_system_base_ / va_converter_base_; + } + + ScalarT toSystemBase(ScalarT value) const + { + return value / toComponentBase(static_cast(ONE)); + } + + ScalarT activeCurrentLowerRateBound(ScalarT ip) const; + ScalarT activeCurrentUpperRateBound(ScalarT ip, ScalarT il) const; + + ScalarT& Vr() + { + return bus_->Vr(); + } + + ScalarT& Vi() + { + return bus_->Vi(); + } + + ScalarT& Ir() + { + return bus_->Ir(); + } + + ScalarT& Ii() + { + return bus_->Ii(); + } + + bus_type* bus_{nullptr}; + + RealT P0_{0}; + RealT Q0_{0}; + RealT mva_base_{0}; + RealT Tg_{0}; + RealT TM_{0}; + RealT Rqmax_{0}; + RealT Rqmin_{0}; + RealT Rpmax_{0}; + bool sL_{false}; + RealT IL1_{0}; + RealT VL0_{0}; + RealT VL1_{0}; + RealT VA0_{0}; + RealT VA1_{0}; + RealT Vhvmax_{0}; + IdxT bus_id_{0}; + + IdxT parameter_error_count_{0}; + RealT Mp_{0}; + RealT va_converter_base_{0}; + RealT use_lvpl_{0}; + RealT bypass_lvpl_{1}; + RealT iq_use_upper_{0}; + RealT iq_use_lower_{1}; + + ScalarT ipcmd_set_{0}; + ScalarT iqcmd_set_{0}; + + ComponentSignals signals_; + std::unique_ptr monitor_; + + std::vector ws_; + std::vector ws_indices_; + }; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp new file mode 100644 index 000000000..49d23b861 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp @@ -0,0 +1,81 @@ +/** + * @file RegcaData.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Modeling data for the REGCA converter model. + */ + +#pragma once + +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + /// Parameter keys for the REGCA converter model. + enum class RegcaParameters + { + P0, ///< Initial active power on system base + Q0, ///< Initial reactive power on system base + mva, ///< MVA base of the REGCA model + Tg, ///< Converter current-control lag time constant + TM, ///< Terminal voltage sensor time constant + Rqmax, ///< Reactive-current recovery positive rate limit + Rqmin, ///< Reactive-current recovery negative rate limit + Rpmax, ///< Active-current magnitude recovery rate limit + sL, ///< LVPL switch + IL1, ///< LVPL upper-current ceiling + VL0, ///< LVPL zero-crossing voltage + VL1, ///< LVPL upper breakpoint voltage + VA0, ///< LVACM lower breakpoint voltage + VA1, ///< LVACM upper breakpoint voltage + Vhvmax ///< Terminal-voltage ceiling for HV reactive management + }; + + /// Ports for the REGCA converter model. + enum class RegcaPorts + { + bus, ///< Terminal bus ID + ipcmd, ///< Optional active-current command signal ID + iqcmd, ///< Optional reactive-current command signal ID + ibranchr, ///< Optional real current measurement output signal ID + ibranchi, ///< Optional imaginary current measurement output signal ID + pbranch, ///< Optional active-power measurement output signal ID + qbranch ///< Optional reactive-power measurement output signal ID + }; + + /// Variables available through the monitor interface. + enum class RegcaMonitorableVariables + { + ir, ///< Real current injection on converter base + ii, ///< Imaginary current injection on converter base + p, ///< Active power injection on system base + q, ///< Reactive power injection on system base + vt, ///< Terminal voltage magnitude + vm, ///< Filtered terminal voltage + ip, ///< Active-current state + iq, ///< Reactive-current state + iqextra, ///< HVRCM extra reactive current + il, ///< LVPL upper-limit current curve + lp, ///< Active-current lower rate bound + up ///< Active-current upper rate bound + }; + + template + struct RegcaData : public ComponentData + { + RegcaData() = default; + + using Parameters = RegcaParameters; + using Ports = RegcaPorts; + using MonitorableVariables = RegcaMonitorableVariables; + }; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaDependencyTracking.cpp new file mode 100644 index 000000000..20d688688 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaDependencyTracking.cpp @@ -0,0 +1,27 @@ +/** + * @file RegcaDependencyTracking.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Dependency-tracking instantiations for the REGCA converter model. + */ + +#include "RegcaImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + template + int Regca::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Regca..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Regca; + template class Regca; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp new file mode 100644 index 000000000..8bfb167d7 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp @@ -0,0 +1,110 @@ +/** + * @file RegcaEnzyme.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Enzyme sparse Jacobian for the REGCA converter model. + */ + +#include + +#include "RegcaImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + template + int Regca::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Regca..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + + J_.zeroMatrix(); + if (J_rows_buffer_ == nullptr) + { + // Reserve space for a dense internal block. Enzyme's sparse storage + // keeps only structural nonzeros for each differentiated block. + 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 RegcaT = GridKit::PhasorDynamics::Converter::Regca; + 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_); + + GridKit::Enzyme::Sparse::DhDy::eval(this, + static_cast(bus_->size()), + y_.size(), + (bus_->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + J_); + return 0; + } + + template class Regca; + template class Regca; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp new file mode 100644 index 000000000..57da6794c --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp @@ -0,0 +1,570 @@ +/** + * @file RegcaImpl.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Definition of the REGCA phasor-dynamics converter model. + */ + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + using Log = ::GridKit::Utilities::Logger; + + template + Regca::Regca(bus_type* bus) + : bus_(bus) + { + size_ = static_cast(RegcaInternalVariables::MAXIMUM); + } + + template + Regca::Regca(bus_type* bus, const model_data_type& data) + : bus_(bus), + monitor_(std::make_unique(data)) + { + initializeParameters(data); + initializeMonitor(); + size_ = static_cast(RegcaInternalVariables::MAXIMUM); + } + + template + Regca::~Regca() + { + } + + template + void Regca::setDerivedParameters() + { + Mp_ = static_cast(100.0) * Rpmax_; + use_lvpl_ = ZERO; + bypass_lvpl_ = ONE; + if (sL_) + { + use_lvpl_ = ONE; + bypass_lvpl_ = ZERO; + } + iq_use_upper_ = ZERO; + iq_use_lower_ = ONE; + va_converter_base_ = mva_base_ * static_cast(1.0e6); + } + + template + ScalarT Regca::activeCurrentLowerRateBound(ScalarT ip) const + { + return -Rpmax_ - (Mp_ - Rpmax_) * Math::sigmoid(ip); + } + + template + ScalarT Regca::activeCurrentUpperRateBound(ScalarT ip, ScalarT il) const + { + const ScalarT sigma_ip = Math::sigmoid(ip); + return Mp_ * (ONE - sigma_ip) + + Rpmax_ * sigma_ip * (bypass_lvpl_ + use_lvpl_ * Math::sigmoid(il - ip)); + } + + template + void Regca::initializeParameters(const model_data_type& data) + { + using Params = typename model_data_type::Parameters; + using Ports = typename model_data_type::Ports; + + parameter_error_count_ = 0; + + auto loadRequiredReal = [&](auto key, RealT& target, const char* name) + { + if (!data.parameters.contains(key)) + { + Log::error() << "Regca: missing required parameter '" << name << "'\n"; + ++parameter_error_count_; + return; + } + + const auto& value = data.parameters.at(key); + if (const auto* real_value = std::get_if(&value)) + { + target = *real_value; + } + else if (const auto* index_value = std::get_if(&value)) + { + target = static_cast(*index_value); + } + else + { + Log::error() << "Regca: parameter '" << name << "' must be numeric\n"; + ++parameter_error_count_; + } + }; + + auto loadRequiredSwitch = [&](auto key, bool& target, const char* name) + { + if (!data.parameters.contains(key)) + { + Log::error() << "Regca: missing required parameter '" << name << "'\n"; + ++parameter_error_count_; + return; + } + + const auto& value = data.parameters.at(key); + if (const auto* bool_value = std::get_if(&value)) + { + target = *bool_value; + } + else if (const auto* index_value = std::get_if(&value); + index_value && (*index_value == 0 || *index_value == 1)) + { + target = (*index_value == 1); + } + else + { + Log::error() << "Regca: parameter '" << name << "' must be bool or 0/1\n"; + ++parameter_error_count_; + } + }; + + loadRequiredReal(Params::P0, P0_, "P0"); + loadRequiredReal(Params::Q0, Q0_, "Q0"); + loadRequiredReal(Params::mva, mva_base_, "mva"); + loadRequiredReal(Params::Tg, Tg_, "Tg"); + loadRequiredReal(Params::TM, TM_, "TM"); + loadRequiredReal(Params::Rqmax, Rqmax_, "Rqmax"); + loadRequiredReal(Params::Rqmin, Rqmin_, "Rqmin"); + loadRequiredReal(Params::Rpmax, Rpmax_, "Rpmax"); + loadRequiredSwitch(Params::sL, sL_, "sL"); + loadRequiredReal(Params::IL1, IL1_, "IL1"); + loadRequiredReal(Params::VL0, VL0_, "VL0"); + loadRequiredReal(Params::VL1, VL1_, "VL1"); + loadRequiredReal(Params::VA0, VA0_, "VA0"); + loadRequiredReal(Params::VA1, VA1_, "VA1"); + loadRequiredReal(Params::Vhvmax, Vhvmax_, "Vhvmax"); + + if (data.ports.contains(Ports::bus)) + { + bus_id_ = data.ports.at(Ports::bus); + } + + setDerivedParameters(); + } + + template + const Model::VariableMonitorBase* Regca::getMonitor() const + { + return monitor_.get(); + } + + template + void Regca::initializeMonitor() + { + using Variable = typename model_data_type::MonitorableVariables; + auto index = [](RegcaInternalVariables variable) + { + return static_cast(variable); + }; + + monitor_->set(Variable::ir, [this, index] + { return y_[index(RegcaInternalVariables::IR)]; }); + monitor_->set(Variable::ii, [this, index] + { return y_[index(RegcaInternalVariables::II)]; }); + monitor_->set(Variable::p, [this, index] + { return y_[index(RegcaInternalVariables::PBR)]; }); + monitor_->set(Variable::q, [this, index] + { return y_[index(RegcaInternalVariables::QBR)]; }); + monitor_->set(Variable::vt, [this, index] + { return y_[index(RegcaInternalVariables::VT)]; }); + monitor_->set(Variable::vm, [this, index] + { return y_[index(RegcaInternalVariables::VM)]; }); + monitor_->set(Variable::ip, [this, index] + { return y_[index(RegcaInternalVariables::IP)]; }); + monitor_->set(Variable::iq, [this, index] + { return y_[index(RegcaInternalVariables::IQ)]; }); + monitor_->set(Variable::iqextra, [this, index] + { return y_[index(RegcaInternalVariables::IQEXTRA)]; }); + monitor_->set(Variable::il, [this, index] + { return y_[index(RegcaInternalVariables::IL)]; }); + monitor_->set(Variable::lp, [this, index] + { return y_[index(RegcaInternalVariables::LP)]; }); + monitor_->set(Variable::up, [this, index] + { return y_[index(RegcaInternalVariables::UP)]; }); + } + + template + int Regca::setGridKitComponentID(IdxT component_id) + { + gridkit_component_id_ = component_id; + return 0; + } + + template + int Regca::allocate() + { + size_ = static_cast(RegcaInternalVariables::MAXIMUM); + auto size = static_cast(size_); + + f_.assign(size, ScalarT{0}); + y_.assign(size, ScalarT{0}); + yp_.assign(size, ScalarT{0}); + tag_.assign(size, false); + variable_indices_.resize(size); + residual_indices_.resize(size); + + wb_.assign(2, ScalarT{0}); + h_.assign(2, ScalarT{0}); + + auto signal_size = static_cast(RegcaExternalVariables::MAXIMUM); + ws_.assign(signal_size, ScalarT{0}); + ws_indices_.assign(signal_size, INVALID_INDEX); + + for (IdxT j = 0; j < size_; ++j) + { + this->setVariableIndex(j, j); + this->setResidualIndex(j, j); + } + + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y_[static_cast(RegcaInternalVariables::IR)], + &(this->getVariableIndex(static_cast(RegcaInternalVariables::IR)))); + } + + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y_[static_cast(RegcaInternalVariables::II)], + &(this->getVariableIndex(static_cast(RegcaInternalVariables::II)))); + } + + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y_[static_cast(RegcaInternalVariables::PBR)], + &(this->getVariableIndex(static_cast(RegcaInternalVariables::PBR)))); + } + + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y_[static_cast(RegcaInternalVariables::QBR)], + &(this->getVariableIndex(static_cast(RegcaInternalVariables::QBR)))); + } + + return 0; + } + + template + int Regca::verify() const + { + int ret = static_cast(parameter_error_count_); + + auto check = [&](bool condition, const char* message) + { + if (!condition) + { + Log::error() << "Regca: " << message << '\n'; + ret += 1; + } + }; + + if (bus_ == nullptr) + { + Log::error() << "Regca: bus pointer is null\n"; + ret += 1; + } + + check(mva_base_ > ZERO, "mva must be positive"); + check(Tg_ > ZERO, "Tg must be positive"); + check(TM_ > ZERO, "TM must be positive"); + check(Rpmax_ > ZERO, "Rpmax must be positive"); + check(Rqmin_ < ZERO && ZERO < Rqmax_, "Rqmin < 0 < Rqmax is required"); + check(IL1_ >= ZERO, "IL1 must be non-negative"); + check(ZERO <= VL0_ && VL0_ < VL1_, "VL0/VL1 must satisfy 0 <= VL0 < VL1"); + check(ZERO <= VA0_ && VA0_ < VA1_, "VA0/VA1 must satisfy 0 <= VA0 < VA1"); + check(Vhvmax_ > ZERO, "Vhvmax must be positive"); + + if (signals_.template isAttached()) + { + if (!signals_.template isLinked()) + { + Log::error() << "Regca: ipcmd signal attached with no linked source\n"; + ret += 1; + } + } + + if (signals_.template isAttached()) + { + if (!signals_.template isLinked()) + { + Log::error() << "Regca: iqcmd signal attached with no linked source\n"; + ret += 1; + } + } + + return ret; + } + + template + int Regca::initialize() + { + if (bus_ == nullptr) + { + Log::error() << "Regca: cannot initialize with null bus\n"; + return 1; + } + + if (parameter_error_count_ > 0 || mva_base_ <= ZERO || Tg_ <= ZERO + || TM_ <= ZERO || Rpmax_ <= ZERO + || !(Rqmin_ < ZERO && ZERO < Rqmax_) + || IL1_ < ZERO || !(ZERO <= VL0_ && VL0_ < VL1_) + || !(ZERO <= VA0_ && VA0_ < VA1_) || Vhvmax_ <= ZERO) + { + Log::error() << "Regca: cannot initialize with invalid parameters\n"; + return 1; + } + + const auto VM = static_cast(RegcaInternalVariables::VM); + const auto IQ = static_cast(RegcaInternalVariables::IQ); + const auto IP = static_cast(RegcaInternalVariables::IP); + const auto VT = static_cast(RegcaInternalVariables::VT); + const auto II = static_cast(RegcaInternalVariables::II); + const auto IQEXTRA = static_cast(RegcaInternalVariables::IQEXTRA); + const auto IL = static_cast(RegcaInternalVariables::IL); + const auto IR = static_cast(RegcaInternalVariables::IR); + const auto LP = static_cast(RegcaInternalVariables::LP); + const auto UP = static_cast(RegcaInternalVariables::UP); + const auto PBR = static_cast(RegcaInternalVariables::PBR); + const auto QBR = static_cast(RegcaInternalVariables::QBR); + + const ScalarT vr = Vr(); + const ScalarT vi = Vi(); + const ScalarT vt = std::sqrt(vr * vr + vi * vi); + + if (vt <= ZERO) + { + Log::error() << "Regca: terminal voltage magnitude must be positive at initialization\n"; + return 1; + } + if (vt >= Vhvmax_) + { + Log::error() << "Regca: terminal voltage magnitude must be below Vhvmax at initialization\n"; + return 1; + } + + // REGCA owns the network terminal and establishes the initial + // converter operating point from the power-flow injection. Controller + // command ports may not have been initialized yet, so initialization + // resolves the commands from P0/Q0 and publishes them to attached + // ports below. P0/Q0 are given on the system base; toComponentBase + // converts them to the converter base the internal states use. + const ScalarT lvacm = Math::linseg(vt, VA0_, VA1_, ONE); + + if (P0_ != ZERO && (vt <= VA0_ || lvacm <= ZERO) ) + { + Log::error() << "Regca: LVACM gain is zero with nonzero initial active power\n"; + return 1; + } + + ScalarT ipcmd0{ZERO}; + if (P0_ != ZERO) + { + ipcmd0 = toComponentBase(static_cast(P0_) / vt) / lvacm; + } + const ScalarT iqcmd0 = toComponentBase(static_cast(Q0_) / vt); + + const ScalarT iqextra0{ZERO}; + const ScalarT qnet0 = iqcmd0 - iqextra0; + + y_[VM] = vt; + y_[VT] = vt; + y_[IP] = ipcmd0; + y_[IQ] = iqcmd0; + y_[IQEXTRA] = iqextra0; + y_[IL] = Math::linseg(vt, VL0_, VL1_, IL1_); + y_[IR] = (qnet0 * vi + lvacm * ipcmd0 * vr) / vt; + y_[II] = (-qnet0 * vr + lvacm * ipcmd0 * vi) / vt; + y_[LP] = activeCurrentLowerRateBound(y_[IP]); + y_[UP] = activeCurrentUpperRateBound(y_[IP], y_[IL]); + y_[PBR] = toSystemBase(vr * y_[IR] + vi * y_[II]); + y_[QBR] = toSystemBase(vi * y_[IR] - vr * y_[II]); + + // Retain the resolved commands as the constant source used during + // residual evaluation when no controller drives the command ports, and + // select the reactive-current rate-limit branch from the command sign. + ipcmd_set_ = ipcmd0; + iqcmd_set_ = iqcmd0; + iq_use_upper_ = ZERO; + iq_use_lower_ = ONE; + if (static_cast(iqcmd0) > ZERO) + { + iq_use_upper_ = ONE; + iq_use_lower_ = ZERO; + } + + // Seed attached command nodes with the steady-state values. Controller + // initialization can use these signal values, and unattached ports fall + // back to the constants stored above. + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(ipcmd0); + } + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(iqcmd0); + } + + std::fill(yp_.begin(), yp_.end(), ZERO); + return 0; + } + + template + int Regca::tagDifferentiable() + { + std::fill(tag_.begin(), tag_.end(), false); + tag_[static_cast(RegcaInternalVariables::VM)] = true; + tag_[static_cast(RegcaInternalVariables::IQ)] = true; + tag_[static_cast(RegcaInternalVariables::IP)] = true; + return 0; + } + + template + __attribute__((always_inline)) inline int Regca::evaluateInternalResidual( + ScalarT* y, + ScalarT* yp, + ScalarT* wb, + ScalarT* ws, + ScalarT* f) + { + const auto VM = static_cast(RegcaInternalVariables::VM); + const auto IQ = static_cast(RegcaInternalVariables::IQ); + const auto IP = static_cast(RegcaInternalVariables::IP); + const auto VT = static_cast(RegcaInternalVariables::VT); + const auto II = static_cast(RegcaInternalVariables::II); + const auto IQEXTRA = static_cast(RegcaInternalVariables::IQEXTRA); + const auto IL = static_cast(RegcaInternalVariables::IL); + const auto IR = static_cast(RegcaInternalVariables::IR); + const auto LP = static_cast(RegcaInternalVariables::LP); + const auto UP = static_cast(RegcaInternalVariables::UP); + const auto PBR = static_cast(RegcaInternalVariables::PBR); + const auto QBR = static_cast(RegcaInternalVariables::QBR); + + const auto IPCMD = static_cast(RegcaExternalVariables::IPCMD); + const auto IQCMD = static_cast(RegcaExternalVariables::IQCMD); + + const ScalarT vm = y[VM]; + const ScalarT iq = y[IQ]; + const ScalarT ip = y[IP]; + const ScalarT vt = y[VT]; + const ScalarT ii = y[II]; + const ScalarT iqextra = y[IQEXTRA]; + const ScalarT il = y[IL]; + const ScalarT ir = y[IR]; + const ScalarT lp = y[LP]; + const ScalarT up = y[UP]; + const ScalarT pbr = y[PBR]; + const ScalarT qbr = y[QBR]; + + const ScalarT vm_dot = yp[VM]; + const ScalarT iq_dot = yp[IQ]; + const ScalarT ip_dot = yp[IP]; + + const ScalarT vr = wb[0]; + const ScalarT vi = wb[1]; + + const ScalarT ipcmd = ws[IPCMD]; + const ScalarT iqcmd = ws[IQCMD]; + + const ScalarT iq_error = iqcmd - iq; + const ScalarT ip_error = ipcmd - ip; + + const ScalarT iq_step = + iq_use_upper_ * (iq_error - Math::ramp(iq_error - Tg_ * Rqmax_)) + + iq_use_lower_ * (iq_error + Math::ramp(Tg_ * Rqmin_ - iq_error)); + + const ScalarT ip_step = + Tg_ * lp + Math::ramp(ip_error - Tg_ * lp) - Math::ramp(ip_error - Tg_ * up); + + f[VM] = -TM_ * vm_dot - vm + vt; + f[IQ] = -Tg_ * iq_dot + iq_step; + f[IP] = -Tg_ * ip_dot + ip_step; + const ScalarT lvacm = Math::linseg(vt, VA0_, VA1_, ONE); + const ScalarT qnet = iq - iqextra; + + f[VT] = -vt * vt + vr * vr + vi * vi; + f[II] = -vt * ii - qnet * vr + lvacm * ip * vi; + f[IQEXTRA] = -iqextra + Math::ramp(vt - Vhvmax_); + f[IL] = -il + Math::linseg(vm, VL0_, VL1_, IL1_); + f[IR] = -vt * ir + qnet * vi + lvacm * ip * vr; + f[LP] = -lp + activeCurrentLowerRateBound(ip); + f[UP] = -up + activeCurrentUpperRateBound(ip, il); + f[PBR] = -pbr + toSystemBase(vr * ir + vi * ii); + f[QBR] = -qbr + toSystemBase(vi * ir - vr * ii); + + return 0; + } + + template + __attribute__((always_inline)) inline int Regca::evaluateBusResidual( + ScalarT* y, + [[maybe_unused]] ScalarT* yp, + [[maybe_unused]] ScalarT* wb, + ScalarT* h) + { + const auto II = static_cast(RegcaInternalVariables::II); + const auto IR = static_cast(RegcaInternalVariables::IR); + + h[0] = toSystemBase(y[IR]); + h[1] = toSystemBase(y[II]); + return 0; + } + + template + int Regca::evaluateResidual() + { + const auto IPCMD = static_cast(RegcaExternalVariables::IPCMD); + const auto IQCMD = static_cast(RegcaExternalVariables::IQCMD); + + ws_[IPCMD] = ipcmd_set_; + ws_[IQCMD] = iqcmd_set_; + std::fill(ws_indices_.begin(), ws_indices_.end(), INVALID_INDEX); + + if (signals_.template isAttached()) + { + ws_[IPCMD] = signals_.template readExternalVariable(); + ws_indices_[IPCMD] = + signals_.template readExternalVariableIndex(); + } + + if (signals_.template isAttached()) + { + ws_[IQCMD] = signals_.template readExternalVariable(); + ws_indices_[IQCMD] = + signals_.template readExternalVariableIndex(); + } + + wb_[0] = Vr(); + wb_[1] = Vi(); + + evaluateInternalResidual(y_.data(), yp_.data(), wb_.data(), ws_.data(), f_.data()); + evaluateBusResidual(y_.data(), yp_.data(), wb_.data(), h_.data()); + + Ir() += h_[0]; + Ii() += h_[1]; + + return 0; + } + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Governor/CMakeLists.txt index 7c7269784..fbe71c740 100644 --- a/GridKit/Model/PhasorDynamics/Governor/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/Governor/CMakeLists.txt @@ -4,3 +4,4 @@ # ]] add_subdirectory(Tgov1) +add_subdirectory(HYGOV) diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Governor/HYGOV/CMakeLists.txt new file mode 100644 index 000000000..1719101b0 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/CMakeLists.txt @@ -0,0 +1,54 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +set(_install_headers Hygov.hpp HygovData.hpp) + +if(GRIDKIT_ENABLE_ENZYME) + gridkit_add_library( + phasor_dynamics_governor_hygov + SOURCES HygovEnzyme.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_governor_hygov + SOURCES Hygov.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_governor_hygov_dependency_tracking + SOURCES HygovDependencyTracking.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_governor_hygov) +target_link_libraries( + phasor_dynamics_components_dependency_tracking + INTERFACE GridKit::phasor_dynamics_governor_hygov_dependency_tracking) diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.cpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.cpp new file mode 100644 index 000000000..07d3ea8c1 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.cpp @@ -0,0 +1,27 @@ +/** + * @file Hygov.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Non-Enzyme instantiation for the HYGOV governor model. + */ + +#include "HygovImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + template + int Hygov::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Hygov..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Hygov; + template class Hygov; + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp new file mode 100644 index 000000000..df24d36c4 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp @@ -0,0 +1,152 @@ +/** + * @file Hygov.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Declaration of the HYGOV governor model. + */ + +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + template + class SignalNode; + + namespace Governor + { + /// Internal variables of a `Hygov`. + enum class HygovInternalVariables : size_t + { + XN, ///< Speed lead-lag denominator state + XF, ///< Governor error filter output + C, ///< Desired-gate position + G, ///< Gate position + Q, ///< Turbine flow + OMEGADB, ///< Deadbanded speed deviation + YOMEGA, ///< Lead-lag-conditioned speed signal + EF, ///< Governor error into the filter + FC, ///< Desired-gate derivative target + RC, ///< Rate-limited desired-gate derivative target + PGV, ///< Nonlinear gate-to-power curve output + PGVSAFE, ///< Safe gate-to-power value + H, ///< Turbine head + PMECH, ///< Mechanical-power output + MAXIMUM, + }; + + /// External variables of a `Hygov`. + enum class HygovExternalVariables : size_t + { + OMEGA, ///< Machine speed deviation + PREF, ///< Active-power/load reference + PAUX, ///< Auxiliary power input + MAXIMUM, + }; + + template + class Hygov : 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::variable_indices_; + using Component::wb_; + using Component::y_; + using Component::yp_; + + public: + using RealT = typename Component::RealT; + using signal_type = SignalNode; + using model_data_type = HygovData; + using MonitorT = Model::VariableMonitor; + + Hygov(); + Hygov(const model_data_type& data); + ~Hygov(); + + 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 model_data_type& data); + void initializeMonitor(); + + ScalarT gatePower(ScalarT gate) const; + RealT invertGatePower(RealT pgv) const; + ScalarT toComponentBase(ScalarT value) const; + ScalarT toPmechBase(ScalarT value) const; + + RealT Trate_{0}; + RealT pmech_mva_base_{0}; + RealT Rperm_{0}; + RealT Rtemp_{0}; + RealT Tr_{0}; + RealT Tf_{0}; + RealT Tg_{0}; + RealT Velm_{0}; + RealT Gmax_{0}; + RealT Gmin_{0}; + RealT Tw_{0}; + RealT At_{0}; + RealT Dturb_{0}; + RealT Qnl_{0}; + RealT Tn_{0}; + RealT Tnp_{0}; + RealT leadlag_gain_{0}; + RealT db1_{0}; + RealT db2_{0}; + RealT Hdam_{1}; + std::array Gv_{}; + std::array Pgv_{}; + + IdxT parameter_error_count_{0}; + + ScalarT pref_set_{0}; + ScalarT paux_set_{0}; + + ComponentSignals signals_; + std::unique_ptr monitor_; + + std::vector ws_; + std::vector ws_indices_; + }; + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovData.hpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovData.hpp new file mode 100644 index 000000000..0e026b29e --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovData.hpp @@ -0,0 +1,91 @@ +/** + * @file HygovData.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Modeling data for the HYGOV governor model. + */ + +#pragma once + +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + /// Parameter keys for the HYGOV governor model. + enum class HygovParameters + { + Trate, ///< Turbine-rating power base + mva, ///< MVA base of the pmech signal + Rperm, ///< Permanent droop + Rtemp, ///< Temporary droop + Tr, ///< Temporary-droop reset time constant + Tf, ///< Governor error filter time constant + Tg, ///< Gate servo time constant + Velm, ///< Maximum desired-gate velocity magnitude + Gmax, ///< Maximum desired-gate position + Gmin, ///< Minimum desired-gate position + Tw, ///< Water inertia time constant + At, ///< Turbine gain + Dturb, ///< Turbine damping coefficient + Qnl, ///< No-load flow at nominal head + Tn, ///< Speed lead-lag numerator time constant + Tnp, ///< Speed lead-lag denominator time constant + db1, ///< Type 1 speed deadband threshold + db2, ///< Unsupported mechanical backlash deadband + Hdam, ///< Head available at dam + Gv0, ///< Gate point 0 + Gv1, ///< Gate point 1 + Gv2, ///< Gate point 2 + Gv3, ///< Gate point 3 + Gv4, ///< Gate point 4 + Gv5, ///< Gate point 5 + Pgv0, ///< Power point 0 + Pgv1, ///< Power point 1 + Pgv2, ///< Power point 2 + Pgv3, ///< Power point 3 + Pgv4, ///< Power point 4 + Pgv5 ///< Power point 5 + }; + + /// Ports for the HYGOV governor model. + enum class HygovPorts + { + bus, ///< Optional terminal bus ID for case-format compatibility + speed, ///< Machine speed-deviation signal ID + pmech, ///< Mechanical-power output signal ID + pref, ///< Optional active-power/load reference signal ID + paux ///< Optional auxiliary power input signal ID + }; + + /// Variables available through the monitor interface. + enum class HygovMonitorableVariables + { + pmech, ///< Mechanical power output + leadlag, ///< Speed lead-lag denominator state + filter, ///< Governor error filter output + desiredgate, ///< Desired-gate position + gate, ///< Gate position + flow, ///< Turbine flow + head, ///< Turbine head + pgv ///< Nonlinear gate-to-power curve output + }; + + template + struct HygovData : public ComponentData + { + HygovData() = default; + + using Parameters = HygovParameters; + using Ports = HygovPorts; + using MonitorableVariables = HygovMonitorableVariables; + }; + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovDependencyTracking.cpp new file mode 100644 index 000000000..b77e0d842 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovDependencyTracking.cpp @@ -0,0 +1,27 @@ +/** + * @file HygovDependencyTracking.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Dependency-tracking instantiations for the HYGOV governor model. + */ + +#include "HygovImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + template + int Hygov::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Hygov..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Hygov; + template class Hygov; + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovEnzyme.cpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovEnzyme.cpp new file mode 100644 index 000000000..6e3190771 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovEnzyme.cpp @@ -0,0 +1,110 @@ +/** + * @file HygovEnzyme.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Enzyme sparse Jacobian for the HYGOV governor model. + */ + +#include + +#include "HygovImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + namespace + { + template + RealT deadband1Derivative(RealT x, RealT lower, RealT upper) + { + static constexpr RealT MU = static_cast(240.0); + + const RealT lower_sig = Math::sigmoid(lower - x); + const RealT upper_sig = Math::sigmoid(x - upper); + return lower_sig + upper_sig + + x * MU + * (upper_sig * (ONE - upper_sig) + - lower_sig * (ONE - lower_sig)); + } + } // namespace + + template + int Hygov::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Hygov..." << 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::Governor::Hygov; + 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_); + + const auto OMEGADB = static_cast(HygovInternalVariables::OMEGADB); + const auto EF = static_cast(HygovInternalVariables::EF); + const auto PMECH = static_cast(HygovInternalVariables::PMECH); + const auto G = static_cast(HygovInternalVariables::G); + + const auto OMEGA = static_cast(HygovExternalVariables::OMEGA); + const auto PREF = static_cast(HygovExternalVariables::PREF); + const auto PAUX = static_cast(HygovExternalVariables::PAUX); + + IdxT nnz = 0; + auto storeSignalDerivative = [&](size_t row, size_t col, RealT value) + { + if (ws_indices_[col] == INVALID_INDEX) + { + return; + } + + J_rows_buffer_[nnz] = residual_indices_.at(row); + J_cols_buffer_[nnz] = ws_indices_[col]; + J_vals_buffer_[nnz] = value; + ++nnz; + }; + + storeSignalDerivative(OMEGADB, + OMEGA, + deadband1Derivative(static_cast(ws_[OMEGA]), + -db1_, + db1_)); + storeSignalDerivative(EF, PREF, ONE); + storeSignalDerivative(EF, PAUX, ONE); + storeSignalDerivative(PMECH, + OMEGA, + -Dturb_ * static_cast(y_[G]) * Trate_ / pmech_mva_base_); + + J_.setValues(1.0, J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, nnz); + return 0; + } + + template class Hygov; + template class Hygov; + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp new file mode 100644 index 000000000..0dae89b13 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp @@ -0,0 +1,541 @@ +/** + * @file HygovImpl.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Definition of the HYGOV governor model. + */ + +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + using Log = ::GridKit::Utilities::Logger; + + template + Hygov::Hygov() + { + size_ = static_cast(HygovInternalVariables::MAXIMUM); + } + + template + Hygov::Hygov(const model_data_type& data) + : monitor_(std::make_unique(data)) + { + initModelParams(data); + initializeMonitor(); + size_ = static_cast(HygovInternalVariables::MAXIMUM); + } + + template + Hygov::~Hygov() + { + } + + template + ScalarT Hygov::toComponentBase(ScalarT value) const + { + return value * pmech_mva_base_ / Trate_; + } + + template + ScalarT Hygov::toPmechBase(ScalarT value) const + { + return value * Trate_ / pmech_mva_base_; + } + + template + void Hygov::initModelParams(const model_data_type& data) + { + using Params = typename model_data_type::Parameters; + + parameter_error_count_ = 0; + + auto loadRequiredReal = [&](auto key, RealT& target, const char* name) + { + if (!data.parameters.contains(key)) + { + Log::error() << "Hygov: 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() << "Hygov: parameter '" << name << "' must be numeric\n"; + ++parameter_error_count_; + } + }; + + auto loadOptionalReal = [&](auto key, RealT& target, const char* name) + { + if (!data.parameters.contains(key)) + { + 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() << "Hygov: parameter '" << name << "' must be numeric\n"; + ++parameter_error_count_; + } + }; + + loadRequiredReal(Params::Trate, Trate_, "Trate"); + loadRequiredReal(Params::mva, pmech_mva_base_, "mva"); + loadRequiredReal(Params::Rperm, Rperm_, "Rperm"); + loadRequiredReal(Params::Rtemp, Rtemp_, "Rtemp"); + loadRequiredReal(Params::Tr, Tr_, "Tr"); + loadRequiredReal(Params::Tf, Tf_, "Tf"); + loadRequiredReal(Params::Tg, Tg_, "Tg"); + loadRequiredReal(Params::Velm, Velm_, "Velm"); + loadRequiredReal(Params::Gmax, Gmax_, "Gmax"); + loadRequiredReal(Params::Gmin, Gmin_, "Gmin"); + loadRequiredReal(Params::Tw, Tw_, "Tw"); + loadRequiredReal(Params::At, At_, "At"); + loadRequiredReal(Params::Dturb, Dturb_, "Dturb"); + loadRequiredReal(Params::Qnl, Qnl_, "Qnl"); + loadRequiredReal(Params::Tn, Tn_, "Tn"); + loadRequiredReal(Params::Tnp, Tnp_, "Tnp"); + loadRequiredReal(Params::db1, db1_, "db1"); + loadOptionalReal(Params::db2, db2_, "db2"); + loadRequiredReal(Params::Hdam, Hdam_, "Hdam"); + + loadRequiredReal(Params::Gv0, Gv_[0], "Gv0"); + loadRequiredReal(Params::Gv1, Gv_[1], "Gv1"); + loadRequiredReal(Params::Gv2, Gv_[2], "Gv2"); + loadRequiredReal(Params::Gv3, Gv_[3], "Gv3"); + loadRequiredReal(Params::Gv4, Gv_[4], "Gv4"); + loadRequiredReal(Params::Gv5, Gv_[5], "Gv5"); + loadRequiredReal(Params::Pgv0, Pgv_[0], "Pgv0"); + loadRequiredReal(Params::Pgv1, Pgv_[1], "Pgv1"); + loadRequiredReal(Params::Pgv2, Pgv_[2], "Pgv2"); + loadRequiredReal(Params::Pgv3, Pgv_[3], "Pgv3"); + loadRequiredReal(Params::Pgv4, Pgv_[4], "Pgv4"); + loadRequiredReal(Params::Pgv5, Pgv_[5], "Pgv5"); + + const bool source_default_curve = + std::all_of(Gv_.begin(), Gv_.end(), [](RealT value) + { return value == ZERO; }) + && std::all_of(Pgv_.begin(), Pgv_.end(), [](RealT value) + { return value == ZERO; }); + if (source_default_curve) + { + Gv_ = {ZERO, static_cast(0.2), static_cast(0.4), static_cast(0.6), static_cast(0.8), ONE}; + Pgv_ = Gv_; + } + + leadlag_gain_ = (Tnp_ > ZERO) ? Tn_ / Tnp_ : ZERO; + } + + template + const Model::VariableMonitorBase* Hygov::getMonitor() const + { + return monitor_.get(); + } + + template + void Hygov::initializeMonitor() + { + using Variable = typename model_data_type::MonitorableVariables; + auto index = [](HygovInternalVariables variable) + { + return static_cast(variable); + }; + + monitor_->set(Variable::pmech, [this, index] + { return y_[index(HygovInternalVariables::PMECH)]; }); + monitor_->set(Variable::leadlag, [this, index] + { return y_[index(HygovInternalVariables::XN)]; }); + monitor_->set(Variable::filter, [this, index] + { return y_[index(HygovInternalVariables::XF)]; }); + monitor_->set(Variable::desiredgate, [this, index] + { return y_[index(HygovInternalVariables::C)]; }); + monitor_->set(Variable::gate, [this, index] + { return y_[index(HygovInternalVariables::G)]; }); + monitor_->set(Variable::flow, [this, index] + { return y_[index(HygovInternalVariables::Q)]; }); + monitor_->set(Variable::head, [this, index] + { return y_[index(HygovInternalVariables::H)]; }); + monitor_->set(Variable::pgv, [this, index] + { return y_[index(HygovInternalVariables::PGV)]; }); + } + + template + int Hygov::setGridKitComponentID(IdxT component_id) + { + gridkit_component_id_ = component_id; + return 0; + } + + template + int Hygov::allocate() + { + size_ = static_cast(HygovInternalVariables::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_.clear(); + + auto signal_size = static_cast(HygovExternalVariables::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(HygovInternalVariables::PMECH)], + &(this->getVariableIndex(static_cast(HygovInternalVariables::PMECH)))); + } + + return 0; + } + + template + int Hygov::verify() const + { + int ret = static_cast(parameter_error_count_); + + auto check = [&](bool condition, const char* message) + { + if (!condition) + { + Log::error() << "Hygov: " << message << '\n'; + ret += 1; + } + }; + + check(Trate_ > ZERO, "Trate must be positive"); + check(pmech_mva_base_ > ZERO, "mva must be positive"); + check(Rperm_ > ZERO, "Rperm must be positive"); + check(Rtemp_ > ZERO, "Rtemp must be positive"); + check(Tr_ > ZERO, "Tr must be positive"); + check(Tf_ > ZERO, "Tf must be positive"); + check(Tg_ >= ZERO, "Tg must be non-negative"); + check(Tw_ > ZERO, "Tw must be positive"); + check(Tn_ >= ZERO, "Tn must be non-negative"); + check(Tnp_ >= ZERO, "Tnp must be non-negative"); + check(Tnp_ > ZERO || Tn_ == ZERO, "Tn must be zero when Tnp is zero"); + check(Velm_ > ZERO, "Velm must be positive"); + check(Gmin_ <= Gmax_, "Gmin must be less than or equal to Gmax"); + check(At_ > ZERO, "At must be positive"); + check(Dturb_ >= ZERO, "Dturb must be non-negative"); + check(Qnl_ >= ZERO, "Qnl must be non-negative"); + check(db1_ >= ZERO, "db1 must be non-negative"); + check(db2_ == ZERO, "nonzero db2 is not supported"); + check(Hdam_ > ZERO, "Hdam must be positive"); + + for (size_t i = 1; i < Gv_.size(); ++i) + { + check(Gv_[i - 1] < Gv_[i], "Gv points must be strictly increasing"); + check(Pgv_[i - 1] <= Pgv_[i], "Pgv points must be non-decreasing"); + } + check(Pgv_[0] >= ZERO, "Pgv0 must be non-negative"); + + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "Hygov: omega signal attached with no linked source\n"; + ret += 1; + } + + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "Hygov: pref signal attached with no linked source\n"; + ret += 1; + } + + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "Hygov: paux signal attached with no linked source\n"; + ret += 1; + } + + return ret; + } + + template + ScalarT Hygov::gatePower(ScalarT gate) const + { + return ScalarT{Pgv_[0]} + + Math::linseg(gate, Gv_[0], Gv_[1], Pgv_[1] - Pgv_[0]) + + Math::linseg(gate, Gv_[1], Gv_[2], Pgv_[2] - Pgv_[1]) + + Math::linseg(gate, Gv_[2], Gv_[3], Pgv_[3] - Pgv_[2]) + + Math::linseg(gate, Gv_[3], Gv_[4], Pgv_[4] - Pgv_[3]) + + Math::linseg(gate, Gv_[4], Gv_[5], Pgv_[5] - Pgv_[4]); + } + + template + typename Hygov::RealT Hygov::invertGatePower( + typename Hygov::RealT pgv) const + { + static constexpr RealT tol = static_cast(1.0e-10); + + if (std::abs(pgv - Pgv_[0]) <= tol) + { + return Gv_[0]; + } + + for (size_t i = 0; i < 5; ++i) + { + if (Pgv_[i + 1] <= Pgv_[i]) + { + continue; + } + + if (Pgv_[i] - tol <= pgv && pgv <= Pgv_[i + 1] + tol) + { + const RealT fraction = (pgv - Pgv_[i]) / (Pgv_[i + 1] - Pgv_[i]); + return Gv_[i] + fraction * (Gv_[i + 1] - Gv_[i]); + } + } + + return std::numeric_limits::quiet_NaN(); + } + + template + int Hygov::initialize() + { + if (parameter_error_count_ > 0 || verify() > 0) + { + Log::error() << "Hygov: cannot initialize with invalid configuration\n"; + return 1; + } + + const auto XN = static_cast(HygovInternalVariables::XN); + const auto XF = static_cast(HygovInternalVariables::XF); + const auto C = static_cast(HygovInternalVariables::C); + const auto G = static_cast(HygovInternalVariables::G); + const auto Q = static_cast(HygovInternalVariables::Q); + const auto OMEGADB = static_cast(HygovInternalVariables::OMEGADB); + const auto YOMEGA = static_cast(HygovInternalVariables::YOMEGA); + const auto EF = static_cast(HygovInternalVariables::EF); + const auto FC = static_cast(HygovInternalVariables::FC); + const auto RC = static_cast(HygovInternalVariables::RC); + const auto PGV = static_cast(HygovInternalVariables::PGV); + const auto PGVSAFE = static_cast(HygovInternalVariables::PGVSAFE); + const auto H = static_cast(HygovInternalVariables::H); + const auto PMECH = static_cast(HygovInternalVariables::PMECH); + + ScalarT omega0{ZERO}; + if (signals_.template isAttached()) + { + omega0 = signals_.template readExternalVariable(); + } + + paux_set_ = ScalarT{ZERO}; + if (signals_.template isAttached()) + { + paux_set_ = signals_.template readExternalVariable(); + } + + const ScalarT pmech0 = toComponentBase(y_[PMECH]); + y_[H] = Hdam_; + y_[Q] = Qnl_ + pmech0 / (At_ * y_[H]); + y_[PGV] = y_[Q] / std::sqrt(y_[H]); + y_[PGVSAFE] = Math::max(y_[PGV], static_cast(0.01)); + + const RealT gate0 = invertGatePower(static_cast(y_[PGV])); + if (std::isnan(gate0)) + { + Log::error() << "Hygov: initial Pgv is outside the invertible gate curve\n"; + return 1; + } + + y_[G] = gate0; + y_[C] = y_[G]; + + if (y_[C] < Gmin_ || y_[C] > Gmax_) + { + Log::error() << "Hygov: initialized gate is outside Gmin/Gmax\n"; + return 1; + } + + y_[OMEGADB] = Math::deadband1(omega0, -db1_, db1_); + y_[XN] = y_[OMEGADB]; + y_[YOMEGA] = y_[OMEGADB]; + y_[XF] = ZERO; + y_[EF] = ZERO; + y_[FC] = ZERO; + y_[RC] = ZERO; + y_[PMECH] = toPmechBase(pmech0); + + pref_set_ = y_[EF] - paux_set_ + y_[YOMEGA] + Rperm_ * y_[C]; + if (signals_.template isAttached()) + { + const ScalarT pref0 = signals_.template readExternalVariable(); + const RealT pref_err = static_cast(pref0 - pref_set_); + if (std::abs(pref_err) > static_cast(1.0e-10)) + { + Log::error() << "Hygov: pref initial condition is not steady state\n"; + return 1; + } + pref_set_ = pref0; + } + + std::fill(yp_.begin(), yp_.end(), ZERO); + return 0; + } + + template + int Hygov::tagDifferentiable() + { + std::fill(tag_.begin(), tag_.end(), false); + tag_[static_cast(HygovInternalVariables::XN)] = (Tnp_ > ZERO); + tag_[static_cast(HygovInternalVariables::XF)] = true; + tag_[static_cast(HygovInternalVariables::C)] = true; + tag_[static_cast(HygovInternalVariables::G)] = (Tg_ > ZERO); + tag_[static_cast(HygovInternalVariables::Q)] = true; + return 0; + } + + template + __attribute__((always_inline)) inline int Hygov::evaluateInternalResidual( + ScalarT* y, + ScalarT* yp, + [[maybe_unused]] ScalarT* wb, + ScalarT* ws, + ScalarT* f) + { + const auto XN = static_cast(HygovInternalVariables::XN); + const auto XF = static_cast(HygovInternalVariables::XF); + const auto C = static_cast(HygovInternalVariables::C); + const auto G = static_cast(HygovInternalVariables::G); + const auto Q = static_cast(HygovInternalVariables::Q); + const auto OMEGADB = static_cast(HygovInternalVariables::OMEGADB); + const auto YOMEGA = static_cast(HygovInternalVariables::YOMEGA); + const auto EF = static_cast(HygovInternalVariables::EF); + const auto FC = static_cast(HygovInternalVariables::FC); + const auto RC = static_cast(HygovInternalVariables::RC); + const auto PGV = static_cast(HygovInternalVariables::PGV); + const auto PGVSAFE = static_cast(HygovInternalVariables::PGVSAFE); + const auto H = static_cast(HygovInternalVariables::H); + const auto PMECH = static_cast(HygovInternalVariables::PMECH); + + const auto OMEGA = static_cast(HygovExternalVariables::OMEGA); + const auto PREF = static_cast(HygovExternalVariables::PREF); + const auto PAUX = static_cast(HygovExternalVariables::PAUX); + + const ScalarT xn = y[XN]; + const ScalarT xf = y[XF]; + const ScalarT c = y[C]; + const ScalarT g = y[G]; + const ScalarT q = y[Q]; + const ScalarT omegadb = y[OMEGADB]; + const ScalarT yomega = y[YOMEGA]; + const ScalarT ef = y[EF]; + const ScalarT fc = y[FC]; + const ScalarT rc = y[RC]; + const ScalarT pgv = y[PGV]; + const ScalarT pgv_safe = y[PGVSAFE]; + const ScalarT head = y[H]; + const ScalarT pmech = y[PMECH]; + + const ScalarT omega = ws[OMEGA]; + const ScalarT pref = ws[PREF]; + const ScalarT paux = ws[PAUX]; + + f[XN] = -Tnp_ * yp[XN] - xn + omegadb; + f[XF] = -Tf_ * yp[XF] - xf + ef; + f[FC] = -Rtemp_ * Tr_ * fc + xf + Tr_ * yp[XF]; + f[C] = -yp[C] + Math::antiwindup(c, rc, Gmin_, Gmax_); + f[G] = -Tg_ * yp[G] - g + c; + f[Q] = -Tw_ * yp[Q] + Hdam_ - head; + + f[OMEGADB] = -omegadb + Math::deadband1(omega, -db1_, db1_); + f[YOMEGA] = -yomega + xn + leadlag_gain_ * (omegadb - xn); + f[EF] = -ef + pref + paux - yomega - Rperm_ * c; + f[RC] = -rc + Math::clamp(fc, -Velm_, Velm_); + f[PGV] = -pgv + gatePower(g); + f[PGVSAFE] = -pgv_safe + Math::max(pgv, static_cast(0.01)); + f[H] = -head + (q / pgv_safe) * (q / pgv_safe); + const ScalarT pmech_output = At_ * head * (q - Qnl_) - Dturb_ * omega * g; + f[PMECH] = -pmech + toPmechBase(pmech_output); + + return 0; + } + + template + int Hygov::evaluateResidual() + { + const auto OMEGA = static_cast(HygovExternalVariables::OMEGA); + const auto PREF = static_cast(HygovExternalVariables::PREF); + const auto PAUX = static_cast(HygovExternalVariables::PAUX); + + ws_[OMEGA] = ZERO; + ws_[PREF] = pref_set_; + ws_[PAUX] = paux_set_; + std::fill(ws_indices_.begin(), ws_indices_.end(), INVALID_INDEX); + + if (signals_.template isAttached()) + { + ws_[OMEGA] = signals_.template readExternalVariable(); + ws_indices_[OMEGA] = + signals_.template readExternalVariableIndex(); + } + + if (signals_.template isAttached()) + { + ws_[PREF] = signals_.template readExternalVariable(); + ws_indices_[PREF] = + signals_.template readExternalVariableIndex(); + } + + if (signals_.template isAttached()) + { + ws_[PAUX] = signals_.template readExternalVariable(); + ws_indices_[PAUX] = + signals_.template readExternalVariableIndex(); + } + + evaluateInternalResidual(y_.data(), yp_.data(), wb_.data(), ws_.data(), f_.data()); + return 0; + } + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/README.md b/GridKit/Model/PhasorDynamics/Governor/HYGOV/README.md new file mode 100644 index 000000000..d0957ecdd --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/README.md @@ -0,0 +1,291 @@ +# **Hydro Turbine-Governor Model (HYGOV)** + +HYGOV is a hydro turbine-governor model with a temporary-droop governor, gate +servo, and nonlinear single-penstock water-column turbine. In GridKit it is +represented as a governor model that reads machine speed deviation and supplies +mechanical power to the machine. + +Notes: +- HYGOV equations use the `Trate` base; the `pmech` signal uses `mva` and is converted at the signal boundary. +- GridKit requires `Trate > 0`; for PowerWorld `Trate = 0` records, set `Trate` to the connected machine `mva` before running. +- Aside from `db2`, model parameters are required input values; the table values are typical values, not parser defaults. +- HYGOV uses `db1` with CommonMath `deadband1`; HYGOVD-only `dbL`/`dbH` and nonzero `db2` are not modeled. +- PowerWorld fields `Ttur`, `Eps`, Kaplan blade-servo points, `Bmax`, and `Tblade` are not part of this implementation. +- Permanent-droop feedback uses desired gate $c$ (State 2), not gate $g$. + +## Block Diagram + +Standard HYGOV block diagram. + +
+ + + Figure 1: Governor HYGOV model. Figure courtesy of [PowerWorld](https://www.powerworld.com/WebHelp/) +
+ +## Model Parameters + +Symbol | Units | JSON | Description | Typical Value | Note +--------------------------------|----------|----------|----------------------------------------------|---------------|------ +$P^{\mathrm{rate}}$ | [MW] | `Trate` | Component-base turbine rating | 100.0 | Required positive value +$S^{\mathrm{pmech}}$ | [MVA] | `mva` | Mechanical-power signal base | 100.0 | Use connected machine `mva` +$R_{\mathrm{perm}}$ | [p.u.] | `Rperm` | Permanent droop | 0.04 | Source diagram label: `R` +$R_{\mathrm{temp}}$ | [p.u.] | `Rtemp` | Temporary droop | 0.3 | Source diagram label: `r` +$T_r$ | [sec] | `Tr` | Temporary-droop reset time constant | 5.0 | +$T_f$ | [sec] | `Tf` | Governor error filter time constant | 0.05 | State 1 +$T_g$ | [sec] | `Tg` | Gate servo time constant | 0.5 | State 3; if zero, $g$ is algebraic +$V_{\mathrm{elm}}$ | [p.u./s] | `Velm` | Maximum desired-gate velocity magnitude | 0.2 | Symmetric rate limit on State 2 +$G^{\max}$ | [p.u.] | `Gmax` | Maximum desired-gate position | 1.0 | +$G^{\min}$ | [p.u.] | `Gmin` | Minimum desired-gate position | 0.0 | +$T_w$ | [sec] | `Tw` | Water inertia time constant | 1.0 | State 4 +$A_t$ | [p.u.] | `At` | Turbine gain | 1.2 | +$D_{\mathrm{turb}}$ | [p.u.] | `Dturb` | Turbine damping coefficient | 0.5 | Multiplied by speed deviation and gate +$q_{\mathrm{NL}}$ | [p.u.] | `Qnl` | No-load flow at nominal head | 0.05 | +$T_n$ | [sec] | `Tn` | Speed lead-lag numerator time constant | 0.0 | Source text labels this "lag" +$T_{np}$ | [sec] | `Tnp` | Speed lead-lag denominator time constant | 0.0 | Source text labels this "lead" +$D_{\omega}$ | [p.u.] | `db1` | Type 1 speed deadband threshold | 0.0 | Source data may provide this in Hz +$H_{\mathrm{dam}}$ | [p.u.] | `Hdam` | Head available at dam | 1.0 | + +The nonlinear gate-to-power curve $N_{\mathrm{GV}}$ is represented by six source points: + +Symbol | Units | JSON | Description | Typical Value | Note +--------------------------------|--------|----------|----------------------------------------------|---------------|------ +$G_V^{(k)}$ | [p.u.] | `Gv0`-`Gv5` | Gate point $k$ of the gain curve | 0.0 | $k=0,\ldots,5$ +$P_{\mathrm{GV}}^{(k)}$ | [p.u.] | `Pgv0`-`Pgv5` | Power point $k$ of the gain curve | 0.0 | $k=0,\ldots,5$ + +If all six `Gv` and all six `Pgv` source values are zero, the source data is treated as +omitting the nonlinear curve. The effective curve used by the equations is the identity +curve with points `(0,0), (0.2,0.2), ..., (1,1)`. + +### Parameter Validation + +Invalid HYGOV parameter sets are rejected by the following checks. If source data preprocessing swaps gate limits, expands gate limits to include the initialized gate, adjusts small positive time constants, converts `db1` from Hz to p.u. speed, or replaces an all-zero source gate curve with the identity curve, apply these checks to the effective values used by the equations. Source data with nonzero `db2` is outside the supported equations below. + +The required checks are: + +```math +\begin{aligned} + &P^{\mathrm{rate}} > 0 \\ + &S^{\mathrm{pmech}} > 0 \\ + &R_{\mathrm{perm}} > 0 \\ + &R_{\mathrm{temp}} > 0 \\ + &T_r > 0 \\ + &T_f > 0,\quad T_g \ge 0 \\ + &T_w > 0 \\ + &T_n \ge 0,\quad T_{np} \ge 0,\quad T_{np}=0 \Rightarrow T_n=0 \\ + &V_{\mathrm{elm}} > 0 \\ + &G^{\min} \le G^{\max} \\ + &A_t > 0 \\ + &D_{\mathrm{turb}} \ge 0 \\ + &q_{\mathrm{NL}} \ge 0 \\ + &D_{\omega} \ge 0 \\ + &H_{\mathrm{dam}} > 0 \\ + &G_V^{(0)} < G_V^{(1)} < \cdots < G_V^{(5)} \\ + &0 \le P_{\mathrm{GV}}^{(0)} \le P_{\mathrm{GV}}^{(1)} \le \cdots \le P_{\mathrm{GV}}^{(5)} +\end{aligned} +``` + +Initialization also requires $N_{\mathrm{GV}}^{-1}$ to be single-valued at the solved initial operating point. + +### Model Derived Parameters + +```math +\begin{aligned} + k_n &= + \begin{cases} + T_n/T_{np} & T_{np} > 0 \\ + 0 & T_{np} = 0 + \end{cases} \\ + S_{\mathrm{hygov}}^{\mathrm{base}} + &= P^{\mathrm{rate}} +\end{aligned} +``` + +Here $S^{\mathrm{pmech}}$ is the `mva` parameter. The `pmech` signal uses +that base, while HYGOV turbine equations use the HYGOV component base: + +```math +\begin{aligned} + P^{\mathrm{hygov}} + &= P^{\mathrm{pmech}}\dfrac{S^{\mathrm{pmech}}}{S_{\mathrm{hygov}}^{\mathrm{base}}} \\ + P^{\mathrm{pmech}} + &= P^{\mathrm{hygov}}\dfrac{S_{\mathrm{hygov}}^{\mathrm{base}}}{S^{\mathrm{pmech}}} +\end{aligned} +``` + +The nonlinear gate-to-power function uses GridKit's smooth [Linear Segment](../../../../CommonMath.md#derived-functions) helper: + +```math +\begin{aligned} + N_{\mathrm{GV}}(x) + &= + P_{\mathrm{GV}}^{(0)} + + \sum_{k\in\{0,\ldots,4\}} + \text{linseg}\!\left( + x;\, + G_V^{(k)},\, + G_V^{(k+1)},\, + P_{\mathrm{GV}}^{(k+1)} - P_{\mathrm{GV}}^{(k)} + \right) +\end{aligned} +``` + +## Model Variables + +### Internal Variables + +#### Differential + +Symbol | Units | Description | Note +------------------------|--------|-------------------------------------|------ +$x_n$ | [p.u.] | Speed lead-lag denominator state | Not circled in Fig. 1; realizes the `Tn`/`Tnp` block +$x_f$ | [p.u.] | Governor error filter output | State 1 in Fig. 1 +$c$ | [p.u.] | Desired-gate position | State 2 in Fig. 1 +$g$ | [p.u.] | Gate position | State 3 in Fig. 1; algebraic when $T_g = 0$ +$q$ | [p.u.] | Turbine flow | State 4 in Fig. 1 + +#### Algebraic + +Symbol | Units | Description | Note +----------------------------------|----------|-------------------------------------|------ +$\omega_{\mathrm{db}}$ | [p.u.] | Type 1 deadbanded speed deviation | Defined by CommonMath `deadband1` +$y_{\omega}$ | [p.u.] | Lead-lag-conditioned speed signal | Output of the speed lead-lag +$e_f$ | [p.u.] | Governor error into the filter | Reference path less conditioned speed and permanent-droop feedback +$f_c$ | [p.u./s] | Desired-gate derivative target | Before rate and position limits +$r_c$ | [p.u./s] | Rate-limited desired-gate derivative target | Limited by $\pm V_{\mathrm{elm}}$ +$P_{\mathrm{GV}}$ | [p.u.] | Nonlinear gate-to-power curve output | $N_{\mathrm{GV}}(g)$ +$P_{\mathrm{GV}}^{\mathrm{safe}}$ | [p.u.] | Safe gate-to-power value | Lower bounded by 0.01 +$H$ | [p.u.] | Turbine head | $(q/P_{\mathrm{GV}}^{\mathrm{safe}})^2$ +$P_{\mathrm{mech}}$ | [p.u.] | Mechanical power to generator | `mva` signal; converted at the HYGOV boundary + +### External Variables + +#### Differential +None. + +#### Algebraic + +Symbol | Units | Description | Note +--------------------------------|--------|--------------------------------|------ +$\omega$ | [p.u.] | Machine speed deviation | Read from a machine model +$P_{\mathrm{ref}}$ | [p.u.] | Active-power/load reference | HYGOV component base; external setpoint or constant parameter; source label: `Pref` +$P_{\mathrm{aux}}$ | [p.u.] | Auxiliary power input | HYGOV component base; optional, defaults to zero; source label: `Paux` + +## Model Equations + +The desired-gate transfer block before limits is: + +```math +\begin{aligned} + \dfrac{c(s)}{x_f(s)} + &= \dfrac{1+sT_r}{R_{\mathrm{temp}}T_rs} +\end{aligned} +``` + +The realization below forms the derivative target $f_c$, limits it to $r_c$, and applies anti-windup at the desired-gate limits. + +### Differential Equations + +```math +\begin{aligned} + 0 &= -T_{np}\dot x_n - x_n + \omega_{\mathrm{db}} \\ + 0 &= -T_f\dot x_f - x_f + e_f \\ + 0 &= -R_{\mathrm{temp}}T_r f_c + x_f + T_r\dot x_f \\ + 0 &= + -\dot c + + \text{antiwindup}\!\left( + c, + r_c, + G^{\min}, + G^{\max} + \right) \\ + 0 &= -T_g\dot g - g + c \\ + 0 &= -T_w\dot q + H_{\mathrm{dam}} - H +\end{aligned} +``` + +CommonMath defines the [Anti-Windup](../../../../CommonMath.md#anti-windup-indicator) target and smooth approximation. + +### Algebraic Equations + +```math +\begin{aligned} + 0 &= -\omega_{\mathrm{db}} + + \text{deadband1}\!\left(\omega, -D_{\omega}, D_{\omega}\right) \\ + 0 &= -y_{\omega} + + x_n + + k_n\left(\omega_{\mathrm{db}} - x_n\right) \\ + 0 &= -e_f + P_{\mathrm{ref}} + P_{\mathrm{aux}} - y_{\omega} - R_{\mathrm{perm}}c \\ + 0 &= -r_c + \text{clamp}\!\left(f_c, -V_{\mathrm{elm}}, V_{\mathrm{elm}}\right) \\[0.5ex] + 0 &= -P_{\mathrm{GV}} + N_{\mathrm{GV}}(g) \\ + 0 &= -P_{\mathrm{GV}}^{\mathrm{safe}} + \text{max}\!\left(P_{\mathrm{GV}}, 0.01\right) \\ + 0 &= -H + \left(\dfrac{q}{P_{\mathrm{GV}}^{\mathrm{safe}}}\right)^2 \\ + 0 &= -P_{\mathrm{mech}}^{\mathrm{pmech}} + + \left(A_tH(q - q_{\mathrm{NL}}) - D_{\mathrm{turb}}\omega g\right) + \dfrac{S_{\mathrm{hygov}}^{\mathrm{base}}}{S^{\mathrm{pmech}}} +\end{aligned} +``` + +The $P_{\mathrm{GV}}^{\mathrm{safe}}$ lower bound protects the head divider near a closed gate. A physically consistent operating point satisfies $P_{\mathrm{GV}}>0.01$ so the lower bound is inactive. + +CommonMath defines the helper targets and smooth approximations for [max, clamp, deadband1, and linseg](../../../../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. If optional signals are not connected, use: + +```math +\begin{aligned} + \omega_0 &= 0 \\ + P_{\mathrm{aux},0} &= 0 +\end{aligned} +``` + +The connected power-source model initializes mechanical power first; HYGOV reads +the `pmech` signal $P_{\mathrm{mech},0}^{\mathrm{pmech}}$ and converts it to +the component-base value $P_{\mathrm{mech},0}^{\mathrm{hygov}}$. At synchronous +speed the damping term vanishes: + +```math +\begin{aligned} + P_{\mathrm{mech},0}^{\mathrm{hygov}} + &= P_{\mathrm{mech},0}^{\mathrm{pmech}} + \dfrac{S^{\mathrm{pmech}}}{S_{\mathrm{hygov}}^{\mathrm{base}}} \\ + H_0 &= H_{\mathrm{dam}} \\ + q_0 &= q_{\mathrm{NL}} + \dfrac{P_{\mathrm{mech},0}^{\mathrm{hygov}}}{A_tH_0} \\ + P_{\mathrm{GV},0} &= \dfrac{q_0}{\sqrt{H_0}} \\ + P_{\mathrm{GV},0}^{\mathrm{safe}} &= \text{max}\!\left(P_{\mathrm{GV},0}, 0.01\right) \\ + g_0 &= N_{\mathrm{GV}}^{-1}\!\left(P_{\mathrm{GV},0}\right) \\ + c_0 &= g_0 +\end{aligned} +``` + +Then evaluate the governor chain. With $\dot c_0=0$ and $\dot x_{f,0}=0$, the temporary-droop block requires $x_{f,0}=0$: + +```math +\begin{aligned} + \omega_{\mathrm{db},0} &= \text{deadband1}\!\left(\omega_0, -D_{\omega}, D_{\omega}\right) \\ + x_{n,0} &= \omega_{\mathrm{db},0} \\ + y_{\omega,0} &= \omega_{\mathrm{db},0} \\ + x_{f,0} &= 0 \\ + e_{f,0} &= x_{f,0} \\ + f_{c,0} &= 0 \\ + r_{c,0} &= 0 \\ + P_{\mathrm{ref},0} &= e_{f,0} - P_{\mathrm{aux},0} + y_{\omega,0} + R_{\mathrm{perm}}c_0 +\end{aligned} +``` + +For an unsaturated zero-derivative start, require $G^{\min} \le c_0 \le G^{\max}$, $P_{\mathrm{GV},0}>0.01$, and a single-valued inverse $N_{\mathrm{GV}}^{-1}(P_{\mathrm{GV},0})$. Starts that bind a gate limit or cannot be served at the available head are outside these closed-form initialization formulas. + +## Model Outputs + +Output | Units | Description | Note +---------------|--------|-------------------------------------|------ +`pmech` | [p.u.] | Mechanical power output | `mva` signal; oriented as mechanical input to the machine +`leadlag` | [p.u.] | Speed lead-lag denominator state | Not circled in Fig. 1 +`filter` | [p.u.] | Governor error filter output | State 1 +`desiredgate` | [p.u.] | Desired-gate position | State 2 +`gate` | [p.u.] | Gate position | State 3 +`flow` | [p.u.] | Turbine flow | State 4 +`head` | [p.u.] | Turbine head | $(q/P_{\mathrm{GV}}^{\mathrm{safe}})^2$ +`pgv` | [p.u.] | Nonlinear gate-to-power curve output | diff --git a/GridKit/Model/PhasorDynamics/Governor/README.md b/GridKit/Model/PhasorDynamics/Governor/README.md index 8e6305aa8..0a107bcf7 100644 --- a/GridKit/Model/PhasorDynamics/Governor/README.md +++ b/GridKit/Model/PhasorDynamics/Governor/README.md @@ -9,3 +9,4 @@ A governor models the control system that regulates the output power of a machin There are a few standard Governor models - Turbine Governor (See [TGOV1](Tgov1/README.md)) +- Hydro Turbine Governor (See [HYGOV](HYGOV/README.md)) diff --git a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md index 32df36cdc..43c8db66a 100644 --- a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md +++ b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md @@ -144,17 +144,18 @@ 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` + `Regca` | WECC REGCA renewable generator/converter model | `bus`, `ipcmd`\*, `iqcmd`\*, `ibranchr`\*, `ibranchi`\*, `pbranch`\*, `qbranch`\* | `P0`, `Q0`, `mva`, `Tg`, `TM`, `Rqmax`, `Rqmin`, `Rpmax`, `sL`, `IL1`, `VL0`, `VL1`, `VA0`, `VA1`, `Vhvmax` | `ir`, `ii`, `p`, `q`, `vt`, `vm`, `ip`, `iq`, `iqextra`, `il`, `lp`, `up` `Tgov1 ` | the TGOV1 governor model | `pmech`, `speed` | `R`, `T1`, `T2`, `T3`, `Pvmax`, `Pvmin`, `Dt` | `none` + `Hygov` | the HYGOV hydro turbine-governor model | `pmech`, `speed`\*, `pref`\*, `paux`\* | `Trate`, `mva`, `Rperm`, `Rtemp`, `Tr`, `Tf`, `Tg`, `Velm`, `Gmax`, `Gmin`, `Tw`, `At`, `Dturb`, `Qnl`, `Tn`, `Tnp`, `db1`, `db2`, `Hdam`, `Gv0`-`Gv5`, `Pgv0`-`Pgv5` | `pmech`, `leadlag`, `filter`, `desiredgate`, `gate`, `flow`, `head`, `pgv` `Ieeet1` | the IEEET1 exciter model | `bus`, `speed`, `efd`, `vs`\* | `Tr`, `Ka`, `Ta`, `Ke`, `Te`, `Kf`, `Tf`, `Vrmin`, `Vrmax`, `E1`, `E2`, `Se1`, `Se2`, `Ispdlim` | `efd`, `ksat` `SexsPti` | the SEXS-PTI simplified exciter model | `bus`, `efd`, `vs`\* | `Ta`, `Tb`, `Te`, `K`, `Efdmax`, `Efdmin` | `efd` - `Ieeest` | the IEEEST stabilizer model | `input`, `output` | `A1`, `A2`, `A3`, `A4`, `A5`, `A6`, `T1`, `T2`, `T3`, `T4`, `T5`, `T6`, `Ks`, `Lsmin`, `Lsmax`, `Vcl`, `Vcu`, `Tdelay` | `vss` + `Ieeest` | the IEEEST stabilizer model | `input`, `output` | `A1`, `A2`, `A3`, `A4`, `A5`, `A6`, `T1`, `T2`, `T3`, `T4`, `T5`, `T6`, `Ks`, `Lsmin`, `Lsmax`, `Vcl`, `Vcu`, `Tdelay` | `vss` `BusFault` | simple impedance-based fault at a bus | `bus`, `status`\* | `state0`, `R`, `X` | `state`, `ir`, `ii` `BusToSignalAdapter` | signal adapter component for a bus | `bus`, `vr`, `vi`, `ir`, `ii` | | Ports marked with \* are optional and, if missing, will be assumed to be connected to a constant value. This list is subject to change. - For `Branch`, `tap` and `phase` are optional parameters. If omitted, `tap` defaults to `1.0` and `phase` defaults to `0.0` radians. Bus `bus1` is the tap side for off-nominal transformer branches. diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index 25478738c..f8c087faf 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -82,6 +82,7 @@ namespace GridKit using namespace Governor; using namespace Exciter; using namespace Stabilizer; + using namespace Converter; // Set system model tolerances rel_tol_ = 1e-7; @@ -151,6 +152,64 @@ namespace GridKit addComponent(adapter); } + // Add REGCA converters + for (const auto& regcadata : data.regca) + { + using DataT = typename SystemModelData::RegcaDataT; + + IdxT bus_index = 0; + if (regcadata.ports.contains(DataT::Ports::bus)) + { + bus_index = regcadata.ports.at(DataT::Ports::bus); + } + + auto* regca = new Converter::Regca(getBus(bus_index), regcadata); + + if (regcadata.ports.contains(DataT::Ports::ipcmd)) + { + const IdxT ipcmd = regcadata.ports.at(DataT::Ports::ipcmd); + regca->getSignals().template attachSignalNode( + getSignal(ipcmd)); + } + + if (regcadata.ports.contains(DataT::Ports::iqcmd)) + { + const IdxT iqcmd = regcadata.ports.at(DataT::Ports::iqcmd); + regca->getSignals().template attachSignalNode( + getSignal(iqcmd)); + } + + if (regcadata.ports.contains(DataT::Ports::ibranchr)) + { + const IdxT ibranchr = regcadata.ports.at(DataT::Ports::ibranchr); + regca->getSignals().template assignSignalNode( + getSignal(ibranchr)); + } + + if (regcadata.ports.contains(DataT::Ports::ibranchi)) + { + const IdxT ibranchi = regcadata.ports.at(DataT::Ports::ibranchi); + regca->getSignals().template assignSignalNode( + getSignal(ibranchi)); + } + + if (regcadata.ports.contains(DataT::Ports::pbranch)) + { + const IdxT pbranch = regcadata.ports.at(DataT::Ports::pbranch); + regca->getSignals().template assignSignalNode( + getSignal(pbranch)); + } + + if (regcadata.ports.contains(DataT::Ports::qbranch)) + { + const IdxT qbranch = regcadata.ports.at(DataT::Ports::qbranch); + regca->getSignals().template assignSignalNode( + getSignal(qbranch)); + } + + addComponent(regca); + } + // Add branches for (const auto& branchdata : data.branch) { @@ -295,6 +354,38 @@ namespace GridKit addComponent(gov); } + // Add HYGOV governors + for (const auto& govdata : data.hygov) + { + auto* gov = new Hygov(govdata); + + if (govdata.ports.contains(HygovData::Ports::speed)) + { + IdxT speed = govdata.ports.at(HygovData::Ports::speed); + gov->getSignals().template attachSignalNode(getSignal(speed)); + } + + if (govdata.ports.contains(HygovData::Ports::pmech)) + { + IdxT pmech = govdata.ports.at(HygovData::Ports::pmech); + gov->getSignals().template assignSignalNode(getSignal(pmech)); + } + + if (govdata.ports.contains(HygovData::Ports::pref)) + { + IdxT pref = govdata.ports.at(HygovData::Ports::pref); + gov->getSignals().template attachSignalNode(getSignal(pref)); + } + + if (govdata.ports.contains(HygovData::Ports::paux)) + { + IdxT paux = govdata.ports.at(HygovData::Ports::paux); + gov->getSignals().template attachSignalNode(getSignal(paux)); + } + + addComponent(gov); + } + for (const auto& excitedata : data.exciter) { IdxT bus_index = 0; diff --git a/GridKit/Model/PhasorDynamics/SystemModelData.hpp b/GridKit/Model/PhasorDynamics/SystemModelData.hpp index 06862516c..8943e1f09 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelData.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelData.hpp @@ -10,8 +10,10 @@ #include #include #include +#include #include #include +#include #include #include #include @@ -38,7 +40,9 @@ namespace GridKit using BusDataT = BusData; using BusToSignalAdapterDataT = BusToSignalAdapterData; using BusFaultDataT = BusFaultData; + using RegcaDataT = Converter::RegcaData; using Tgov1DataT = Governor::Tgov1Data; + using HygovDataT = Governor::HygovData; using Ieeet1DataT = Exciter::Ieeet1Data; using SexsPtiDataT = Exciter::SexsPtiData; using IeeestDataT = Stabilizer::IeeestData; @@ -91,12 +95,14 @@ namespace GridKit std::vector adapter; ///< bus-to-signal adapters within the model std::vector branch; ///< Branches within the model std::vector bus_fault; ///< Bus faults within the model + std::vector regca; ///< REGCA converter instances within the model std::vector genrou; ///< GENROU instances within the model std::vector gensal; ///< GENSAL instances within the model std::vector genclassical; ///< Classical generator instances within the model std::vector load; ///< Loads within the model std::vector loadzip; ///< Loads within the model std::vector gov; ///< Governors within the model + std::vector hygov; ///< HYGOV governors within the model std::vector exciter; ///< Exciters within the model std::vector sexspti; ///< SEXS-PTI exciters within the model std::vector stabilizer; ///< Stabilizers within the model diff --git a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp index 00a6278c2..8b0a95cda 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp @@ -136,12 +136,24 @@ namespace GridKit raw_component.get_to(loadzip); sm.loadzip.push_back(loadzip); } + else if (kind == "Regca") + { + typename SystemModelData::RegcaDataT regca; + raw_component.get_to(regca); + sm.regca.push_back(regca); + } else if (kind == "Tgov1") { typename SystemModelData::Tgov1DataT gov; raw_component.get_to(gov); sm.gov.push_back(gov); } + else if (kind == "Hygov") + { + typename SystemModelData::HygovDataT gov; + raw_component.get_to(gov); + sm.hygov.push_back(gov); + } else if (kind == "Ieeet1") { typename SystemModelData::Ieeet1DataT exciter; diff --git a/docs/Figures/PhasorDynamics/HYGOV_diagram.png b/docs/Figures/PhasorDynamics/HYGOV_diagram.png new file mode 100755 index 000000000..fc7d7916d Binary files /dev/null and b/docs/Figures/PhasorDynamics/HYGOV_diagram.png differ diff --git a/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp b/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp index 41f6d2ca1..0f7f659b3 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..74783aca7 100644 --- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt +++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt @@ -78,6 +78,14 @@ target_link_libraries( GridKit::phasor_dynamics_bus_dependency_tracking GridKit::testing) +add_executable(test_phasor_governor_hygov runGovernorHygovTests.cpp) +target_link_libraries( + test_phasor_governor_hygov + GridKit::definitions + GridKit::phasor_dynamics_components + GridKit::phasor_dynamics_components_dependency_tracking + GridKit::testing) + add_executable(test_phasor_exciter_ieeet1 runExciterIeeet1Tests.cpp) target_link_libraries( test_phasor_exciter_ieeet1 @@ -94,6 +102,14 @@ target_link_libraries( GridKit::phasor_dynamics_components_dependency_tracking GridKit::testing) +add_executable(test_phasor_converter_regca runConverterRegcaTests.cpp) +target_link_libraries( + test_phasor_converter_regca + 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 @@ -135,9 +151,11 @@ add_test(NAME PhasorDynamicsBusToSignalAdapterTest COMMAND test_phasor_bustosign add_test(NAME PhasorDynamicsBranchTest COMMAND test_phasor_branch) add_test(NAME PhasorDynamicsGenrouTest COMMAND test_phasor_genrou) add_test(NAME PhasorDynamicsGovernorTgov1Test COMMAND test_phasor_governor_tgov1) +add_test(NAME PhasorDynamicsGovernorHygovTest COMMAND test_phasor_governor_hygov) add_test(NAME PhasorDynamicsExciterIeeet1Test COMMAND test_phasor_exciter_ieeet1) add_test(NAME PhasorDynamicsGensalTest COMMAND test_phasor_gensal) add_test(NAME PhasorDynamicsExciterSexsPtiTest COMMAND test_phasor_exciter_sexspti) +add_test(NAME PhasorDynamicsConverterRegcaTest COMMAND test_phasor_converter_regca) add_test(NAME PhasorDynamicsStabilizerIeeestTest COMMAND test_phasor_stabilizer_ieeest) add_test(NAME PhasorDynamicsGenClassicalTest COMMAND test_phasor_gen_classical) add_test(NAME PhasorDynamicsLoadTest COMMAND test_phasor_load) @@ -156,9 +174,11 @@ install( test_phasor_loadzip test_phasor_genrou test_phasor_governor_tgov1 + test_phasor_governor_hygov test_phasor_exciter_ieeet1 test_phasor_gensal test_phasor_exciter_sexspti + test_phasor_converter_regca test_phasor_stabilizer_ieeest test_phasor_gen_classical test_phasor_system diff --git a/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp new file mode 100644 index 000000000..90c77bd03 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp @@ -0,0 +1,932 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + template + class ConverterRegcaTests + { + public: + using RealT = typename PhasorDynamics::Component::RealT; + + ConverterRegcaTests() = default; + ~ConverterRegcaTests() = default; + + static constexpr ScalarT kTol = static_cast(1.0e-9); + + TestOutcome constructionAndValidation() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + + PhasorDynamics::Converter::Regca minimal(&bus); + success *= (minimal.size() == static_cast(Vars::MAXIMUM)); + success *= (minimal.getMonitor() == nullptr); + success *= (minimal.verify() > 0); + + PhasorDynamics::Converter::Regca configured(&bus, makeData()); + success *= (configured.size() == static_cast(Vars::MAXIMUM)); + success *= (configured.getMonitor() != nullptr); + success *= (configured.verify() == 0); + + return success.report(__func__); + } + + TestOutcome parameterValidation() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + + auto missing = makeData(); + missing.parameters.erase(Params::Tg); + PhasorDynamics::Converter::Regca missing_model(&bus, missing); + success *= (missing_model.verify() > 0); + + auto missing_p0 = makeData(); + missing_p0.parameters.erase(Params::P0); + PhasorDynamics::Converter::Regca missing_p0_model(&bus, missing_p0); + success *= (missing_p0_model.verify() > 0); + + auto missing_q0 = makeData(); + missing_q0.parameters.erase(Params::Q0); + PhasorDynamics::Converter::Regca missing_q0_model(&bus, missing_q0); + success *= (missing_q0_model.verify() > 0); + + auto bad_switch = makeData(); + bad_switch.parameters[Params::sL] = static_cast(2); + PhasorDynamics::Converter::Regca bad_switch_model(&bus, bad_switch); + success *= (bad_switch_model.verify() > 0); + + success *= invalidParameterCase(bus, Params::mva, static_cast(0.0)); + success *= invalidParameterCase(bus, Params::Tg, static_cast(0.0)); + success *= invalidParameterCase(bus, Params::TM, static_cast(0.0)); + success *= invalidParameterCase(bus, Params::Rpmax, static_cast(0.0)); + success *= invalidParameterCase(bus, Params::Rqmin, static_cast(0.0)); + success *= invalidParameterCase(bus, Params::IL1, static_cast(-0.1)); + success *= invalidParameterCase(bus, Params::VL1, static_cast(0.3)); + success *= invalidParameterCase(bus, Params::VA1, static_cast(0.3)); + success *= invalidParameterCase(bus, Params::Vhvmax, static_cast(0.0)); + + return success.report(__func__); + } + + TestOutcome initializesFromPowerFlowAndPublishesSignals() + { + TestStatus success = true; + + const ScalarT vr{0.8}; + const ScalarT vi{0.6}; + const ScalarT p0{0.4}; + const ScalarT q0{-0.1}; + const RealT mva{50.0}; + + PhasorDynamics::Bus bus(vr, vi); + bus.allocate(); + bus.initialize(); + + auto data = makeData(); + data.parameters[Params::mva] = mva; + data.parameters[Params::P0] = static_cast(p0); + data.parameters[Params::Q0] = static_cast(q0); + + PhasorDynamics::Converter::Regca regca(&bus, data); + + ScalarT ipcmd_value{-1.0}; + ScalarT iqcmd_value{1.0}; + IdxT ipcmd_index = 20; + IdxT iqcmd_index = 21; + + PhasorDynamics::SignalNode ipcmd_node; + PhasorDynamics::SignalNode iqcmd_node; + PhasorDynamics::SignalNode ir_node; + PhasorDynamics::SignalNode ii_node; + PhasorDynamics::SignalNode pbranch_node; + PhasorDynamics::SignalNode qbranch_node; + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + + regca.getSignals().template attachSignalNode(&ipcmd_node); + regca.getSignals().template attachSignalNode(&iqcmd_node); + regca.getSignals().template assignSignalNode(&ir_node); + regca.getSignals().template assignSignalNode(&ii_node); + regca.getSignals().template assignSignalNode(&pbranch_node); + regca.getSignals().template assignSignalNode(&qbranch_node); + + success *= (regca.allocate() == 0); + success *= (regca.verify() == 0); + success *= (regca.initialize() == 0); + success *= (regca.evaluateResidual() == 0); + + const ScalarT vt = std::sqrt(vr * vr + vi * vi); + const ScalarT lvacm = Math::linseg(vt, static_cast(0.4), static_cast(0.9), ONE); + const ScalarT ipcmd0 = toComponentBase(p0 / vt, mva) / lvacm; + const ScalarT iqcmd0 = toComponentBase(q0 / vt, mva); + const ScalarT qnet0 = iqcmd0; + const ScalarT ir0 = (lvacm * ipcmd0 * vr + qnet0 * vi) / vt; + const ScalarT ii0 = (lvacm * ipcmd0 * vi - qnet0 * vr) / vt; + + success *= scalarMatches(regca.y()[index(Vars::VM)], vt, "VM"); + success *= scalarMatches(regca.y()[index(Vars::VT)], vt, "VT"); + success *= scalarMatches(regca.y()[index(Vars::IP)], ipcmd0, "IP"); + success *= scalarMatches(regca.y()[index(Vars::IQ)], iqcmd0, "IQ"); + success *= scalarMatches(regca.y()[index(Vars::IR)], ir0, "IR"); + success *= scalarMatches(regca.y()[index(Vars::II)], ii0, "II"); + success *= scalarMatches(regca.y()[index(Vars::PBR)], p0, "PBR"); + success *= scalarMatches(regca.y()[index(Vars::QBR)], q0, "QBR"); + success *= scalarMatches(ipcmd_node.read(), ipcmd0, "ipcmd signal"); + success *= scalarMatches(iqcmd_node.read(), iqcmd0, "iqcmd signal"); + success *= scalarMatches(ir_node.read(), ir0, "ir signal"); + success *= scalarMatches(ii_node.read(), ii0, "ii signal"); + success *= scalarMatches(pbranch_node.read(), p0, "pbranch signal"); + success *= scalarMatches(qbranch_node.read(), q0, "qbranch signal"); + + success *= allResidualsZero(regca); + return success.report(__func__); + } + + TestOutcome unconnectedCommandsRemainConstant() + { + TestStatus success = true; + + auto data = makeData(); + data.parameters[Params::P0] = static_cast(0.6); + data.parameters[Params::Q0] = static_cast(0.2); + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + PhasorDynamics::Converter::Regca regca(&bus, data); + + success *= (regca.allocate() == 0); + success *= (regca.verify() == 0); + success *= (regca.initialize() == 0); + success *= (regca.evaluateResidual() == 0); + + success *= isEqual(regca.y()[index(Vars::IP)], static_cast(0.6), kTol); + success *= isEqual(regca.y()[index(Vars::IQ)], static_cast(0.2), kTol); + success *= allResidualsZero(regca); + + return success.report(__func__); + } + + TestOutcome externalCommandsDriveRuntimeResidual() + { + TestStatus success = true; + + auto data = makeData(); + data.parameters[Params::P0] = static_cast(0.5); + data.parameters[Params::Q0] = static_cast(-0.1); + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + ScalarT ipcmd_value{0.0}; + ScalarT iqcmd_value{0.0}; + IdxT ipcmd_index = 22; + IdxT iqcmd_index = 23; + + PhasorDynamics::SignalNode ipcmd_node; + PhasorDynamics::SignalNode iqcmd_node; + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + + PhasorDynamics::Converter::Regca regca(&bus, data); + regca.getSignals().template attachSignalNode(&ipcmd_node); + regca.getSignals().template attachSignalNode(&iqcmd_node); + + success *= (regca.allocate() == 0); + success *= (regca.verify() == 0); + success *= (regca.initialize() == 0); + + ipcmd_value = static_cast(0.55); + iqcmd_value = static_cast(-0.05); + bus.evaluateResidual(); + success *= (regca.evaluateResidual() == 0); + + const ScalarT tg = static_cast(0.02); + const ScalarT ip_error = ipcmd_value - regca.y()[index(Vars::IP)]; + const ScalarT iq_error = iqcmd_value - regca.y()[index(Vars::IQ)]; + + const ScalarT expected_ip_residual = + tg * regca.y()[index(Vars::LP)] + + Math::ramp(ip_error - tg * regca.y()[index(Vars::LP)]) + - Math::ramp(ip_error - tg * regca.y()[index(Vars::UP)]); + const ScalarT expected_iq_residual = + iq_error + Math::ramp(tg * static_cast(-999.0) - iq_error); + + success *= isEqual(regca.getResidual()[index(Vars::IP)], expected_ip_residual, kTol); + success *= isEqual(regca.getResidual()[index(Vars::IQ)], expected_iq_residual, kTol); + success *= (regca.getResidual()[index(Vars::IP)] > ZERO); + success *= (regca.getResidual()[index(Vars::IQ)] > ZERO); + + return success.report(__func__); + } + + TestOutcome invalidInitialization() + { + TestStatus success = true; + + auto data = makeData(); + + { + PhasorDynamics::Bus bus(1.3, 0.0); + bus.allocate(); + bus.initialize(); + + PhasorDynamics::Converter::Regca regca(&bus, data); + success *= (regca.allocate() == 0); + success *= (regca.initialize() > 0); + } + + { + PhasorDynamics::Bus bus(0.0, 0.0); + bus.allocate(); + bus.initialize(); + + PhasorDynamics::Converter::Regca regca(&bus, data); + success *= (regca.allocate() == 0); + success *= (regca.initialize() > 0); + } + + { + auto low_voltage_data = data; + low_voltage_data.parameters[Params::P0] = static_cast(0.5); + + PhasorDynamics::Bus bus(0.2, 0.0); + bus.allocate(); + bus.initialize(); + + PhasorDynamics::Converter::Regca regca(&bus, low_voltage_data); + success *= (regca.allocate() == 0); + success *= (regca.initialize() > 0); + } + + { + auto low_voltage_data = data; + low_voltage_data.parameters[Params::P0] = static_cast(0.0); + low_voltage_data.parameters[Params::Q0] = static_cast(0.1); + + PhasorDynamics::Bus bus(0.2, 0.0); + bus.allocate(); + bus.initialize(); + + PhasorDynamics::Converter::Regca regca(&bus, low_voltage_data); + success *= (regca.allocate() == 0); + success *= (regca.initialize() == 0); + } + + return success.report(__func__); + } + + TestOutcome signalVerification() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + PhasorDynamics::Converter::Regca regca(&bus, makeData()); + + PhasorDynamics::SignalNode ipcmd_node; + PhasorDynamics::SignalNode iqcmd_node; + ScalarT ipcmd_value{0.25}; + ScalarT iqcmd_value{-0.10}; + IdxT ipcmd_index = 24; + IdxT iqcmd_index = 25; + + regca.getSignals().template attachSignalNode(&ipcmd_node); + success *= (regca.verify() > 0); + + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + success *= (regca.verify() == 0); + + regca.getSignals().template attachSignalNode(&iqcmd_node); + success *= (regca.verify() > 0); + + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + success *= (regca.verify() == 0); + + return success.report(__func__); + } + + TestOutcome nullBusVerification() + { + TestStatus success = true; + + PhasorDynamics::Converter::Regca regca(nullptr, makeData()); + success *= (regca.verify() > 0); + success *= (regca.initialize() > 0); + + return success.report(__func__); + } + + TestOutcome busInjectionUsesSystemBase() + { + TestStatus success = true; + + const ScalarT vr{0.8}; + const ScalarT vi{0.6}; + const ScalarT p0{0.4}; + const ScalarT q0{-0.1}; + const RealT mva{50.0}; + + PhasorDynamics::Bus bus(vr, vi); + bus.allocate(); + bus.initialize(); + + auto data = makeData(); + data.parameters[Params::mva] = mva; + data.parameters[Params::P0] = static_cast(p0); + data.parameters[Params::Q0] = static_cast(q0); + + PhasorDynamics::Converter::Regca regca(&bus, data); + success *= (regca.allocate() == 0); + success *= (regca.initialize() == 0); + + bus.evaluateResidual(); + success *= (regca.evaluateResidual() == 0); + + success *= isEqual(bus.Ir(), toSystemBase(regca.y()[index(Vars::IR)], mva), kTol); + success *= isEqual(bus.Ii(), toSystemBase(regca.y()[index(Vars::II)], mva), kTol); + success *= isEqual(bus.Ir(), static_cast(0.26), kTol); + success *= isEqual(bus.Ii(), static_cast(0.32), kTol); + + return success.report(__func__); + } + + TestOutcome residualEquations() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(0.95, 0.25); + bus.allocate(); + bus.initialize(); + + ScalarT ipcmd_value{0.9}; + ScalarT iqcmd_value{0.1}; + IdxT ipcmd_index = 26; + IdxT iqcmd_index = 27; + + PhasorDynamics::SignalNode ipcmd_node; + PhasorDynamics::SignalNode iqcmd_node; + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + + PhasorDynamics::Converter::Regca regca(&bus, makeDynamicData()); + regca.getSignals().template attachSignalNode(&ipcmd_node); + regca.getSignals().template attachSignalNode(&iqcmd_node); + + success *= (regca.allocate() == 0); + success *= (regca.initialize() == 0); + + setResidualState(regca); + bus.evaluateResidual(); + success *= (regca.evaluateResidual() == 0); + + const auto expected = expectedResidualForState(ipcmd_value, iqcmd_value); + success *= vectorMatches(regca.getResidual(), expected, "REGCA residual"); + + return success.report(__func__); + } + + TestOutcome highVoltageReactiveCurrentRoot() + { + TestStatus success = true; + + auto data = makeDynamicData(); + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + PhasorDynamics::Converter::Regca regca(&bus, data); + success *= (regca.allocate() == 0); + success *= (regca.initialize() == 0); + + regca.y()[index(Vars::VT)] = static_cast(1.35); + regca.y()[index(Vars::IQEXTRA)] = + Math::ramp(regca.y()[index(Vars::VT)] - static_cast(1.3)); + + bus.evaluateResidual(); + success *= (regca.evaluateResidual() == 0); + success *= scalarMatches(regca.getResidual()[index(Vars::IQEXTRA)], + static_cast(0.0), + "HVRCM residual root"); + + return success.report(__func__); + } + + TestOutcome limiterBranchCoverage() + { + TestStatus success = true; + + { + auto data = makeDynamicData(); + data.parameters[Params::Q0] = static_cast(0.2); + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + ScalarT iqcmd_value{0.2}; + IdxT iqcmd_index = 28; + PhasorDynamics::SignalNode iqcmd_node; + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + + PhasorDynamics::Converter::Regca regca(&bus, data); + regca.getSignals().template attachSignalNode(&iqcmd_node); + + success *= (regca.allocate() == 0); + success *= (regca.verify() == 0); + success *= (regca.initialize() == 0); + + iqcmd_value = static_cast(0.4); + bus.evaluateResidual(); + success *= (regca.evaluateResidual() == 0); + + const ScalarT iq_error = iqcmd_value - regca.y()[index(Vars::IQ)]; + const ScalarT expected = + iq_error - Math::ramp(iq_error - static_cast(0.2) * static_cast(0.5)); + + success *= scalarMatches(regca.getResidual()[index(Vars::IQ)], + expected, + "positive IQ branch"); + } + + { + auto data = makeDynamicData(); + data.parameters[Params::sL] = false; + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + PhasorDynamics::Converter::Regca regca(&bus, data); + success *= (regca.allocate() == 0); + success *= (regca.initialize() == 0); + + bus.evaluateResidual(); + success *= (regca.evaluateResidual() == 0); + + const ScalarT ip = regca.y()[index(Vars::IP)]; + const ScalarT sigma_ip = Math::sigmoid(ip); + const RealT rpmax = static_cast(0.7); + const RealT mp = static_cast(100.0) * rpmax; + const ScalarT expected = mp * (ONE - sigma_ip) + rpmax * sigma_ip; + + success *= scalarMatches(regca.y()[index(Vars::UP)], expected, "UP without LVPL"); + success *= scalarMatches(regca.getResidual()[index(Vars::UP)], + static_cast(0.0), + "UP residual without LVPL"); + } + + return success.report(__func__); + } + + TestOutcome jsonParseAndSystemAssembly() + { + TestStatus success = true; + + std::istringstream input(R"json( +{ + "header": { + "format_version": 0, + "format_revision": 1, + "case_name": "REGCA full model", + "case_description": "REGCA parser behavior test", + "case_comments": "", + "freq_base": 60.0, + "va_base": 100000000.0 + }, + "buses": [ + { "number": 1, "class": "bus", "name": "Bus 1", "init": { "Vr": 1.0, "Vi": 0.0 }, "v_base": 1.0 } + ], + "devices": [ + { + "class": "Regca", + "ports": { "bus": 1 }, + "id": "CV1", + "params": { + "P0": 0.0, + "Q0": 0.0, + "mva": 100, + "Tg": 0.02, + "TM": 0.02, + "Rqmax": 999.0, + "Rqmin": -999.0, + "Rpmax": 999.0, + "sL": true, + "IL1": 1.1, + "VL0": 0.4, + "VL1": 0.9, + "VA0": 0.4, + "VA1": 0.9, + "Vhvmax": 1.2 + }, + "mon": ["ir", "ii", "p", "q"] + } + ] +} +)json"); + + auto data = PhasorDynamics::parseSystemModelData(input); + success *= (data.regca.size() == 1); + const auto& regca_data = data.regca[0]; + success *= (regca_data.device_class == "Regca"); + success *= (regca_data.ports.at(PhasorDynamics::Converter::RegcaPorts::bus) == 1); + success *= (std::get_if(®ca_data.parameters.at(Params::P0)) != nullptr); + success *= (std::get_if(®ca_data.parameters.at(Params::Q0)) != nullptr); + success *= (std::get_if(®ca_data.parameters.at(Params::mva)) != nullptr); + success *= (std::get_if(®ca_data.parameters.at(Params::sL)) != nullptr); + + PhasorDynamics::SystemModel system(data); + success *= (system.allocate() == 0); + success *= (system.initialize() == 0); + success *= (system.tagDifferentiable() == 0); + success *= (system.evaluateResidual() == 0); + success *= (system.evaluateJacobian() == 0); + success *= (system.size() == 14); + success *= isEqual(system.getResidual()[0], 0.0, static_cast(kTol)); + success *= isEqual(system.getResidual()[1], 0.0, static_cast(kTol)); + for (size_t i = 2; i < system.getResidual().size(); ++i) + { + success *= isEqual(system.getResidual()[i], 0.0, static_cast(kTol)); + } + + return success.report(__func__); + } + +#ifdef GRIDKIT_ENABLE_ENZYME + TestOutcome jacobian() + { + TestStatus success = true; + + auto dependency_tracking_jacobian = dependencyTrackingJacobian(); + auto enzyme_jacobian = enzymeJacobian(); + + success *= (dependency_tracking_jacobian.size() == enzyme_jacobian.size()); + const auto nrows = std::min(dependency_tracking_jacobian.size(), enzyme_jacobian.size()); + for (size_t i = 0; i < nrows; ++i) + { + success *= isEqual(dependency_tracking_jacobian[i], enzyme_jacobian[i], static_cast(1.0e-8)); + } + + return success.report(__func__); + } +#endif + + private: + using Params = PhasorDynamics::Converter::RegcaParameters; + using Vars = PhasorDynamics::Converter::RegcaInternalVariables; + using Ext = PhasorDynamics::Converter::RegcaExternalVariables; + + static size_t index(Vars variable) + { + return static_cast(variable); + } + + auto makeData() -> PhasorDynamics::Converter::RegcaData + { + using Ports = PhasorDynamics::Converter::RegcaPorts; + using Mon = PhasorDynamics::Converter::RegcaMonitorableVariables; + + PhasorDynamics::Converter::RegcaData data; + data.device_class = "Regca"; + data.disambiguation_string = "regca_test"; + data.ports[Ports::bus] = 1; + data.monitored_variables.insert(Mon::ir); + data.monitored_variables.insert(Mon::ii); + data.monitored_variables.insert(Mon::p); + data.monitored_variables.insert(Mon::q); + + data.parameters[Params::P0] = static_cast(0.0); + data.parameters[Params::Q0] = static_cast(0.0); + data.parameters[Params::mva] = static_cast(100.0); + data.parameters[Params::Tg] = static_cast(0.02); + data.parameters[Params::TM] = static_cast(0.02); + data.parameters[Params::Rqmax] = static_cast(999.0); + data.parameters[Params::Rqmin] = static_cast(-999.0); + data.parameters[Params::Rpmax] = static_cast(999.0); + data.parameters[Params::sL] = true; + data.parameters[Params::IL1] = static_cast(1.1); + data.parameters[Params::VL0] = static_cast(0.4); + data.parameters[Params::VL1] = static_cast(0.9); + data.parameters[Params::VA0] = static_cast(0.4); + data.parameters[Params::VA1] = static_cast(0.9); + data.parameters[Params::Vhvmax] = static_cast(1.2); + + return data; + } + + auto makeDynamicData() -> PhasorDynamics::Converter::RegcaData + { + auto data = makeData(); + data.parameters[Params::P0] = static_cast(0.6); + data.parameters[Params::Q0] = static_cast(-0.2); + data.parameters[Params::Tg] = static_cast(0.2); + data.parameters[Params::TM] = static_cast(0.4); + data.parameters[Params::Rqmax] = static_cast(0.5); + data.parameters[Params::Rqmin] = static_cast(-0.6); + data.parameters[Params::Rpmax] = static_cast(0.7); + data.parameters[Params::IL1] = static_cast(1.1); + data.parameters[Params::Vhvmax] = static_cast(1.3); + return data; + } + + bool invalidParameterCase(PhasorDynamics::Bus& bus, Params param, RealT value) + { + auto data = makeData(); + data.parameters[param] = value; + PhasorDynamics::Converter::Regca model(&bus, data); + return model.verify() > 0; + } + + ScalarT toComponentBase(ScalarT value, RealT mva) const + { + return value * static_cast(100.0) / mva; + } + + ScalarT toSystemBase(ScalarT value, RealT mva) const + { + return value * mva / static_cast(100.0); + } + + ScalarT activeLower(ScalarT ip, RealT rpmax) const + { + const RealT mp = static_cast(100.0) * rpmax; + return -rpmax - (mp - rpmax) * Math::sigmoid(ip); + } + + ScalarT activeUpper(ScalarT ip, ScalarT il, RealT rpmax, bool use_lvpl) const + { + const RealT mp = static_cast(100.0) * rpmax; + const ScalarT sigma_ip = Math::sigmoid(ip); + ScalarT output = mp * (ONE - sigma_ip) + rpmax * sigma_ip; + if (use_lvpl) + { + output = mp * (ONE - sigma_ip) + rpmax * sigma_ip * Math::sigmoid(il - ip); + } + return output; + } + + bool allResidualsZero(PhasorDynamics::Converter::Regca& regca) const + { + bool success = true; + for (size_t i = 0; i < static_cast(regca.size()); ++i) + { + if (!isEqual(regca.getResidual()[i], static_cast(0.0), kTol)) + { + std::cout << "REGCA residual row " << i << " is " << regca.getResidual()[i] << "\n"; + success = false; + } + if (!isEqual(regca.yp()[i], static_cast(0.0), kTol)) + { + std::cout << "REGCA derivative row " << i << " is " << regca.yp()[i] << "\n"; + success = false; + } + } + return success; + } + + bool vectorMatches(const std::vector& actual, + const std::vector& expected, + const char* label) const + { + bool success = (actual.size() == expected.size()); + const auto n = std::min(actual.size(), expected.size()); + for (size_t i = 0; i < n; ++i) + { + if (!isEqual(actual[i], expected[i], kTol)) + { + std::cout << label << " mismatch at row " << i << ": " + << actual[i] << " != " << expected[i] << "\n"; + success = false; + } + } + return success; + } + + bool scalarMatches(ScalarT actual, ScalarT expected, const char* label) const + { + if (isEqual(actual, expected, kTol)) + { + return true; + } + + std::cout << label << " mismatch: " << actual << " != " << expected << "\n"; + return false; + } + + void setResidualState(PhasorDynamics::Converter::Regca& regca) + { + regca.y()[index(Vars::VM)] = static_cast(0.86); + regca.y()[index(Vars::IQ)] = static_cast(-0.2); + regca.y()[index(Vars::IP)] = static_cast(0.85); + regca.y()[index(Vars::VT)] = static_cast(0.98); + regca.y()[index(Vars::II)] = static_cast(0.18); + regca.y()[index(Vars::IQEXTRA)] = static_cast(0.03); + regca.y()[index(Vars::IL)] = static_cast(0.72); + regca.y()[index(Vars::IR)] = static_cast(0.5); + regca.y()[index(Vars::LP)] = static_cast(-0.4); + regca.y()[index(Vars::UP)] = static_cast(0.3); + regca.y()[index(Vars::PBR)] = static_cast(0.52); + regca.y()[index(Vars::QBR)] = static_cast(-0.046); + + regca.yp()[index(Vars::VM)] = static_cast(0.01); + regca.yp()[index(Vars::IQ)] = static_cast(-0.02); + regca.yp()[index(Vars::IP)] = static_cast(0.03); + } + + std::vector expectedResidualForState(ScalarT ipcmd, ScalarT iqcmd) const + { + constexpr RealT tg = static_cast(0.2); + constexpr RealT tm = static_cast(0.4); + constexpr RealT rqmin = static_cast(-0.6); + constexpr RealT rpmax = static_cast(0.7); + constexpr RealT vl0 = static_cast(0.4); + constexpr RealT vl1 = static_cast(0.9); + constexpr RealT va0 = static_cast(0.4); + constexpr RealT va1 = static_cast(0.9); + constexpr RealT il1 = static_cast(1.1); + constexpr ScalarT vr = static_cast(0.95); + constexpr ScalarT vi = static_cast(0.25); + constexpr ScalarT vm = static_cast(0.86); + constexpr ScalarT iq = static_cast(-0.2); + constexpr ScalarT ip = static_cast(0.85); + constexpr ScalarT vt = static_cast(0.98); + constexpr ScalarT ii = static_cast(0.18); + constexpr ScalarT iqext = static_cast(0.03); + constexpr ScalarT il = static_cast(0.72); + constexpr ScalarT ir = static_cast(0.5); + constexpr ScalarT lp = static_cast(-0.4); + constexpr ScalarT up = static_cast(0.3); + constexpr ScalarT pbr = static_cast(0.52); + constexpr ScalarT qbr = static_cast(-0.046); + constexpr ScalarT vmdot = static_cast(0.01); + constexpr ScalarT iqdot = static_cast(-0.02); + constexpr ScalarT ipdot = static_cast(0.03); + + const ScalarT iq_error = iqcmd - iq; + const ScalarT ip_error = ipcmd - ip; + const ScalarT lvacm = Math::linseg(vt, va0, va1, ONE); + const ScalarT qnet = iq - iqext; + + std::vector expected(static_cast(Vars::MAXIMUM), static_cast(0.0)); + expected[index(Vars::VM)] = -tm * vmdot - vm + vt; + expected[index(Vars::IQ)] = -tg * iqdot + iq_error + Math::ramp(tg * rqmin - iq_error); + expected[index(Vars::IP)] = + -tg * ipdot + tg * lp + Math::ramp(ip_error - tg * lp) - Math::ramp(ip_error - tg * up); + expected[index(Vars::VT)] = -vt * vt + vr * vr + vi * vi; + expected[index(Vars::II)] = -vt * ii + lvacm * ip * vi - qnet * vr; + expected[index(Vars::IQEXTRA)] = -iqext + Math::ramp(vt - static_cast(1.3)); + expected[index(Vars::IL)] = -il + Math::linseg(vm, vl0, vl1, il1); + expected[index(Vars::IR)] = -vt * ir + lvacm * ip * vr + qnet * vi; + expected[index(Vars::LP)] = -lp + activeLower(ip, rpmax); + expected[index(Vars::UP)] = -up + activeUpper(ip, il, rpmax, true); + expected[index(Vars::PBR)] = -pbr + vr * ir + vi * ii; + expected[index(Vars::QBR)] = -qbr + vi * ir - vr * ii; + return expected; + } + +#ifdef GRIDKIT_ENABLE_ENZYME + void setResidualStateDep(PhasorDynamics::Converter::Regca& regca, + PhasorDynamics::Bus& bus) + { + bus.y()[0].setValue(0.95); + bus.y()[1].setValue(0.25); + + regca.y()[index(Vars::VM)].setValue(0.86); + regca.y()[index(Vars::IQ)].setValue(-0.2); + regca.y()[index(Vars::IP)].setValue(0.85); + regca.y()[index(Vars::VT)].setValue(0.98); + regca.y()[index(Vars::II)].setValue(0.18); + regca.y()[index(Vars::IQEXTRA)].setValue(0.03); + regca.y()[index(Vars::IL)].setValue(0.72); + regca.y()[index(Vars::IR)].setValue(0.5); + regca.y()[index(Vars::LP)].setValue(-0.4); + regca.y()[index(Vars::UP)].setValue(0.3); + regca.y()[index(Vars::PBR)].setValue(0.52); + regca.y()[index(Vars::QBR)].setValue(-0.046); + + regca.yp()[index(Vars::VM)].setValue(0.01); + regca.yp()[index(Vars::IQ)].setValue(-0.02); + regca.yp()[index(Vars::IP)].setValue(0.03); + } + + std::vector dependencyTrackingJacobian() + { + using DepVar = DependencyTracking::Variable; + + auto data = makeDynamicData(); + + PhasorDynamics::Bus bus(DepVar{0.95}, DepVar{0.25}); + PhasorDynamics::Converter::Regca regca(&bus, data); + + PhasorDynamics::SignalNode ipcmd_node; + PhasorDynamics::SignalNode iqcmd_node; + DepVar ipcmd_value{0.9}; + DepVar iqcmd_value{0.1}; + IdxT ipcmd_index = static_cast(regca.size() + bus.size()); + IdxT iqcmd_index = static_cast(regca.size() + bus.size() + 1); + + bus.allocate(); + regca.allocate(); + bus.initialize(); + setResidualStateDep(regca, bus); + + for (IdxT i = 0; i < regca.size(); ++i) + { + regca.y()[static_cast(i)].setVariableNumber(i); + regca.yp()[static_cast(i)].setVariableNumber(i); + } + for (IdxT i = 0; i < bus.size(); ++i) + { + bus.y()[static_cast(i)].setVariableNumber(i + regca.size()); + } + ipcmd_value.setVariableNumber(ipcmd_index); + iqcmd_value.setVariableNumber(iqcmd_index); + + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + regca.getSignals().template attachSignalNode(&ipcmd_node); + regca.getSignals().template attachSignalNode(&iqcmd_node); + + bus.evaluateResidual(); + regca.evaluateResidual(); + + std::vector dependencies( + static_cast(regca.size() + bus.size())); + for (IdxT i = 0; i < regca.size(); ++i) + { + dependencies[static_cast(i)] = regca.getResidual()[static_cast(i)].getDependencies(); + } + dependencies[static_cast(regca.size())] = bus.Ir().getDependencies(); + dependencies[static_cast(regca.size() + 1)] = bus.Ii().getDependencies(); + + return dependencies; + } + + std::vector enzymeJacobian() + { + auto data = makeDynamicData(); + + PhasorDynamics::Bus bus(0.95, 0.25); + PhasorDynamics::Converter::Regca regca(&bus, data); + + PhasorDynamics::SignalNode ipcmd_node; + PhasorDynamics::SignalNode iqcmd_node; + ScalarT ipcmd_value{0.9}; + ScalarT iqcmd_value{0.1}; + IdxT ipcmd_index = static_cast(regca.size() + bus.size()); + IdxT iqcmd_index = static_cast(regca.size() + bus.size() + 1); + + bus.allocate(); + regca.allocate(); + for (IdxT i = 0; i < bus.size(); ++i) + { + bus.setVariableIndex(i, i + regca.size()); + bus.setResidualIndex(i, i + regca.size()); + } + + bus.initialize(); + setResidualState(regca); + regca.updateTime(0.0, 1.0); + + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + regca.getSignals().template attachSignalNode(&ipcmd_node); + regca.getSignals().template attachSignalNode(&iqcmd_node); + + bus.evaluateResidual(); + regca.evaluateResidual(); + regca.evaluateJacobian(); + + auto model_jacobian = regca.getJacobian(); + model_jacobian.deduplicate(); + return MapFromCOO(model_jacobian); + } +#endif + }; + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/GovernorHygovTests.hpp b/tests/UnitTests/PhasorDynamics/GovernorHygovTests.hpp new file mode 100644 index 000000000..ef13a92f9 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/GovernorHygovTests.hpp @@ -0,0 +1,345 @@ +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + template + class GovernorHygovTests + { + public: + using RealT = typename PhasorDynamics::Component::RealT; + + static constexpr ScalarT kTol = static_cast(1.0e-8); + + TestOutcome constructionAndValidation() + { + TestStatus success = true; + + PhasorDynamics::Governor::Hygov hygov(makeHygovData()); + success *= (hygov.size() + == static_cast(PhasorDynamics::Governor::HygovInternalVariables::MAXIMUM)); + success *= (hygov.getMonitor() != nullptr); + success *= (hygov.verify() == 0); + + auto source_default_hygov = makeHygovSourceDefaultData(); + PhasorDynamics::Governor::Hygov source_default_hygov_model(source_default_hygov); + success *= (source_default_hygov_model.verify() == 0); + + auto bad_hygov = makeHygovData(); + bad_hygov.parameters[PhasorDynamics::Governor::HygovParameters::db2] = static_cast(0.01); + PhasorDynamics::Governor::Hygov bad_hygov_model(bad_hygov); + success *= (bad_hygov_model.verify() > 0); + + auto bad_hygov_leadlag = makeHygovData(); + bad_hygov_leadlag.parameters[PhasorDynamics::Governor::HygovParameters::Tn] = + static_cast(0.1); + bad_hygov_leadlag.parameters[PhasorDynamics::Governor::HygovParameters::Tnp] = + static_cast(0.0); + PhasorDynamics::Governor::Hygov bad_hygov_leadlag_model(bad_hygov_leadlag); + success *= (bad_hygov_leadlag_model.verify() > 0); + + auto bad_hygov_curve = makeHygovData(); + bad_hygov_curve.parameters[PhasorDynamics::Governor::HygovParameters::Gv2] = static_cast(0.2); + PhasorDynamics::Governor::Hygov bad_hygov_curve_model(bad_hygov_curve); + success *= (bad_hygov_curve_model.verify() > 0); + + auto bad_hygov_trate = makeHygovData(); + bad_hygov_trate.parameters[PhasorDynamics::Governor::HygovParameters::Trate] = + static_cast(0.0); + PhasorDynamics::Governor::Hygov bad_hygov_trate_model(bad_hygov_trate); + success *= (bad_hygov_trate_model.verify() > 0); + + auto bad_hygov_base = makeHygovData(); + bad_hygov_base.parameters[PhasorDynamics::Governor::HygovParameters::mva] = + static_cast(0.0); + PhasorDynamics::Governor::Hygov bad_hygov_base_model(bad_hygov_base); + success *= (bad_hygov_base_model.verify() > 0); + + return success.report(__func__); + } + + TestOutcome signals() + { + using Var = PhasorDynamics::Governor::HygovInternalVariables; + + TestStatus success = true; + + PhasorDynamics::SignalNode pmech_node; + PhasorDynamics::SignalNode omega_node; + ScalarT pmech_value{0.40}; + ScalarT omega_value{0.0}; + IdxT pmech_index = 4; + IdxT omega_index = 5; + pmech_node.set(&pmech_value, &pmech_index); + omega_node.set(&omega_value, &omega_index); + + PhasorDynamics::Governor::Hygov hygov(makeHygovData()); + auto& hygov_signals = hygov.getSignals(); + hygov_signals.template assignSignalNode( + &pmech_node); + hygov_signals.template attachSignalNode( + &omega_node); + + success *= (hygov.allocate() == 0); + hygov.y()[index(Var::PMECH)] = pmech_value; + success *= (hygov.verify() == 0); + success *= (hygov.initialize() == 0); + success *= (hygov.tagDifferentiable() == 0); + success *= (hygov.evaluateResidual() == 0); + + success *= isEqual(hygov.y()[index(Var::Q)], static_cast(0.50), kTol); + success *= isEqual(hygov.y()[index(Var::G)], static_cast(0.50), kTol); + success *= isEqual(hygov.y()[index(Var::C)], static_cast(0.50), kTol); + success *= (hygov.tag()[index(Var::G)] == false); + + for (size_t i = 0; i < hygov.getResidual().size(); ++i) + { + success *= isEqual(hygov.getResidual()[i], static_cast(0.0), kTol); + success *= isEqual(hygov.yp()[i], static_cast(0.0), kTol); + } + + return success.report(__func__); + } + + TestOutcome sourceDefault() + { + using Var = PhasorDynamics::Governor::HygovInternalVariables; + + TestStatus success = true; + + PhasorDynamics::SignalNode pmech_node; + ScalarT pmech_value{0.40}; + IdxT pmech_index = 4; + pmech_node.set(&pmech_value, &pmech_index); + + PhasorDynamics::Governor::Hygov hygov(makeHygovSourceDefaultData()); + hygov.getSignals().template assignSignalNode( + &pmech_node); + + success *= (hygov.allocate() == 0); + hygov.y()[index(Var::PMECH)] = pmech_value; + success *= (hygov.verify() == 0); + success *= (hygov.initialize() == 0); + success *= (hygov.tagDifferentiable() == 0); + success *= (hygov.evaluateResidual() == 0); + + success *= isEqual(hygov.y()[index(Var::Q)], static_cast(0.50), kTol); + success *= isEqual(hygov.y()[index(Var::PGV)], static_cast(0.50), kTol); + success *= isEqual(hygov.y()[index(Var::G)], static_cast(0.50), kTol); + success *= (hygov.tag()[index(Var::XN)] == false); + + for (size_t i = 0; i < hygov.getResidual().size(); ++i) + { + success *= isEqual(hygov.getResidual()[i], static_cast(0.0), kTol); + success *= isEqual(hygov.yp()[i], static_cast(0.0), kTol); + } + + return success.report(__func__); + } + + TestOutcome baseConversion() + { + using Var = PhasorDynamics::Governor::HygovInternalVariables; + using Params = PhasorDynamics::Governor::HygovParameters; + + TestStatus success = true; + + PhasorDynamics::SignalNode pmech_node; + ScalarT pmech_value{0.40}; + IdxT pmech_index = 6; + pmech_node.set(&pmech_value, &pmech_index); + + auto data = makeHygovData(); + data.parameters[Params::Trate] = static_cast(50.0); + + PhasorDynamics::Governor::Hygov hygov(data); + hygov.getSignals().template assignSignalNode(&pmech_node); + + success *= (hygov.allocate() == 0); + hygov.y()[index(Var::PMECH)] = pmech_value; + success *= (hygov.verify() == 0); + success *= (hygov.initialize() == 0); + success *= (hygov.evaluateResidual() == 0); + + success *= isEqual(hygov.y()[index(Var::Q)], static_cast(0.90), kTol); + success *= isEqual(hygov.y()[index(Var::G)], static_cast(0.90), kTol); + success *= isEqual(hygov.y()[index(Var::PMECH)], static_cast(0.40), kTol); + success *= isEqual(pmech_node.read(), static_cast(0.40), kTol); + + for (size_t i = 0; i < hygov.getResidual().size(); ++i) + { + success *= isEqual(hygov.getResidual()[i], static_cast(0.0), kTol); + success *= isEqual(hygov.yp()[i], static_cast(0.0), kTol); + } + + return success.report(__func__); + } + + TestOutcome jsonParseAndSystemAssembly() + { + TestStatus success = true; + + std::istringstream input(R"json( +{ + "header": { + "format_version": 0, + "format_revision": 1, + "case_name": "hydro governor", + "case_description": "HYGOV 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": 10, "name": "Pmech" } + ], + "devices": [ + { + "class": "Genrou", + "ports": { "bus": 1, "pmech": 10 }, + "id": "GEN1", + "params": { + "p0": 0.3, "q0": 0.0, "H": 3.0, "D": 0.0, "Ra": 0.0, + "Tdop": 7.0, "Tdopp": 0.04, "Tqop": 0.75, "Tqopp": 0.05, + "Xd": 2.1, "Xdp": 0.2, "Xdpp": 0.18, "Xq": 0.5, "Xqp": 0.5, + "Xqpp": 0.18, "Xl": 0.15, "S10": 0.0, "S12": 0.0, + "mva": 100.0 + } + }, + { + "class": "Hygov", + "ports": { "pmech": 10 }, + "id": "HYG1", + "params": { + "Trate": 50.0, "mva": 100.0, "Rperm": 0.05, "Rtemp": 0.4, "Tr": 5.0, + "Tf": 0.2, "Tg": 0.0, "Velm": 0.5, "Gmax": 1.0, "Gmin": 0.0, + "Tw": 1.0, "At": 1.0, "Dturb": 0.0, "Qnl": 0.1, + "Tn": 0.0, "Tnp": 1.0, "db1": 0.0, "db2": 0.0, "Hdam": 1.0, + "Gv0": 0.0, "Gv1": 0.2, "Gv2": 0.4, "Gv3": 0.6, "Gv4": 0.8, "Gv5": 1.0, + "Pgv0": 0.0, "Pgv1": 0.2, "Pgv2": 0.4, "Pgv3": 0.6, "Pgv4": 0.8, "Pgv5": 1.0 + } + } + ] +} +)json"); + + auto data = PhasorDynamics::parseSystemModelData(input); + success *= (data.hygov.size() == 1); + success *= (std::get( + data.hygov[0].parameters.at(PhasorDynamics::Governor::HygovParameters::Trate)) + == static_cast(50.0)); + success *= (std::get( + data.hygov[0].parameters.at(PhasorDynamics::Governor::HygovParameters::mva)) + == static_cast(100.0)); + PhasorDynamics::SystemModel system(data); + success *= (system.allocate() == 0); + success *= (system.initialize() == 0); + const auto hygov_offset = + static_cast(system.size() + - static_cast(PhasorDynamics::Governor::HygovInternalVariables::MAXIMUM)); + success *= isEqual(system.y()[hygov_offset + index(PhasorDynamics::Governor::HygovInternalVariables::Q)], + static_cast(0.70), + kTol); + success *= isEqual(system.y()[hygov_offset + index(PhasorDynamics::Governor::HygovInternalVariables::G)], + static_cast(0.70), + kTol); + success *= (system.evaluateResidual() == 0); + success *= (system.size() == 35); + + return success.report(__func__); + } + + private: + static size_t index(PhasorDynamics::Governor::HygovInternalVariables variable) + { + return static_cast(variable); + } + + auto makeHygovData() -> PhasorDynamics::Governor::HygovData + { + using Params = PhasorDynamics::Governor::HygovParameters; + using Mon = PhasorDynamics::Governor::HygovMonitorableVariables; + + PhasorDynamics::Governor::HygovData data; + data.device_class = "Hygov"; + data.disambiguation_string = "hygov_test"; + data.monitored_variables.insert(Mon::pmech); + data.monitored_variables.insert(Mon::gate); + + data.parameters[Params::Trate] = static_cast(100.0); + data.parameters[Params::mva] = static_cast(100.0); + data.parameters[Params::Rperm] = static_cast(0.05); + data.parameters[Params::Rtemp] = static_cast(0.4); + data.parameters[Params::Tr] = static_cast(5.0); + data.parameters[Params::Tf] = static_cast(0.2); + data.parameters[Params::Tg] = static_cast(0.0); + data.parameters[Params::Velm] = static_cast(0.5); + data.parameters[Params::Gmax] = static_cast(1.0); + data.parameters[Params::Gmin] = static_cast(0.0); + data.parameters[Params::Tw] = static_cast(1.0); + data.parameters[Params::At] = static_cast(1.0); + data.parameters[Params::Dturb] = static_cast(0.0); + data.parameters[Params::Qnl] = static_cast(0.1); + data.parameters[Params::Tn] = static_cast(0.0); + data.parameters[Params::Tnp] = static_cast(1.0); + data.parameters[Params::db1] = static_cast(0.0); + data.parameters[Params::db2] = static_cast(0.0); + data.parameters[Params::Hdam] = static_cast(1.0); + data.parameters[Params::Gv0] = static_cast(0.0); + data.parameters[Params::Gv1] = static_cast(0.2); + data.parameters[Params::Gv2] = static_cast(0.4); + data.parameters[Params::Gv3] = static_cast(0.6); + data.parameters[Params::Gv4] = static_cast(0.8); + data.parameters[Params::Gv5] = static_cast(1.0); + data.parameters[Params::Pgv0] = static_cast(0.0); + data.parameters[Params::Pgv1] = static_cast(0.2); + data.parameters[Params::Pgv2] = static_cast(0.4); + data.parameters[Params::Pgv3] = static_cast(0.6); + data.parameters[Params::Pgv4] = static_cast(0.8); + data.parameters[Params::Pgv5] = static_cast(1.0); + + return data; + } + + auto makeHygovSourceDefaultData() -> PhasorDynamics::Governor::HygovData + { + using Params = PhasorDynamics::Governor::HygovParameters; + + auto data = makeHygovData(); + data.parameters[Params::Tn] = static_cast(0.0); + data.parameters[Params::Tnp] = static_cast(0.0); + data.parameters[Params::Gv0] = static_cast(0.0); + data.parameters[Params::Gv1] = static_cast(0.0); + data.parameters[Params::Gv2] = static_cast(0.0); + data.parameters[Params::Gv3] = static_cast(0.0); + data.parameters[Params::Gv4] = static_cast(0.0); + data.parameters[Params::Gv5] = static_cast(0.0); + data.parameters[Params::Pgv0] = static_cast(0.0); + data.parameters[Params::Pgv1] = static_cast(0.0); + data.parameters[Params::Pgv2] = static_cast(0.0); + data.parameters[Params::Pgv3] = static_cast(0.0); + data.parameters[Params::Pgv4] = static_cast(0.0); + data.parameters[Params::Pgv5] = static_cast(0.0); + + return data; + } + }; + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp b/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp index c32df9f58..4d67d39f5 100644 --- a/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp +++ b/tests/UnitTests/PhasorDynamics/SystemSingleComponentTests.hpp @@ -166,6 +166,57 @@ namespace GridKit return success.report(__func__); } + TestOutcome regca() + { + TestStatus success = true; + + PhasorDynamics::SystemModel* system = new PhasorDynamics::SystemModel(); + + PhasorDynamics::BusInfinite bus(1.0, 0.0); + system->addBus(&bus); + + PhasorDynamics::Converter::Regca regca(&bus, makeRegcaData()); + system->addComponent(®ca); + + success *= system->allocate() == 0; + success *= system->initialize() == 0; + success *= system->evaluateResidual() == 0; + success *= system->evaluateJacobian() == 0; + success *= system->size() == regca.size(); + + delete system; + system = nullptr; + + return success.report(__func__); + } + + private: + auto makeRegcaData() -> PhasorDynamics::Converter::RegcaData + { + using Params = PhasorDynamics::Converter::RegcaParameters; + + PhasorDynamics::Converter::RegcaData data; + data.device_class = "Regca"; + data.disambiguation_string = "regca_test"; + data.parameters[Params::P0] = static_cast(1.0); + data.parameters[Params::Q0] = static_cast(0.0); + data.parameters[Params::mva] = static_cast(100.0); + data.parameters[Params::Tg] = static_cast(0.02); + data.parameters[Params::TM] = static_cast(0.02); + data.parameters[Params::Rqmax] = static_cast(999.0); + data.parameters[Params::Rqmin] = static_cast(-999.0); + data.parameters[Params::Rpmax] = static_cast(999.0); + data.parameters[Params::sL] = true; + data.parameters[Params::IL1] = static_cast(1.1); + data.parameters[Params::VL0] = static_cast(0.4); + data.parameters[Params::VL1] = static_cast(0.9); + data.parameters[Params::VA0] = static_cast(0.4); + data.parameters[Params::VA1] = static_cast(0.9); + data.parameters[Params::Vhvmax] = static_cast(1.2); + return data; + } + + public: TestOutcome genrou() { TestStatus success = true; diff --git a/tests/UnitTests/PhasorDynamics/runConverterRegcaTests.cpp b/tests/UnitTests/PhasorDynamics/runConverterRegcaTests.cpp new file mode 100644 index 000000000..2e9c1c482 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/runConverterRegcaTests.cpp @@ -0,0 +1,27 @@ +#include "ConverterRegcaTests.hpp" + +int main() +{ + GridKit::Testing::TestingResults result; + + GridKit::Testing::ConverterRegcaTests test; + + result += test.constructionAndValidation(); + result += test.parameterValidation(); + result += test.initializesFromPowerFlowAndPublishesSignals(); + result += test.unconnectedCommandsRemainConstant(); + result += test.externalCommandsDriveRuntimeResidual(); + result += test.invalidInitialization(); + result += test.signalVerification(); + result += test.nullBusVerification(); + result += test.busInjectionUsesSystemBase(); + result += test.residualEquations(); + result += test.highVoltageReactiveCurrentRoot(); + result += test.limiterBranchCoverage(); + result += test.jsonParseAndSystemAssembly(); +#ifdef GRIDKIT_ENABLE_ENZYME + result += test.jacobian(); +#endif + + return result.summary(); +} diff --git a/tests/UnitTests/PhasorDynamics/runGovernorHygovTests.cpp b/tests/UnitTests/PhasorDynamics/runGovernorHygovTests.cpp new file mode 100644 index 000000000..45154c841 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/runGovernorHygovTests.cpp @@ -0,0 +1,16 @@ +#include "GovernorHygovTests.hpp" + +int main() +{ + GridKit::Testing::TestingResults result; + + GridKit::Testing::GovernorHygovTests test; + + result += test.constructionAndValidation(); + result += test.signals(); + result += test.sourceDefault(); + result += test.baseConversion(); + result += test.jsonParseAndSystemAssembly(); + + return result.summary(); +} diff --git a/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp b/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp index fc4f07e80..2eb8cecb8 100644 --- a/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp +++ b/tests/UnitTests/PhasorDynamics/runSystemSingleComponentTests.cpp @@ -14,6 +14,7 @@ int main() result += test.ieeet1(); result += test.load(); result += test.loadZIP(); + result += test.regca(); result += test.genrou(); result += test.genClassical(); result += test.tgov1();