diff --git a/CHANGELOG.md b/CHANGELOG.md index 066786e09..6911fb4fd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -62,6 +62,7 @@ - Added cmake-format hooks, including in pre-commit. - Added off-nominal tap ratio and phase shift support to the PhasorDynamics `Branch` model. - Added portable Vector class to GridKit +- Added `REGCA` converter model implementation for PhasorDynamics. ## v0.1 diff --git a/GridKit/CommonMath.hpp b/GridKit/CommonMath.hpp index e4a7ec87f..8f9f8c5cf 100644 --- a/GridKit/CommonMath.hpp +++ b/GridKit/CommonMath.hpp @@ -162,7 +162,31 @@ namespace GridKit } /** - * @brief Smooth two-sided deadband function + * @brief Smooth Type 1 no-offset two-sided deadband function + * + * Smooth approximation to a deadband that returns zero inside the band and + * passes the input through unchanged outside the band. + * + * @tparam ScalarT - scalar data type + * @tparam RealT - Real data type (see GridKit::ScalarTraits::RealT) + * + * @param[in] x - Input signal + * @param[in] lower - Lower breakpoint + * @param[in] upper - Upper breakpoint + * @return Smooth no-offset deadbanded value + */ + template + __attribute__((always_inline)) inline ScalarT deadband1( + const ScalarT x, + const RealT lower, + const RealT upper) + { + assert(lower <= upper); + return x * (sigmoid(lower - x) + sigmoid(x - upper)); + } + + /** + * @brief Smooth Type 2 offset two-sided deadband function * * Smooth approximation to x - min(max(x, lower), upper), composed from the * smooth ramp function. @@ -173,10 +197,10 @@ namespace GridKit * @param[in] x - Input signal * @param[in] lower - Lower breakpoint * @param[in] upper - Upper breakpoint - * @return Smooth deadbanded value + * @return Smooth offset deadbanded value */ template - __attribute__((always_inline)) inline ScalarT deadband( + __attribute__((always_inline)) inline ScalarT deadband2( const ScalarT x, const RealT lower, const RealT upper) diff --git a/GridKit/Model/PhasorDynamics/CMakeLists.txt b/GridKit/Model/PhasorDynamics/CMakeLists.txt index 90ddb2abc..73d80e510 100644 --- a/GridKit/Model/PhasorDynamics/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/CMakeLists.txt @@ -33,6 +33,7 @@ add_subdirectory(Branch) add_subdirectory(Bus) add_subdirectory(BusFault) add_subdirectory(BusToSignalAdapter) +add_subdirectory(Converter) add_subdirectory(Exciter) add_subdirectory(Governor) add_subdirectory(Load) diff --git a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp index e2ae075d2..8d989fede 100644 --- a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp +++ b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include diff --git a/GridKit/Model/PhasorDynamics/Converter/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Converter/CMakeLists.txt new file mode 100644 index 000000000..cafb2cc36 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/CMakeLists.txt @@ -0,0 +1,6 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +add_subdirectory(REGCA) diff --git a/GridKit/Model/PhasorDynamics/Converter/REECA/README.md b/GridKit/Model/PhasorDynamics/Converter/REECA/README.md index 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..0f6a2aca6 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.cpp @@ -0,0 +1,27 @@ +/** + * @file Regca.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Non-Enzyme instantiation for the REGCA converter model. + */ + +#include "RegcaImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + template + int Regca::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Regca..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Regca; + template class Regca; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp new file mode 100644 index 000000000..0a91d5b0c --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/Regca.hpp @@ -0,0 +1,197 @@ +/** + * @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 ScalarT = scalar_type; + using IdxT = index_type; + using RealT = typename Component::RealT; + using BusT = BusBase; + using SignalT = SignalNode; + using ModelDataT = RegcaData; + using MonitorT = Model::VariableMonitor; + + Regca(BusT* bus); + Regca(BusT* bus, const ModelDataT& data); + ~Regca(); + + int setGridKitComponentID(IdxT) override final; + int allocate() override final; + int verify() const override final; + int initialize() override final; + int tagDifferentiable() override final; + int 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 ModelDataT& data); + void initializeMonitor(); + void setDerivedParameters(); + + ScalarT toComponentBase(ScalarT value) const + { + return value * va_system_base_ / va_converter_base_; + } + + ScalarT toSystemBase(ScalarT value) const + { + return value / toComponentBase(static_cast(ONE)); + } + + ScalarT 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(); + } + + BusT* bus_{nullptr}; + + RealT P0_{0}; + RealT Q0_{0}; + RealT mva_base_{0}; + RealT Tg_{0}; + RealT TM_{0}; + RealT Rqmax_{0}; + RealT Rqmin_{0}; + RealT Rpmax_{0}; + bool sL_{false}; + RealT IL1_{0}; + RealT VL0_{0}; + RealT VL1_{0}; + RealT VA0_{0}; + RealT VA1_{0}; + RealT Vhvmax_{0}; + IdxT bus_id_{0}; + + IdxT parameter_error_count_{0}; + RealT Mp_{0}; + RealT va_converter_base_{0}; + RealT use_lvpl_{0}; + RealT bypass_lvpl_{1}; + RealT iq_use_upper_{0}; + RealT iq_use_lower_{1}; + + ScalarT ipcmd_set_{0}; + ScalarT iqcmd_set_{0}; + + ComponentSignals signals_; + std::unique_ptr monitor_; + + std::vector ws_; + std::vector ws_indices_; + }; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaData.hpp new file mode 100644 index 000000000..998049475 --- /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..bf470cc26 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaDependencyTracking.cpp @@ -0,0 +1,27 @@ +/** + * @file RegcaDependencyTracking.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Dependency-tracking instantiations for the REGCA converter model. + */ + +#include "RegcaImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + template + int Regca::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Regca..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Regca; + template class Regca; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaEnzyme.cpp new file mode 100644 index 000000000..41b4a51b6 --- /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..0a195e0c9 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REGCA/RegcaImpl.hpp @@ -0,0 +1,575 @@ +/** + * @file RegcaImpl.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Definition of the REGCA phasor-dynamics converter model. + */ + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + using Log = ::GridKit::Utilities::Logger; + + template + Regca::Regca(BusT* bus) + : bus_(bus) + { + size_ = static_cast(RegcaInternalVariables::MAXIMUM); + } + + template + Regca::Regca(BusT* bus, const ModelDataT& data) + : bus_(bus), + monitor_(std::make_unique(data)) + { + initializeParameters(data); + initializeMonitor(); + size_ = static_cast(RegcaInternalVariables::MAXIMUM); + } + + template + Regca::~Regca() + { + } + + template + void Regca::setDerivedParameters() + { + Mp_ = static_cast(100.0) * Rpmax_; + use_lvpl_ = ZERO; + bypass_lvpl_ = ONE; + if (sL_) + { + use_lvpl_ = ONE; + bypass_lvpl_ = ZERO; + } + iq_use_upper_ = ZERO; + iq_use_lower_ = ONE; + va_converter_base_ = mva_base_ * static_cast(1.0e6); + } + + template + scalar_type Regca::activeCurrentLowerRateBound( + scalar_type ip) const + { + return -Rpmax_ - (Mp_ - Rpmax_) * Math::sigmoid(ip); + } + + template + scalar_type Regca::activeCurrentUpperRateBound( + scalar_type ip, + scalar_type il) const + { + const ScalarT sigma_ip = Math::sigmoid(ip); + return Mp_ * (ONE - sigma_ip) + + Rpmax_ * sigma_ip * (bypass_lvpl_ + use_lvpl_ * Math::sigmoid(il - ip)); + } + + template + void Regca::initializeParameters(const ModelDataT& data) + { + using Params = typename ModelDataT::Parameters; + using Ports = typename ModelDataT::Ports; + + parameter_error_count_ = 0; + + auto load_required_real = [&](auto key, RealT& target, const char* name) + { + if (!data.parameters.contains(key)) + { + Log::error() << "Regca: missing required parameter '" << name << "'\n"; + ++parameter_error_count_; + return; + } + + const auto& value = data.parameters.at(key); + if (const auto* real_value = std::get_if(&value)) + { + target = *real_value; + } + else if (const auto* index_value = std::get_if(&value)) + { + target = static_cast(*index_value); + } + else + { + Log::error() << "Regca: parameter '" << name << "' must be numeric\n"; + ++parameter_error_count_; + } + }; + + auto load_required_switch = [&](auto key, bool& target, const char* name) + { + if (!data.parameters.contains(key)) + { + Log::error() << "Regca: missing required parameter '" << name << "'\n"; + ++parameter_error_count_; + return; + } + + const auto& value = data.parameters.at(key); + if (const auto* bool_value = std::get_if(&value)) + { + target = *bool_value; + } + else if (const auto* index_value = std::get_if(&value); + index_value && (*index_value == 0 || *index_value == 1)) + { + target = (*index_value == 1); + } + else + { + Log::error() << "Regca: parameter '" << name << "' must be bool or 0/1\n"; + ++parameter_error_count_; + } + }; + + load_required_real(Params::P0, P0_, "P0"); + load_required_real(Params::Q0, Q0_, "Q0"); + load_required_real(Params::mva, mva_base_, "mva"); + load_required_real(Params::Tg, Tg_, "Tg"); + load_required_real(Params::TM, TM_, "TM"); + load_required_real(Params::Rqmax, Rqmax_, "Rqmax"); + load_required_real(Params::Rqmin, Rqmin_, "Rqmin"); + load_required_real(Params::Rpmax, Rpmax_, "Rpmax"); + load_required_switch(Params::sL, sL_, "sL"); + load_required_real(Params::IL1, IL1_, "IL1"); + load_required_real(Params::VL0, VL0_, "VL0"); + load_required_real(Params::VL1, VL1_, "VL1"); + load_required_real(Params::VA0, VA0_, "VA0"); + load_required_real(Params::VA1, VA1_, "VA1"); + load_required_real(Params::Vhvmax, Vhvmax_, "Vhvmax"); + + if (data.ports.contains(Ports::bus)) + { + bus_id_ = data.ports.at(Ports::bus); + } + + setDerivedParameters(); + } + + template + const Model::VariableMonitorBase* Regca::getMonitor() const + { + return monitor_.get(); + } + + template + void Regca::initializeMonitor() + { + using Variable = typename ModelDataT::MonitorableVariables; + auto index = [](RegcaInternalVariables variable) + { + return static_cast(variable); + }; + + monitor_->set(Variable::ir, [this, index] + { return y_[index(RegcaInternalVariables::IR)]; }); + monitor_->set(Variable::ii, [this, index] + { return y_[index(RegcaInternalVariables::II)]; }); + monitor_->set(Variable::p, [this, index] + { return y_[index(RegcaInternalVariables::PBR)]; }); + monitor_->set(Variable::q, [this, index] + { return y_[index(RegcaInternalVariables::QBR)]; }); + 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/INPUT_FORMAT.md b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md index 32df36cdc..b33849140 100644 --- a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md +++ b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md @@ -144,17 +144,17 @@ are specified: `Genrou` | 6th order machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Tdop`, `Tdopp`, `Tqop`, `Tqopp`, `Xd`, `Xdp`, `Xdpp`, `Xq`, `Xqp`, `Xqpp`, `Xl`, `S10`, `S12`, `mva` | `ir`, `ii`, `p`, `q`, `delta`, `omega`, `speed` `Gensal` | 5th order salient-pole machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Tdop`, `Tdopp`, `Tqopp`, `Xd`, `Xdp`, `Xdpp`, `Xq`, `Xl`, `S10`, `S12`, `mva` | `ir`, `ii`, `p`, `q`, `delta`, `omega`, `speed`, `Eqp`, `psidp`, `psiqpp`, `psidpp`, `vd`, `vq`, `te`, `id`, `iq` `GenClassical` | the classical machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Xdp`, `mva` | `ir`, `ii`, `p`, `q`, `delta`, `omega` + `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` `Ieeet1` | the IEEET1 exciter model | `bus`, `speed`, `efd`, `vs`\* | `Tr`, `Ka`, `Ta`, `Ke`, `Te`, `Kf`, `Tf`, `Vrmin`, `Vrmax`, `E1`, `E2`, `Se1`, `Se2`, `Ispdlim` | `efd`, `ksat` `SexsPti` | the SEXS-PTI simplified exciter model | `bus`, `efd`, `vs`\* | `Ta`, `Tb`, `Te`, `K`, `Efdmax`, `Efdmin` | `efd` - `Ieeest` | the IEEEST stabilizer model | `input`, `output` | `A1`, `A2`, `A3`, `A4`, `A5`, `A6`, `T1`, `T2`, `T3`, `T4`, `T5`, `T6`, `Ks`, `Lsmin`, `Lsmax`, `Vcl`, `Vcu`, `Tdelay` | `vss` + `Ieeest` | the IEEEST stabilizer model | `input`, `output` | `A1`, `A2`, `A3`, `A4`, `A5`, `A6`, `T1`, `T2`, `T3`, `T4`, `T5`, `T6`, `Ks`, `Lsmin`, `Lsmax`, `Vcl`, `Vcu`, `Tdelay` | `vss` `BusFault` | simple impedance-based fault at a bus | `bus`, `status`\* | `state0`, `R`, `X` | `state`, `ir`, `ii` `BusToSignalAdapter` | signal adapter component for a bus | `bus`, `vr`, `vi`, `ir`, `ii` | | Ports marked with \* are optional and, if missing, will be assumed to be connected to a constant value. This list is subject to change. - For `Branch`, `tap` and `phase` are optional parameters. If omitted, `tap` defaults to `1.0` and `phase` defaults to `0.0` radians. Bus `bus1` is the tap side for off-nominal transformer branches. diff --git a/GridKit/Model/PhasorDynamics/SystemModelData.hpp b/GridKit/Model/PhasorDynamics/SystemModelData.hpp index d4deef66e..7cb2765c5 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelData.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelData.hpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -40,6 +41,7 @@ namespace GridKit using BusDataT = BusData; using BusToSignalAdapterDataT = BusToSignalAdapterData; using BusFaultDataT = BusFaultData; + using RegcaDataT = Converter::RegcaData; using Tgov1DataT = Governor::Tgov1Data; using Ieeet1DataT = Exciter::Ieeet1Data; using SexsPtiDataT = Exciter::SexsPtiData; @@ -93,6 +95,7 @@ namespace GridKit std::vector adapter; ///< bus-to-signal adapters within the model std::vector branch; ///< Branches within the model std::vector bus_fault; ///< Bus faults within the model + std::vector 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 diff --git a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp index 00a6278c2..b215c0c86 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp @@ -136,6 +136,12 @@ namespace GridKit raw_component.get_to(loadzip); sm.loadzip.push_back(loadzip); } + else if (kind == "Regca") + { + typename SystemModelData::RegcaDataT regca; + raw_component.get_to(regca); + sm.regca.push_back(regca); + } else if (kind == "Tgov1") { typename SystemModelData::Tgov1DataT gov; diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index 891c2c57d..c28e7864f 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -46,6 +46,7 @@ namespace GridKit using namespace Governor; using namespace Exciter; using namespace Stabilizer; + using namespace Converter; // Set system model tolerances rel_tol_ = 1e-7; @@ -115,6 +116,62 @@ namespace GridKit addComponent(adapter); } + // Add REGCA converters + for (const auto& regcadata : data.regca) + { + IdxT bus_index = 0; + if (regcadata.ports.contains(RegcaPorts::bus)) + { + bus_index = regcadata.ports.at(RegcaPorts::bus); + } + + auto* regca = new Regca(getBus(bus_index), regcadata); + + if (regcadata.ports.contains(RegcaPorts::ipcmd)) + { + IdxT ipcmd = regcadata.ports.at(RegcaPorts::ipcmd); + constexpr auto IPCMD = RegcaExternalVariables::IPCMD; + regca->getSignals().template attachSignalNode(getSignal(ipcmd)); + } + + if (regcadata.ports.contains(RegcaPorts::iqcmd)) + { + IdxT iqcmd = regcadata.ports.at(RegcaPorts::iqcmd); + constexpr auto IQCMD = RegcaExternalVariables::IQCMD; + regca->getSignals().template attachSignalNode(getSignal(iqcmd)); + } + + if (regcadata.ports.contains(RegcaPorts::ibranchr)) + { + IdxT ibranchr = regcadata.ports.at(RegcaPorts::ibranchr); + constexpr auto IR = RegcaInternalVariables::IR; + regca->getSignals().template assignSignalNode(getSignal(ibranchr)); + } + + if (regcadata.ports.contains(RegcaPorts::ibranchi)) + { + IdxT ibranchi = regcadata.ports.at(RegcaPorts::ibranchi); + constexpr auto II = RegcaInternalVariables::II; + regca->getSignals().template assignSignalNode(getSignal(ibranchi)); + } + + if (regcadata.ports.contains(RegcaPorts::pbranch)) + { + IdxT pbranch = regcadata.ports.at(RegcaPorts::pbranch); + constexpr auto PBR = RegcaInternalVariables::PBR; + regca->getSignals().template assignSignalNode(getSignal(pbranch)); + } + + if (regcadata.ports.contains(RegcaPorts::qbranch)) + { + IdxT qbranch = regcadata.ports.at(RegcaPorts::qbranch); + constexpr auto QBR = RegcaInternalVariables::QBR; + regca->getSignals().template assignSignalNode(getSignal(qbranch)); + } + + addComponent(regca); + } + // Add branches for (const auto& branchdata : data.branch) { diff --git a/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp b/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp index 1ce2d8b1f..9c85c6388 100644 --- a/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp +++ b/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp @@ -55,7 +55,31 @@ namespace GridKit return success.report(__func__); } - TestOutcome deadband() + TestOutcome deadband1() + { + TestStatus success = true; + + const ScalarT lower = scalar(-0.05); + const ScalarT upper = scalar(0.10); + const ScalarT tol = scalar(kSmoothTolerance); + + // Type 1 (no-offset) deadband returns ~0 inside the band ... + const ScalarT inside_input = scalar(0.02); + success *= (std::abs(Math::deadband1(inside_input, lower, upper)) < tol * tol); + + // ... and passes the input through unchanged (no offset) outside the band. + const ScalarT far_above_input = scalar(4.0); + const ScalarT far_below_input = scalar(-4.0); + + success *= std::isfinite(Math::deadband1(far_above_input, lower, upper)); + success *= within(Math::deadband1(far_above_input, lower, upper), far_above_input, tol); + success *= std::isfinite(Math::deadband1(far_below_input, lower, upper)); + success *= within(Math::deadband1(far_below_input, lower, upper), far_below_input, tol); + + return success.report(__func__); + } + + TestOutcome deadband2() { TestStatus success = true; @@ -69,12 +93,12 @@ namespace GridKit const ScalarT expected_below = below_input - lower; const ScalarT expected_above = above_input - upper; - success *= within(Math::deadband(below_input, lower, upper), expected_below, tol); - success *= (std::abs(Math::deadband(inside_input, lower, upper)) < tol * tol); - success *= within(Math::deadband(above_input, lower, upper), expected_above, tol); + success *= within(Math::deadband2(below_input, lower, upper), expected_below, tol); + success *= (std::abs(Math::deadband2(inside_input, lower, upper)) < tol * tol); + success *= within(Math::deadband2(above_input, lower, upper), expected_above, tol); - const ScalarT lower_breakpoint = Math::deadband(lower, lower, upper); - const ScalarT upper_breakpoint = Math::deadband(upper, lower, upper); + const ScalarT lower_breakpoint = Math::deadband2(lower, lower, upper); + const ScalarT upper_breakpoint = Math::deadband2(upper, lower, upper); success *= (lower_breakpoint < scalar(0.0)); success *= (upper_breakpoint > scalar(0.0)); @@ -82,7 +106,7 @@ namespace GridKit success *= (std::abs(upper_breakpoint) < tol); const ScalarT x = scalar(-0.4); - success *= (std::abs(Math::deadband(x, lower, upper) + success *= (std::abs(Math::deadband2(x, lower, upper) - (x - Math::clamp(x, lower, upper))) < scalar(kRoundoffTolerance)); @@ -91,17 +115,17 @@ namespace GridKit const ScalarT expected_far_above = far_above_input - upper; const ScalarT expected_far_below = far_below_input - lower; - success *= std::isfinite(Math::deadband(far_above_input, lower, upper)); - success *= within(Math::deadband(far_above_input, lower, upper), expected_far_above, tol); - success *= std::isfinite(Math::deadband(far_below_input, lower, upper)); - success *= within(Math::deadband(far_below_input, lower, upper), expected_far_below, tol); + success *= std::isfinite(Math::deadband2(far_above_input, lower, upper)); + success *= within(Math::deadband2(far_above_input, lower, upper), expected_far_above, tol); + success *= std::isfinite(Math::deadband2(far_below_input, lower, upper)); + success *= within(Math::deadband2(far_below_input, lower, upper), expected_far_below, tol); const ScalarT point = scalar(0.25); const ScalarT above_point = scalar(0.75); const ScalarT below_point = scalar(-0.25); - success *= (std::abs(Math::deadband(above_point, point, point) - (above_point - point)) + success *= (std::abs(Math::deadband2(above_point, point, point) - (above_point - point)) < scalar(kRoundoffTolerance)); - success *= (std::abs(Math::deadband(below_point, point, point) - (below_point - point)) + success *= (std::abs(Math::deadband2(below_point, point, point) - (below_point - point)) < scalar(kRoundoffTolerance)); return success.report(__func__); diff --git a/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp b/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp index 233c36ddd..aecdc7ae8 100644 --- a/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp +++ b/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp @@ -7,7 +7,8 @@ int main() GridKit::Testing::SmoothnessIndicatorTests test; result += test.clamp(); - result += test.deadband(); + result += test.deadband1(); + result += test.deadband2(); result += test.limitIndicators(); result += test.slew(); result += test.linseg(); diff --git a/tests/UnitTests/PhasorDynamics/CMakeLists.txt b/tests/UnitTests/PhasorDynamics/CMakeLists.txt index bac53ec7b..b09a2f450 100644 --- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt +++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt @@ -94,6 +94,14 @@ target_link_libraries( GridKit::phasor_dynamics_systemmodel_dependency_tracking GridKit::testing) +add_executable(test_phasor_converter_regca runConverterRegcaTests.cpp) +target_link_libraries( + test_phasor_converter_regca + GridKit::definitions + GridKit::phasor_dynamics_systemmodel + GridKit::phasor_dynamics_systemmodel_dependency_tracking + GridKit::testing) + add_executable(test_phasor_stabilizer_ieeest runStabilizerIeeestTests.cpp) target_link_libraries( test_phasor_stabilizer_ieeest @@ -138,6 +146,7 @@ add_test(NAME PhasorDynamicsGovernorTgov1Test COMMAND test_phasor_governor_tgov1 add_test(NAME PhasorDynamicsExciterIeeet1Test COMMAND test_phasor_exciter_ieeet1) add_test(NAME PhasorDynamicsGensalTest COMMAND test_phasor_gensal) add_test(NAME PhasorDynamicsExciterSexsPtiTest COMMAND test_phasor_exciter_sexspti) +add_test(NAME 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) @@ -159,6 +168,7 @@ install( test_phasor_exciter_ieeet1 test_phasor_gensal test_phasor_exciter_sexspti + test_phasor_converter_regca test_phasor_stabilizer_ieeest test_phasor_gen_classical test_phasor_system diff --git a/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp new file mode 100644 index 000000000..7ef9e971e --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/ConverterRegcaTests.hpp @@ -0,0 +1,951 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + template + class ConverterRegcaTests + { + public: + using ScalarT = scalar_type; + using IdxT = index_type; + using RealT = typename PhasorDynamics::Component::RealT; + + ConverterRegcaTests() = default; + ~ConverterRegcaTests() = default; + + static constexpr ScalarT kTol = static_cast(1.0e-9); + + TestOutcome constructionAndValidation() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + + PhasorDynamics::Converter::Regca minimal(&bus); + success *= (minimal.size() == static_cast(Vars::MAXIMUM)); + success *= (minimal.getMonitor() == nullptr); + success *= (minimal.verify() > 0); + + PhasorDynamics::Converter::Regca configured(&bus, makeData()); + success *= (configured.size() == static_cast(Vars::MAXIMUM)); + success *= (configured.getMonitor() != nullptr); + success *= (configured.verify() == 0); + + return success.report(__func__); + } + + TestOutcome parameterValidation() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + + auto missing = makeData(); + missing.parameters.erase(Params::Tg); + PhasorDynamics::Converter::Regca missing_model(&bus, missing); + success *= (missing_model.verify() > 0); + + auto missing_p0 = makeData(); + missing_p0.parameters.erase(Params::P0); + PhasorDynamics::Converter::Regca missing_p0_model(&bus, missing_p0); + success *= (missing_p0_model.verify() > 0); + + auto missing_q0 = makeData(); + missing_q0.parameters.erase(Params::Q0); + PhasorDynamics::Converter::Regca missing_q0_model(&bus, missing_q0); + success *= (missing_q0_model.verify() > 0); + + auto bad_switch = makeData(); + bad_switch.parameters[Params::sL] = static_cast(2); + PhasorDynamics::Converter::Regca bad_switch_model(&bus, bad_switch); + success *= (bad_switch_model.verify() > 0); + + success *= invalidParameterCase(bus, Params::mva, static_cast(0.0)); + success *= invalidParameterCase(bus, Params::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/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/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();