diff --git a/CHANGELOG.md b/CHANGELOG.md index ba0f1dde6..155919048 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -67,6 +67,7 @@ - Added `ConstantSignalSource` component as first use case for `BusToSignalAdapter`. - Added IDA option to suppress algebraic variables in local error tests. - Removed `COO_Matrix` class. +- Added `REPCA` 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 baf76ab90..91bb8d04b 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/Component.hpp b/GridKit/Model/PhasorDynamics/Component.hpp index 20dbd0403..8d2a33117 100644 --- a/GridKit/Model/PhasorDynamics/Component.hpp +++ b/GridKit/Model/PhasorDynamics/Component.hpp @@ -112,6 +112,15 @@ namespace GridKit return abs_tol_; } + int setAbsoluteTolerance(RealT rel_tol) override + { + for (IdxT i = 0; i < size_; ++i) + { + abs_tol_[static_cast(i)] = rel_tol; + } + return 0; + } + std::vector& getResidual() override { return f_; diff --git a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp index d8fe3b2dc..1e0a85124 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..be3287ed4 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/CMakeLists.txt @@ -0,0 +1,6 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +add_subdirectory(REPCA) diff --git a/GridKit/Model/PhasorDynamics/Converter/README.md b/GridKit/Model/PhasorDynamics/Converter/README.md index ad38ba19a..f14960a80 100644 --- a/GridKit/Model/PhasorDynamics/Converter/README.md +++ b/GridKit/Model/PhasorDynamics/Converter/README.md @@ -12,3 +12,4 @@ The GridKit converter documentation includes: - Renewable Energy Generator/Converter Model REGCA (See [REGCA](REGCA/README.md)) - Renewable Energy Generator/Converter Model REGCB (See [REGCB](REGCB/README.md)) - Renewable Energy Electrical Control Model REECA (See [REECA](REECA/README.md)) +- Renewable Energy Plant Control Model REPCA (See [REPCA](REPCA/README.md)) diff --git a/GridKit/Model/PhasorDynamics/Converter/REPCA/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Converter/REPCA/CMakeLists.txt new file mode 100644 index 000000000..c315b9c8d --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REPCA/CMakeLists.txt @@ -0,0 +1,54 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +set(_install_headers Repca.hpp RepcaData.hpp) + +if(GRIDKIT_ENABLE_ENZYME) + gridkit_add_library( + phasor_dynamics_converter_repca + SOURCES RepcaEnzyme.cpp + HEADERS ${_install_headers} + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal + PRIVATE + ClangEnzymeFlags + COMPILE_OPTIONS + PRIVATE + -mllvm + -enzyme-auto-sparsity=1 + -fno-math-errno) +else() + gridkit_add_library( + phasor_dynamics_converter_repca + SOURCES Repca.cpp + HEADERS ${_install_headers} + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal) +endif() + +gridkit_add_library( + phasor_dynamics_converter_repca_dependency_tracking + SOURCES RepcaDependencyTracking.cpp + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal_dependency_tracking) + +target_link_libraries( + phasor_dynamics_components + INTERFACE GridKit::phasor_dynamics_converter_repca) +target_link_libraries( + phasor_dynamics_components_dependency_tracking + INTERFACE GridKit::phasor_dynamics_converter_repca_dependency_tracking) diff --git a/GridKit/Model/PhasorDynamics/Converter/REPCA/README.md b/GridKit/Model/PhasorDynamics/Converter/REPCA/README.md new file mode 100644 index 000000000..7cc170d38 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REPCA/README.md @@ -0,0 +1,336 @@ +# **Renewable Energy Plant Control Model (REPCA)** + +REPCA is a WECC renewable energy plant control model for inverter-coupled +resources. In GridKit it is represented as a plant-level signal-control model +that computes active- and reactive-power commands for downstream electrical +control models. + +## Block Diagram + +Standard REPCA block diagram. + +![](../../../../../docs/Figures/PhasorDynamics_REPCA_Diagram.png) + +Figure 1: REPCA block diagram. Figure courtesy of [PowerWorld](https://www.powerworld.com/WebHelp/) + +## Model Parameters + +Symbol | Units | JSON | Description | Typical Value | Note +------------------------------------|----------|-------------|---------------------------------------------------------|---------------|------ +$S^\mathrm{base}$ | [MVA] | `mva` | REPCA component power base | 100.0 | +$s_{\mathrm{comp}}$ | [binary] | `VcompFlag` | Voltage-compensation mode flag | 1.0 | +$s_{\mathrm{ref}}$ | [binary] | `RefFlag` | Reactive-loop reference flag | 1.0 | +$s_{\mathrm{freq}}$ | [binary] | `Freqflag` | Active-power control output flag | 0.0 | +$T_{\mathrm{fltr}}$ | [sec] | `Tfltr` | Voltage and reactive-power measurement filter time constant | 0.05 | +$V^\mathrm{frz}$ | [p.u.] | `Vfrz` | Regulated-voltage threshold below which the reactive PI state freezes | 0.7 | +$R_c$ | [p.u.] | `Rc` | Line-drop compensation resistance | 0.0 | +$X_c$ | [p.u.] | `Xc` | Line-drop compensation reactance | 0.0 | +$K_c$ | [p.u.] | `Kc` | Reactive-current compensation gain | 1.0 | +$D_{\mathrm{bd1}}$ | [p.u.] | `dbdlow` | Lower reactive-loop deadband threshold | 0.0 | +$D_{\mathrm{bd2}}$ | [p.u.] | `dbdupper` | Upper reactive-loop deadband threshold | 0.0 | +$e^{\max}$ | [p.u.] | `emax` | Maximum reactive-loop error limit | 1.0 | +$e^{\min}$ | [p.u.] | `emin` | Minimum reactive-loop error limit | -1.0 | +$K_{\mathrm{p}}$ | [p.u.] | `Kp` | Reactive-power controller proportional gain | 10.0 | +$K_{\mathrm{i}}$ | [p.u./s] | `Ki` | Reactive-power controller integral gain | 10.0 | +$Q^{\max}$ | [p.u.] | `Qmax` | Maximum reactive-power command | 1.0 | +$Q^{\min}$ | [p.u.] | `Qmin` | Minimum reactive-power command | -1.0 | +$T_{\mathrm{ft}}$ | [sec] | `Tft` | Reactive-command lead time constant | 0.0 | +$T_{\mathrm{fv}}$ | [sec] | `Tfv` | Reactive-command lag time constant | 3.0 | +$T_{\mathrm{p}}$ | [sec] | `Tp` | Active-power measurement filter time constant | 0.0 | +$D_{\mathrm{bd1}}^f$ | [p.u.] | `fdbd1` | Lower frequency-error deadband threshold | 0.0 | +$D_{\mathrm{bd2}}^f$ | [p.u.] | `fdbd2` | Upper frequency-error deadband threshold | 0.0 | +$D_{\mathrm{dn}}$ | [p.u.] | `Ddn` | Down-frequency droop response gain | 20.0 | +$D_{\mathrm{up}}$ | [p.u.] | `Dup` | Up-frequency droop response gain | 0.0 | +$e_P^{\max}$ | [p.u.] | `femax` | Maximum active-power error limit | 1.0 | +$e_P^{\min}$ | [p.u.] | `femin` | Minimum active-power error limit | -1.0 | +$K_{\mathrm{pg}}$ | [p.u.] | `Kpg` | Active-power controller proportional gain | 10.0 | +$K_{\mathrm{ig}}$ | [p.u./s] | `Kig` | Active-power controller integral gain | 10.0 | +$P^{\max}$ | [p.u.] | `Pmax` | Maximum active-power command | 2.0 | +$P^{\min}$ | [p.u.] | `Pmin` | Minimum active-power command | 0.0 | +$T_{\mathrm{lag}}$ | [sec] | `Tlag` | Active-power command lag time constant | 3.0 | + +### Parameter Validation + +Invalid REPCA parameter sets are rejected by the following checks. Let $\epsilon_T=10^{-3}$. + +```math +\begin{aligned} + T &\leftarrow \max\!\left(T, \epsilon_T\right) + \quad T\in\{T_{\mathrm{fltr}},T_{\mathrm{p}},T_{\mathrm{lag}}\} \\ + S^\mathrm{base} + &> 0 \\ + s_{\mathrm{comp}}, s_{\mathrm{ref}}, s_{\mathrm{freq}} + &\in \{0,1\} \\ + T_{\mathrm{fv}} + &> 0 \\ + D_{\mathrm{bd1}} + &\le 0 \le D_{\mathrm{bd2}} \\ + e^{\min} + &\le 0 \le e^{\max} \\ + Q^{\min} + &\le Q^{\max} \\ + D_{\mathrm{bd1}}^f + &\le 0 \le D_{\mathrm{bd2}}^f \\ + e_P^{\min} + &\le 0 \le e_P^{\max} \\ + P^{\min} + &\le P^{\max} +\end{aligned} +``` + +### Model Derived Parameters + +```math +\begin{aligned} + s_{\mathrm{comp}}^\mathrm{off} &= 1 - s_{\mathrm{comp}} \\ + s_{\mathrm{ref}}^\mathrm{off} &= 1 - s_{\mathrm{ref}} \\ + k_{\mathrm{base}} &= \dfrac{S^\mathrm{sys}}{S^\mathrm{base}} +\end{aligned} +``` + +## Model Variables + +### Internal Variables + +#### Differential + +Symbol | Units | Description | Note +------------------------|--------|-------------------------------------|------ +$V^\mathrm{meas}$ | [p.u.] | Filtered regulated voltage | State 1 in Fig. 1; source label: `Vmeas` +$Q^\mathrm{meas}$ | [p.u.] | Filtered reactive-power signal | State 2 in Fig. 1; source label: `Qmeas` +$x_Q^\mathrm{PI}$ | [p.u.] | Reactive PI controller state | State 3 in Fig. 1; source label: `Reactive PI` +$x_Q^\mathrm{lag}$ | [p.u.] | Reactive-command lead-lag state | State 4 in Fig. 1; source label: `Qext` +$P^\mathrm{meas}$ | [p.u.] | Filtered active-power signal | State 5 in Fig. 1; source label: `Pmeas` +$x_P^\mathrm{PI}$ | [p.u.] | Active-power PI controller state | State 6 in Fig. 1; source label: `Power PI` +$P^\mathrm{ref}$ | [p.u.] | Active-power command lag state | State 7 in Fig. 1; source label: `Pref` + +#### Algebraic + +Symbol | Units | Description | Note +--------------------------------|--------|-------------------------------------|------ +$V$ | [p.u.] | Regulated-bus voltage magnitude | +$V^\mathrm{ldc}$ | [p.u.] | Line-drop compensated voltage magnitude | +$V^\mathrm{droop}$ | [p.u.] | Reactive-droop-compensated voltage | +$V^\mathrm{ctrl}$ | [p.u.] | Selected voltage-measurement input | +$s_{\mathrm{frz}}$ | [binary] | Reactive-PI voltage-enable indicator | +$e_{\mathrm{RQ}}$ | [p.u.] | Selected reactive-loop error | +$e_{\mathrm{RQ}}^\mathrm{db}$ | [p.u.] | Deadbanded reactive-loop error | +$e_{\mathrm{RQ}}^\mathrm{lim}$ | [p.u.] | Limited reactive-loop error | +$Q^\mathrm{PI}$ | [p.u.] | Reactive PI output | +$Q^\mathrm{ext}$ | [p.u.] | Reactive-power command output | System base +$e_f$ | [p.u.] | Frequency error after deadband | +$e_P$ | [p.u.] | Active-power control error | +$e_P^\mathrm{lim}$ | [p.u.] | Limited active-power control error | +$P^\mathrm{PI}$ | [p.u.] | Active-power PI output | +$P^\mathrm{ext}$ | [p.u.] | Active-power command output | System base + +### External Variables + +#### Differential +None. + +#### Algebraic + +Symbol | Units | Type | Description | Note +-------------------------------------|--------|--------|-----------------------------------|------ +$V_{\mathrm{r}}$ | [p.u.] | Known | Regulated-bus voltage, real component | Source label: `Vreg` +$V_{\mathrm{i}}$ | [p.u.] | Known | Regulated-bus voltage, imaginary component | Source label: `Vreg` +$I_{\mathrm{r}}$ | [p.u.] | Known | Branch current real component | Source label: `Ibranch` +$I_{\mathrm{i}}$ | [p.u.] | Known | Branch current imaginary component | Source label: `Ibranch` +$P^\mathrm{br}$ | [p.u.] | Known | Branch active power | +$Q^\mathrm{br}$ | [p.u.] | Known | Branch reactive power | +$f$ | [p.u.] | Known | Frequency input | Source label: `Freq` +$f^\mathrm{ref}$ | [p.u.] | Unknown | Frequency reference | Source label: `Freq_ref` +$V^\mathrm{ref}$ | [p.u.] | Unknown | Voltage-control reference | +$Q^\mathrm{ref}$ | [p.u.] | Unknown | Reactive-power reference | System base +$P_\mathrm{plant}^\mathrm{ref}$ | [p.u.] | Unknown | Plant active-power reference | System base + +## Model Equations + +### Differential Equations + +```math +\begin{aligned} + 0 &= + -\dot{V}^\mathrm{meas} + + \dfrac{1}{T_{\mathrm{fltr}}} + \left(V^\mathrm{ctrl} - V^\mathrm{meas}\right) \\ + 0 &= + -\dot{Q}^\mathrm{meas} + + \dfrac{1}{T_{\mathrm{fltr}}} + \left(k_{\mathrm{base}}Q^\mathrm{br} - Q^\mathrm{meas}\right) \\ + 0 &= + -\dot{x}_Q^\mathrm{PI} + + s_{\mathrm{frz}} + \text{antiwindup}\!\left( + Q^\mathrm{PI}, + K_{\mathrm{i}}e_{\mathrm{RQ}}^\mathrm{lim}, + Q^{\min}, + Q^{\max} + \right) \\ + 0 &= + -\dot{x}_Q^\mathrm{lag} + + \dfrac{1}{T_{\mathrm{fv}}} + \left(Q^\mathrm{PI} - x_Q^\mathrm{lag}\right) \\ + 0 &= + -\dot{P}^\mathrm{meas} + + \dfrac{1}{T_{\mathrm{p}}} + \left(k_{\mathrm{base}}P^\mathrm{br} - P^\mathrm{meas}\right) \\ + 0 &= + -\dot{x}_P^\mathrm{PI} + + \text{antiwindup}\!\left( + P^\mathrm{PI}, + K_{\mathrm{ig}}e_P^\mathrm{lim}, + P^{\min}, + P^{\max} + \right) \\ + 0 &= + -\dot{P}^\mathrm{ref} + + \dfrac{1}{T_{\mathrm{lag}}} + \left(P^\mathrm{PI} - P^\mathrm{ref}\right) +\end{aligned} +``` + +CommonMath defines the [Anti-Windup](../../../../CommonMath.md#anti-windup-indicator) target and smooth approximation. + +### Algebraic Equations + +```math +\begin{aligned} + 0 &= -V^2 + + V_{\mathrm{r}}^2 + + V_{\mathrm{i}}^2 \\ + 0 &= + -\left(V^\mathrm{ldc}\right)^2 + + \left(V_{\mathrm{r}} - R_c I_{\mathrm{r}} + X_c I_{\mathrm{i}}\right)^2 + + \left(V_{\mathrm{i}} - R_c I_{\mathrm{i}} - X_c I_{\mathrm{r}}\right)^2 \\ + 0 &= -V^\mathrm{droop} + V + K_c k_{\mathrm{base}}Q^\mathrm{br} \\ + 0 &= -V^\mathrm{ctrl} + s_{\mathrm{comp}}V^\mathrm{ldc} + s_{\mathrm{comp}}^\mathrm{off}V^\mathrm{droop} \\ + 0 &= -s_{\mathrm{frz}} + \text{above}(V, V^\mathrm{frz}) \\[0.5ex] + 0 &= -e_{\mathrm{RQ}} + + s_{\mathrm{ref}}\left(V^\mathrm{ref} - V^\mathrm{meas}\right) + + s_{\mathrm{ref}}^\mathrm{off}\left(k_{\mathrm{base}}Q^\mathrm{ref} - Q^\mathrm{meas}\right) \\ + 0 &= -e_{\mathrm{RQ}}^\mathrm{db} + + \text{deadband2}(e_{\mathrm{RQ}}, D_{\mathrm{bd1}}, D_{\mathrm{bd2}}) \\ + 0 &= -e_{\mathrm{RQ}}^\mathrm{lim} + + \text{clamp}(e_{\mathrm{RQ}}^\mathrm{db}, e^{\min}, e^{\max}) \\ + 0 &= -Q^\mathrm{PI} + \text{clamp}(K_{\mathrm{p}} e_{\mathrm{RQ}}^\mathrm{lim} + x_Q^\mathrm{PI}, Q^{\min}, Q^{\max}) \\ + 0 &= -T_{\mathrm{fv}}(k_{\mathrm{base}}Q^\mathrm{ext} - x_Q^\mathrm{lag}) + + T_{\mathrm{ft}}(Q^\mathrm{PI} - x_Q^\mathrm{lag}) \\[0.5ex] + 0 &= -e_f + \text{deadband2}(f^\mathrm{ref} - f, D_{\mathrm{bd1}}^f, D_{\mathrm{bd2}}^f) \\ + 0 &= -e_P + + k_{\mathrm{base}}P_\mathrm{plant}^\mathrm{ref} + - P^\mathrm{meas} + + D_{\mathrm{dn}}\rho(e_f) + - D_{\mathrm{up}}\rho(-e_f) \\ + 0 &= -e_P^\mathrm{lim} + \text{clamp}(e_P, e_P^{\min}, e_P^{\max}) \\ + 0 &= -P^\mathrm{PI} + \text{clamp}(K_{\mathrm{pg}} e_P^\mathrm{lim} + x_P^\mathrm{PI}, P^{\min}, P^{\max}) \\ + 0 &= -k_{\mathrm{base}}P^\mathrm{ext} + s_{\mathrm{freq}}P^\mathrm{ref} +\end{aligned} +``` + +CommonMath defines [above](../../../../CommonMath.md#derived-functions), +[clamp](../../../../CommonMath.md#derived-functions), +[deadband2](../../../../CommonMath.md#derived-functions), and +[ramp](../../../../CommonMath.md#primitives) $\rho$. + +## Initialization + +### External Priors + +```math +\begin{aligned} + V_{\mathrm{r},0}, V_{\mathrm{i},0} + &\leftarrow \text{regulated-bus voltage} \\ + I_{\mathrm{r},0}, I_{\mathrm{i},0} + &\leftarrow \text{branch current} \\ + P_0^\mathrm{br}, Q_0^\mathrm{br} + &\leftarrow \text{branch power} \\ + f_0 + &\leftarrow \text{frequency input} \\ + Q_0^\mathrm{ext} + &\leftarrow \text{reactive-power command start} \\ + P_0^\mathrm{ext} + &\leftarrow \text{active-power command start} +\end{aligned} +``` + +### Internal Initialization + +```math +\begin{aligned} + V_0 + &= \sqrt{ + V_{\mathrm{r},0}^2 + + V_{\mathrm{i},0}^2 + } \\ + V_0^\mathrm{ldc} + &= \sqrt{ + \left(V_{\mathrm{r},0} - R_c I_{\mathrm{r},0} + X_c I_{\mathrm{i},0}\right)^2 + + \left(V_{\mathrm{i},0} - R_c I_{\mathrm{i},0} - X_c I_{\mathrm{r},0}\right)^2 + } \\ + V_0^\mathrm{droop} + &= V_0 + + K_c k_{\mathrm{base}}Q_0^\mathrm{br} \\ + V_0^\mathrm{ctrl} + &= s_{\mathrm{comp}}V_0^\mathrm{ldc} + + s_{\mathrm{comp}}^\mathrm{off}V_0^\mathrm{droop} \\ + e_{f,0} &= 0 \\ + P_0^\mathrm{freq} + &= D_{\mathrm{dn}}\rho(e_{f,0}) + - D_{\mathrm{up}}\rho(-e_{f,0}) \\ + V_0^\mathrm{meas} &= V_0^\mathrm{ctrl} \\ + Q_0^\mathrm{meas} &= k_{\mathrm{base}}Q_0^\mathrm{br} \\ + P_0^\mathrm{meas} &= k_{\mathrm{base}}P_0^\mathrm{br} \\ + s_{\mathrm{frz},0} + &= \text{above}\left(V_0, V^\mathrm{frz}\right) \\ + e_{\mathrm{RQ},0} &= 0 \\ + e_{\mathrm{RQ},0}^\mathrm{db} + &= \text{deadband2}\left(e_{\mathrm{RQ},0}, D_{\mathrm{bd1}}, D_{\mathrm{bd2}}\right) \\ + e_{\mathrm{RQ},0}^\mathrm{lim} + &= \text{clamp}\left(e_{\mathrm{RQ},0}^\mathrm{db}, e^{\min}, e^{\max}\right) \\ + Q_0^\mathrm{PI} &= k_{\mathrm{base}}Q_0^\mathrm{ext} \\ + x_{Q,0}^\mathrm{lag} &= k_{\mathrm{base}}Q_0^\mathrm{ext} \\ + x_{Q,0}^\mathrm{PI} + &= k_{\mathrm{base}}Q_0^\mathrm{ext} - K_{\mathrm{p}}e_{\mathrm{RQ},0}^\mathrm{lim} \\ + e_{P,0} &= 0 \\ + e_{P,0}^\mathrm{lim} + &= \text{clamp}\left(e_{P,0}, e_P^{\min}, e_P^{\max}\right) \\ + P_0^\mathrm{ref} + &= + \begin{cases} + k_{\mathrm{base}}P_0^\mathrm{ext} & s_{\mathrm{freq}} = 1 \\ + \text{clamp} + \left(k_{\mathrm{base}}P_0^\mathrm{br}, P^{\min}, P^{\max}\right) + & s_{\mathrm{freq}} = 0 + \end{cases} \\ + P_0^\mathrm{PI} &= P_0^\mathrm{ref} \\ + x_{P,0}^\mathrm{PI} + &= P_0^\mathrm{ref} - K_{\mathrm{pg}}e_{P,0}^\mathrm{lim} +\end{aligned} +``` + +### External Solved + +```math +\begin{aligned} + f_0^\mathrm{ref} &\leftarrow f_0 \\ + V_0^\mathrm{ref} &\leftarrow V_0^\mathrm{meas} \\ + Q_0^\mathrm{ref} &\leftarrow \dfrac{Q_0^\mathrm{meas}}{k_{\mathrm{base}}} \\ + P_{\mathrm{plant},0}^\mathrm{ref} + &\leftarrow \dfrac{P_0^\mathrm{meas} - P_0^\mathrm{freq}}{k_{\mathrm{base}}} +\end{aligned} +``` + +Initialization rejects nonzero reactive-power or active-power PI antiwindup +rates. + +## Monitorable Outputs + +Output | Units | Description | Note +----------------|--------|-------------------------------------|------ +`qext` | [p.u.] | Reactive-power command output | $Q^\mathrm{ext}$ (system base) +`pext` | [p.u.] | Active-power command output | $P^\mathrm{ext}$ (system base) +`vmeas` | [p.u.] | Filtered regulated voltage | $V^\mathrm{meas}$ +`qmeas` | [p.u.] | Filtered reactive-power signal | $Q^\mathrm{meas}$ (component base) +`pmeas` | [p.u.] | Filtered active-power signal | $P^\mathrm{meas}$ (component base) diff --git a/GridKit/Model/PhasorDynamics/Converter/REPCA/Repca.cpp b/GridKit/Model/PhasorDynamics/Converter/REPCA/Repca.cpp new file mode 100644 index 000000000..56a711282 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REPCA/Repca.cpp @@ -0,0 +1,27 @@ +/** + * @file Repca.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Template instantiations for the REPCA plant-control model. + */ + +#include "RepcaImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + template + int Repca::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Repca..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Repca; + template class Repca; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REPCA/Repca.hpp b/GridKit/Model/PhasorDynamics/Converter/REPCA/Repca.hpp new file mode 100644 index 000000000..79affc7f9 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REPCA/Repca.hpp @@ -0,0 +1,204 @@ +/** + * @file Repca.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Declaration of the REPCA phasor-dynamics plant-control model. + */ + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + template + class BusBase; + + template + class SignalNode; + + namespace Converter + { + /// Internal variables of a `Repca`. + enum class RepcaInternalVariables : size_t + { + VMEAS, ///< Filtered regulated voltage + QMEAS, ///< Filtered reactive-power signal + XQPI, ///< Reactive PI controller state + XQLAG, ///< Reactive-command lead-lag state + PMEAS, ///< Filtered active-power signal + XPPI, ///< Active-power PI controller state + PREF, ///< Active-power command lag state + V, ///< Regulated-bus voltage magnitude + VLDC, ///< Line-drop compensated voltage magnitude + VDROOP, ///< Reactive-droop-compensated voltage + VCTRL, ///< Selected voltage-measurement input + SFRZ, ///< Reactive-PI voltage-enable indicator + ERQ, ///< Selected reactive-loop error + ERQDB, ///< Deadbanded reactive-loop error + ERQLIM, ///< Limited reactive-loop error + QPI, ///< Reactive PI output + QEXT, ///< Reactive-power command output + EF, ///< Frequency error after deadband + EP, ///< Active-power control error + EPLIM, ///< Limited active-power control error + PPI, ///< Active-power PI output + PEXT, ///< Active-power command output + MAXIMUM, + }; + + /// External variables of a `Repca`. + enum class RepcaExternalVariables : size_t + { + IBRANCHR, ///< Branch current real component + IBRANCHI, ///< Branch current imaginary component + PBRANCH, ///< Branch active power + QBRANCH, ///< Branch reactive power + FREQ, ///< Frequency input + FREQREF, ///< Frequency reference + VREF, ///< Voltage-control reference + QREF, ///< Reactive-power reference + PPLANTREF, ///< Plant active-power reference + MAXIMUM, + }; + + template + class Repca : public Component + { + using Component::abs_tol_; + using Component::alpha_; + using Component::f_; + using Component::gridkit_component_id_; + using Component::J_cols_buffer_; + using Component::J_rows_buffer_; + using Component::J_vals_buffer_; + using Component::nnz_; + using Component::residual_indices_; + using Component::size_; + using Component::tag_; + using Component::va_system_base_; + using Component::variable_indices_; + using Component::wb_; + using Component::y_; + using Component::yp_; + + public: + using ScalarT = scalar_type; + using IdxT = index_type; + using RealT = typename Component::RealT; + using BusT = BusBase; + using SignalT = SignalNode; + using ModelDataT = RepcaData; + using MonitorT = Model::VariableMonitor; + + Repca(BusT* bus); + Repca(BusT* bus, const ModelDataT& data); + ~Repca(); + + int setGridKitComponentID(IdxT) override final; + int allocate() override final; + int verify() const override final; + int initialize() override final; + int tagDifferentiable() override final; + int setAbsoluteTolerance(RealT) 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( + const ScalarT*, const ScalarT*, const ScalarT*, const ScalarT*, ScalarT*); + + private: + void initializeParameters(const ModelDataT& data); + void initializeMonitor(); + void setDerivedParameters(); + + ScalarT toComponentBase(ScalarT value) const; + ScalarT toSystemBase(ScalarT value) const; + + ScalarT& Vr() + { + return bus_->Vr(); + } + + ScalarT& Vi() + { + return bus_->Vi(); + } + + static constexpr RealT TIME_CONSTANT_MINIMUM = static_cast(1.0e-3); + + BusT* bus_{nullptr}; + + RealT mva_{static_cast(100.0)}; + bool VcompFlag_{true}; + bool RefFlag_{true}; + bool Freqflag_{false}; + RealT Tfltr_{static_cast(0.05)}; + RealT Vfrz_{static_cast(0.7)}; + RealT Rc_{ZERO}; + RealT Xc_{ZERO}; + RealT Kc_{ONE}; + RealT dbdlow_{ZERO}; + RealT dbdupper_{ZERO}; + RealT emax_{ONE}; + RealT emin_{-ONE}; + RealT Kp_{static_cast(10.0)}; + RealT Ki_{static_cast(10.0)}; + RealT Qmax_{ONE}; + RealT Qmin_{-ONE}; + RealT Tft_{ZERO}; + RealT Tfv_{static_cast(3.0)}; + RealT Tp_{ZERO}; + RealT fdbd1_{ZERO}; + RealT fdbd2_{ZERO}; + RealT Ddn_{static_cast(20.0)}; + RealT Dup_{ZERO}; + RealT femax_{ONE}; + RealT femin_{-ONE}; + RealT Kpg_{static_cast(10.0)}; + RealT Kig_{static_cast(10.0)}; + RealT Pmax_{static_cast(2.0)}; + RealT Pmin_{ZERO}; + RealT Tlag_{static_cast(3.0)}; + + IdxT parameter_error_count_{0}; + RealT va_component_base_{ZERO}; + RealT vcomp_on_{ONE}; + RealT vcomp_off_{ZERO}; + RealT ref_on_{ONE}; + RealT ref_off_{ZERO}; + RealT freq_on_{ZERO}; + + ScalarT freqref_set_{ONE}; + ScalarT vref_set_{ONE}; + ScalarT qref_set_{ZERO}; + ScalarT pplantref_set_{ZERO}; + + 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/REPCA/RepcaData.hpp b/GridKit/Model/PhasorDynamics/Converter/REPCA/RepcaData.hpp new file mode 100644 index 000000000..6deb36ece --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REPCA/RepcaData.hpp @@ -0,0 +1,95 @@ +/** + * @file RepcaData.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Modeling data for the REPCA converter plant-control model. + */ + +#pragma once + +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + /// Parameter keys for the REPCA plant-control model. + enum class RepcaParameters + { + mva, ///< REPCA component power base + VcompFlag, ///< Voltage-compensation mode flag + RefFlag, ///< Reactive-loop reference flag + Freqflag, ///< Active-power control output flag + Tfltr, ///< Voltage and reactive-power measurement filter time constant + Vfrz, ///< Reactive-PI freeze voltage threshold + Rc, ///< Line-drop compensation resistance + Xc, ///< Line-drop compensation reactance + Kc, ///< Reactive-current compensation gain + dbdlow, ///< Lower reactive-loop deadband threshold + dbdupper, ///< Upper reactive-loop deadband threshold + emax, ///< Maximum reactive-loop error limit + emin, ///< Minimum reactive-loop error limit + Kp, ///< Reactive-power controller proportional gain + Ki, ///< Reactive-power controller integral gain + Qmax, ///< Maximum reactive-power command + Qmin, ///< Minimum reactive-power command + Tft, ///< Reactive-command lead time constant + Tfv, ///< Reactive-command lag time constant + Tp, ///< Active-power measurement filter time constant + fdbd1, ///< Lower frequency-error deadband threshold + fdbd2, ///< Upper frequency-error deadband threshold + Ddn, ///< Down-frequency droop response gain + Dup, ///< Up-frequency droop response gain + femax, ///< Maximum active-power error limit + femin, ///< Minimum active-power error limit + Kpg, ///< Active-power controller proportional gain + Kig, ///< Active-power controller integral gain + Pmax, ///< Maximum active-power command + Pmin, ///< Minimum active-power command + Tlag ///< Active-power command lag time constant + }; + + /// Ports for the REPCA plant-control model. + enum class RepcaPorts + { + bus, ///< Regulated bus ID + ibranchr, ///< Branch current real-component signal ID + ibranchi, ///< Branch current imaginary-component signal ID + pbranch, ///< Branch active-power signal ID + qbranch, ///< Branch reactive-power signal ID + freq, ///< Frequency input signal ID + freqref, ///< Frequency-reference signal ID + vref, ///< Voltage-reference signal ID + qref, ///< Reactive-power reference signal ID + pplantref, ///< Plant active-power reference signal ID + qext, ///< Reactive-power command output signal ID + pext ///< Active-power command output signal ID + }; + + /// Variables available through the monitor interface. + enum class RepcaMonitorableVariables + { + qext, ///< Reactive-power command output + pext, ///< Active-power command output + vmeas, ///< Filtered regulated voltage + qmeas, ///< Filtered reactive-power signal + pmeas ///< Filtered active-power signal + }; + + template + struct RepcaData : public ComponentData + { + RepcaData() = default; + + using Parameters = RepcaParameters; + using Ports = RepcaPorts; + using MonitorableVariables = RepcaMonitorableVariables; + }; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REPCA/RepcaDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Converter/REPCA/RepcaDependencyTracking.cpp new file mode 100644 index 000000000..7b8365122 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REPCA/RepcaDependencyTracking.cpp @@ -0,0 +1,27 @@ +/** + * @file RepcaDependencyTracking.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Dependency-tracking instantiations for the REPCA plant-control model. + */ + +#include "RepcaImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + template + int Repca::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Repca..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Repca; + template class Repca; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REPCA/RepcaEnzyme.cpp b/GridKit/Model/PhasorDynamics/Converter/REPCA/RepcaEnzyme.cpp new file mode 100644 index 000000000..9e7cb3abf --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REPCA/RepcaEnzyme.cpp @@ -0,0 +1,109 @@ +/** + * @file RepcaEnzyme.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Enzyme sparse Jacobian for the REPCA plant-control model. + */ + +#include + +#include "RepcaImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + template + int Repca::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Repca..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + + if (J_rows_buffer_ == nullptr) + { + auto size = static_cast(size_); + auto bus_size = static_cast(bus_->size()); + auto signal_size = static_cast(ws_.size()); + auto buffer_size = 2 * size * size + size * bus_size + size * signal_size; + J_rows_buffer_ = new IdxT[buffer_size]; + J_cols_buffer_ = new IdxT[buffer_size]; + J_vals_buffer_ = new RealT[buffer_size]; + } + + using RepcaT = GridKit::PhasorDynamics::Converter::Repca; + using Fn = GridKit::Enzyme::Sparse::MemberFunctions; + + nnz_ = 0; + + GridKit::Enzyme::Sparse::DfDy::eval( + this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + GridKit::Enzyme::Sparse::DfDyp::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_, + nnz_); + + 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_, + nnz_); + + 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_, + nnz_); + + this->constructCoo(); + + return 0; + } + + template class Repca; + template class Repca; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REPCA/RepcaImpl.hpp b/GridKit/Model/PhasorDynamics/Converter/REPCA/RepcaImpl.hpp new file mode 100644 index 000000000..a623b5646 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REPCA/RepcaImpl.hpp @@ -0,0 +1,703 @@ +/** + * @file RepcaImpl.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Definition of the REPCA phasor-dynamics plant-control model. + */ + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + using Log = ::GridKit::Utilities::Logger; + + template + Repca::Repca(BusT* bus) + : bus_(bus) + { + size_ = static_cast(RepcaInternalVariables::MAXIMUM); + setDerivedParameters(); + } + + template + Repca::Repca(BusT* bus, const ModelDataT& data) + : bus_(bus), + monitor_(std::make_unique(data)) + { + initializeParameters(data); + initializeMonitor(); + size_ = static_cast(RepcaInternalVariables::MAXIMUM); + setDerivedParameters(); + } + + template + Repca::~Repca() + { + } + + template + void Repca::initializeParameters(const ModelDataT& data) + { + using Params = typename ModelDataT::Parameters; + + parameter_error_count_ = 0; + + auto load_real = [&](auto key, RealT& target, const char* name) + { + if (!data.parameters.contains(key)) + { + return; + } + + const auto& value = data.parameters.at(key); + 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() << "Repca: parameter '" << name << "' must be numeric\n"; + ++parameter_error_count_; + } + }; + + auto load_switch = [&](auto key, bool& target, const char* name) + { + if (!data.parameters.contains(key)) + { + return; + } + + const auto& value = data.parameters.at(key); + if (const auto* bool_value = std::get_if(&value)) + { + target = *bool_value; + } + else if (const auto* index_value = std::get_if(&value); + index_value && (*index_value == 0 || *index_value == 1)) + { + target = (*index_value == 1); + } + else if (const auto* real_value = std::get_if(&value); + real_value && (*real_value == ZERO || *real_value == ONE) ) + { + target = (*real_value == ONE); + } + else + { + Log::error() << "Repca: parameter '" << name << "' must be bool or 0/1\n"; + ++parameter_error_count_; + } + }; + + load_real(Params::mva, mva_, "mva"); + load_switch(Params::VcompFlag, VcompFlag_, "VcompFlag"); + load_switch(Params::RefFlag, RefFlag_, "RefFlag"); + load_switch(Params::Freqflag, Freqflag_, "Freqflag"); + load_real(Params::Tfltr, Tfltr_, "Tfltr"); + load_real(Params::Vfrz, Vfrz_, "Vfrz"); + load_real(Params::Rc, Rc_, "Rc"); + load_real(Params::Xc, Xc_, "Xc"); + load_real(Params::Kc, Kc_, "Kc"); + load_real(Params::dbdlow, dbdlow_, "dbdlow"); + load_real(Params::dbdupper, dbdupper_, "dbdupper"); + load_real(Params::emax, emax_, "emax"); + load_real(Params::emin, emin_, "emin"); + load_real(Params::Kp, Kp_, "Kp"); + load_real(Params::Ki, Ki_, "Ki"); + load_real(Params::Qmax, Qmax_, "Qmax"); + load_real(Params::Qmin, Qmin_, "Qmin"); + load_real(Params::Tft, Tft_, "Tft"); + load_real(Params::Tfv, Tfv_, "Tfv"); + load_real(Params::Tp, Tp_, "Tp"); + load_real(Params::fdbd1, fdbd1_, "fdbd1"); + load_real(Params::fdbd2, fdbd2_, "fdbd2"); + load_real(Params::Ddn, Ddn_, "Ddn"); + load_real(Params::Dup, Dup_, "Dup"); + load_real(Params::femax, femax_, "femax"); + load_real(Params::femin, femin_, "femin"); + load_real(Params::Kpg, Kpg_, "Kpg"); + load_real(Params::Kig, Kig_, "Kig"); + load_real(Params::Pmax, Pmax_, "Pmax"); + load_real(Params::Pmin, Pmin_, "Pmin"); + load_real(Params::Tlag, Tlag_, "Tlag"); + } + + template + void Repca::setDerivedParameters() + { + Tfltr_ = std::max(Tfltr_, TIME_CONSTANT_MINIMUM); + Tp_ = std::max(Tp_, TIME_CONSTANT_MINIMUM); + Tlag_ = std::max(Tlag_, TIME_CONSTANT_MINIMUM); + + va_component_base_ = mva_ * static_cast(1.0e6); + + vcomp_on_ = VcompFlag_ ? ONE : ZERO; + vcomp_off_ = ONE - vcomp_on_; + ref_on_ = RefFlag_ ? ONE : ZERO; + ref_off_ = ONE - ref_on_; + freq_on_ = Freqflag_ ? ONE : ZERO; + } + + template + scalar_type Repca::toComponentBase(scalar_type value) const + { + return value * va_system_base_ / va_component_base_; + } + + template + scalar_type Repca::toSystemBase(scalar_type value) const + { + return value / toComponentBase(static_cast(ONE)); + } + + template + const Model::VariableMonitorBase* Repca::getMonitor() const + { + return monitor_.get(); + } + + template + void Repca::initializeMonitor() + { + using Variable = typename ModelDataT::MonitorableVariables; + auto index = [](RepcaInternalVariables variable) + { + return static_cast(variable); + }; + + monitor_->set(Variable::qext, [this, index] + { return y_[index(RepcaInternalVariables::QEXT)]; }); + monitor_->set(Variable::pext, [this, index] + { return y_[index(RepcaInternalVariables::PEXT)]; }); + monitor_->set(Variable::vmeas, [this, index] + { return y_[index(RepcaInternalVariables::VMEAS)]; }); + monitor_->set(Variable::qmeas, [this, index] + { return y_[index(RepcaInternalVariables::QMEAS)]; }); + monitor_->set(Variable::pmeas, [this, index] + { return y_[index(RepcaInternalVariables::PMEAS)]; }); + } + + template + int Repca::setGridKitComponentID(IdxT component_id) + { + gridkit_component_id_ = component_id; + return 0; + } + + template + int Repca::allocate() + { + size_ = static_cast(RepcaInternalVariables::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); + abs_tol_.assign(size, ScalarT{0}); + variable_indices_.resize(size); + residual_indices_.resize(size); + + wb_.assign(2, ScalarT{0}); + + auto signal_size = static_cast(RepcaExternalVariables::MAXIMUM); + ws_.assign(signal_size, ScalarT{0}); + ws_indices_.assign(signal_size, INVALID_INDEX); + + for (IdxT j = 0; j < size_; ++j) + { + this->setVariableIndex(j, j); + this->setResidualIndex(j, j); + } + + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y_[static_cast(RepcaInternalVariables::QEXT)], + &(this->getVariableIndex(static_cast(RepcaInternalVariables::QEXT)))); + } + + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y_[static_cast(RepcaInternalVariables::PEXT)], + &(this->getVariableIndex(static_cast(RepcaInternalVariables::PEXT)))); + } + + return 0; + } + + template + int Repca::verify() const + { + int ret = static_cast(parameter_error_count_); + + auto check = [&](bool condition, const char* message) + { + if (!condition) + { + Log::error() << "Repca: " << message << '\n'; + ret += 1; + } + }; + + if (bus_ == nullptr) + { + Log::error() << "Repca: bus pointer is null\n"; + ret += 1; + } + + check(mva_ > ZERO, "mva must be positive"); + check(Tfv_ > ZERO, "Tfv must be positive"); + check(dbdlow_ <= ZERO && ZERO <= dbdupper_, + "dbdlow <= 0 <= dbdupper is required"); + check(emin_ <= ZERO && ZERO <= emax_, + "emin <= 0 <= emax is required"); + check(Qmin_ <= Qmax_, "Qmin must be less than or equal to Qmax"); + check(fdbd1_ <= ZERO && ZERO <= fdbd2_, + "fdbd1 <= 0 <= fdbd2 is required"); + check(femin_ <= ZERO && ZERO <= femax_, + "femin <= 0 <= femax is required"); + check(Pmin_ <= Pmax_, "Pmin must be less than or equal to Pmax"); + + auto check_required_signal = [&](bool attached, bool linked, const char* name) + { + if (!attached) + { + Log::error() << "Repca: " << name << " signal is required\n"; + ret += 1; + } + else if (!linked) + { + Log::error() << "Repca: " << name << " signal attached with no linked source\n"; + ret += 1; + } + }; + + check_required_signal( + signals_.template isAttached(), + signals_.template isAttached() + && signals_.template isLinked(), + "ibranchr"); + check_required_signal( + signals_.template isAttached(), + signals_.template isAttached() + && signals_.template isLinked(), + "ibranchi"); + check_required_signal( + signals_.template isAttached(), + signals_.template isAttached() + && signals_.template isLinked(), + "pbranch"); + check_required_signal( + signals_.template isAttached(), + signals_.template isAttached() + && signals_.template isLinked(), + "qbranch"); + check_required_signal( + signals_.template isAttached(), + signals_.template isAttached() + && signals_.template isLinked(), + "freq"); + + auto check_optional_signal = [&](bool attached, bool linked, const char* name) + { + if (attached && !linked) + { + Log::error() << "Repca: " << name << " signal attached with no linked source\n"; + ret += 1; + } + }; + + check_optional_signal( + signals_.template isAttached(), + signals_.template isAttached() + && signals_.template isLinked(), + "freqref"); + check_optional_signal( + signals_.template isAttached(), + signals_.template isAttached() + && signals_.template isLinked(), + "vref"); + check_optional_signal( + signals_.template isAttached(), + signals_.template isAttached() + && signals_.template isLinked(), + "qref"); + check_optional_signal( + signals_.template isAttached(), + signals_.template isAttached() + && signals_.template isLinked(), + "pplantref"); + + return ret; + } + + template + int Repca::initialize() + { + if (verify() > 0) + { + Log::error() << "Repca: cannot initialize with invalid configuration\n"; + return 1; + } + + const auto VMEAS = static_cast(RepcaInternalVariables::VMEAS); + const auto QMEAS = static_cast(RepcaInternalVariables::QMEAS); + const auto XQPI = static_cast(RepcaInternalVariables::XQPI); + const auto XQLAG = static_cast(RepcaInternalVariables::XQLAG); + const auto PMEAS = static_cast(RepcaInternalVariables::PMEAS); + const auto XPPI = static_cast(RepcaInternalVariables::XPPI); + const auto PREF = static_cast(RepcaInternalVariables::PREF); + const auto V = static_cast(RepcaInternalVariables::V); + const auto VLDC = static_cast(RepcaInternalVariables::VLDC); + const auto VDROOP = static_cast(RepcaInternalVariables::VDROOP); + const auto VCTRL = static_cast(RepcaInternalVariables::VCTRL); + const auto SFRZ = static_cast(RepcaInternalVariables::SFRZ); + const auto ERQ = static_cast(RepcaInternalVariables::ERQ); + const auto ERQDB = static_cast(RepcaInternalVariables::ERQDB); + const auto ERQLIM = static_cast(RepcaInternalVariables::ERQLIM); + const auto QPI = static_cast(RepcaInternalVariables::QPI); + const auto QEXT = static_cast(RepcaInternalVariables::QEXT); + const auto EF = static_cast(RepcaInternalVariables::EF); + const auto EP = static_cast(RepcaInternalVariables::EP); + const auto EPLIM = static_cast(RepcaInternalVariables::EPLIM); + const auto PPI = static_cast(RepcaInternalVariables::PPI); + const auto PEXT = static_cast(RepcaInternalVariables::PEXT); + + const ScalarT vr = Vr(); + const ScalarT vi = Vi(); + const ScalarT ibranchr = + signals_.template readExternalVariable(); + const ScalarT ibranchi = + signals_.template readExternalVariable(); + const ScalarT pbranch = + signals_.template readExternalVariable(); + const ScalarT qbranch = + signals_.template readExternalVariable(); + const ScalarT freq = + signals_.template readExternalVariable(); + + const ScalarT qext0 = toComponentBase(y_[QEXT]); + const ScalarT pext0 = toComponentBase(y_[PEXT]); + + const ScalarT vldc_r = vr - Rc_ * ibranchr + Xc_ * ibranchi; + const ScalarT vldc_i = vi - Rc_ * ibranchi - Xc_ * ibranchr; + + y_[V] = std::sqrt(vr * vr + vi * vi); + y_[VLDC] = std::sqrt(vldc_r * vldc_r + vldc_i * vldc_i); + y_[VDROOP] = y_[V] + Kc_ * toComponentBase(qbranch); + y_[VCTRL] = vcomp_on_ * y_[VLDC] + vcomp_off_ * y_[VDROOP]; + + y_[EF] = ZERO; + const ScalarT pfreq0 = Ddn_ * Math::ramp(y_[EF]) - Dup_ * Math::ramp(-y_[EF]); + + y_[VMEAS] = y_[VCTRL]; + y_[QMEAS] = toComponentBase(qbranch); + y_[PMEAS] = toComponentBase(pbranch); + y_[SFRZ] = Math::above(y_[V], Vfrz_); + + y_[ERQ] = ZERO; + y_[ERQDB] = Math::deadband2(y_[ERQ], dbdlow_, dbdupper_); + y_[ERQLIM] = Math::clamp(y_[ERQDB], emin_, emax_); + y_[QPI] = qext0; + y_[XQLAG] = qext0; + y_[QEXT] = toSystemBase(qext0); + y_[XQPI] = qext0 - Kp_ * y_[ERQLIM]; + + y_[EP] = ZERO; + y_[EPLIM] = Math::clamp(y_[EP], femin_, femax_); + y_[PREF] = Freqflag_ ? pext0 : Math::clamp(y_[PMEAS], Pmin_, Pmax_); + y_[PPI] = y_[PREF]; + y_[XPPI] = y_[PREF] - Kpg_ * y_[EPLIM]; + y_[PEXT] = toSystemBase(freq_on_ * y_[PREF]); + + const ScalarT q_aw = Math::antiwindup(y_[QPI], Ki_ * y_[ERQLIM], Qmin_, Qmax_); + const ScalarT p_aw = Math::antiwindup(y_[PPI], Kig_ * y_[EPLIM], Pmin_, Pmax_); + if (std::abs(static_cast(q_aw)) > static_cast(1.0e-10)) + { + Log::error() << "Repca: reactive PI antiwindup rate is nonzero at initialization\n"; + return 1; + } + if (std::abs(static_cast(p_aw)) > static_cast(1.0e-10)) + { + Log::error() << "Repca: active PI antiwindup rate is nonzero at initialization\n"; + return 1; + } + + freqref_set_ = freq; + vref_set_ = y_[VMEAS]; + qref_set_ = toSystemBase(y_[QMEAS]); + pplantref_set_ = toSystemBase(y_[PMEAS] - pfreq0); + + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable( + freqref_set_); + } + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(vref_set_); + } + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(qref_set_); + } + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable( + pplantref_set_); + } + + std::fill(yp_.begin(), yp_.end(), ZERO); + return 0; + } + + template + int Repca::tagDifferentiable() + { + std::fill(tag_.begin(), tag_.end(), false); + tag_[static_cast(RepcaInternalVariables::VMEAS)] = true; + tag_[static_cast(RepcaInternalVariables::QMEAS)] = true; + tag_[static_cast(RepcaInternalVariables::XQPI)] = true; + tag_[static_cast(RepcaInternalVariables::XQLAG)] = true; + tag_[static_cast(RepcaInternalVariables::PMEAS)] = true; + tag_[static_cast(RepcaInternalVariables::XPPI)] = true; + tag_[static_cast(RepcaInternalVariables::PREF)] = true; + return 0; + } + + template + int Repca::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + + template + __attribute__((always_inline)) inline int Repca::evaluateInternalResidual( + const ScalarT* y, + const ScalarT* yp, + const ScalarT* wb, + const ScalarT* ws, + ScalarT* f) + { + const auto VMEAS = static_cast(RepcaInternalVariables::VMEAS); + const auto QMEAS = static_cast(RepcaInternalVariables::QMEAS); + const auto XQPI = static_cast(RepcaInternalVariables::XQPI); + const auto XQLAG = static_cast(RepcaInternalVariables::XQLAG); + const auto PMEAS = static_cast(RepcaInternalVariables::PMEAS); + const auto XPPI = static_cast(RepcaInternalVariables::XPPI); + const auto PREF = static_cast(RepcaInternalVariables::PREF); + const auto V = static_cast(RepcaInternalVariables::V); + const auto VLDC = static_cast(RepcaInternalVariables::VLDC); + const auto VDROOP = static_cast(RepcaInternalVariables::VDROOP); + const auto VCTRL = static_cast(RepcaInternalVariables::VCTRL); + const auto SFRZ = static_cast(RepcaInternalVariables::SFRZ); + const auto ERQ = static_cast(RepcaInternalVariables::ERQ); + const auto ERQDB = static_cast(RepcaInternalVariables::ERQDB); + const auto ERQLIM = static_cast(RepcaInternalVariables::ERQLIM); + const auto QPI = static_cast(RepcaInternalVariables::QPI); + const auto QEXT = static_cast(RepcaInternalVariables::QEXT); + const auto EF = static_cast(RepcaInternalVariables::EF); + const auto EP = static_cast(RepcaInternalVariables::EP); + const auto EPLIM = static_cast(RepcaInternalVariables::EPLIM); + const auto PPI = static_cast(RepcaInternalVariables::PPI); + const auto PEXT = static_cast(RepcaInternalVariables::PEXT); + + const auto IBRANCHR = static_cast(RepcaExternalVariables::IBRANCHR); + const auto IBRANCHI = static_cast(RepcaExternalVariables::IBRANCHI); + const auto PBRANCH = static_cast(RepcaExternalVariables::PBRANCH); + const auto QBRANCH = static_cast(RepcaExternalVariables::QBRANCH); + const auto FREQ = static_cast(RepcaExternalVariables::FREQ); + const auto FREQREF = static_cast(RepcaExternalVariables::FREQREF); + const auto VREF = static_cast(RepcaExternalVariables::VREF); + const auto QREF = static_cast(RepcaExternalVariables::QREF); + const auto PPLANTREF = static_cast(RepcaExternalVariables::PPLANTREF); + + const ScalarT vmeas = y[VMEAS]; + const ScalarT qmeas = y[QMEAS]; + const ScalarT xqpi = y[XQPI]; + const ScalarT xqlag = y[XQLAG]; + const ScalarT pmeas = y[PMEAS]; + const ScalarT xppi = y[XPPI]; + const ScalarT pref = y[PREF]; + const ScalarT v = y[V]; + const ScalarT vldc = y[VLDC]; + const ScalarT vdroop = y[VDROOP]; + const ScalarT vctrl = y[VCTRL]; + const ScalarT sfrz = y[SFRZ]; + const ScalarT erq = y[ERQ]; + const ScalarT erqdb = y[ERQDB]; + const ScalarT erqlim = y[ERQLIM]; + const ScalarT qpi = y[QPI]; + const ScalarT qext = toComponentBase(y[QEXT]); + const ScalarT ef = y[EF]; + const ScalarT ep = y[EP]; + const ScalarT eplim = y[EPLIM]; + const ScalarT ppi = y[PPI]; + const ScalarT pext = toComponentBase(y[PEXT]); + + const ScalarT vmeas_dot = yp[VMEAS]; + const ScalarT qmeas_dot = yp[QMEAS]; + const ScalarT xqpi_dot = yp[XQPI]; + const ScalarT xqlag_dot = yp[XQLAG]; + const ScalarT pmeas_dot = yp[PMEAS]; + const ScalarT xppi_dot = yp[XPPI]; + const ScalarT pref_dot = yp[PREF]; + + const ScalarT vr = wb[0]; + const ScalarT vi = wb[1]; + + const ScalarT ibranchr = ws[IBRANCHR]; + const ScalarT ibranchi = ws[IBRANCHI]; + const ScalarT pbranch = ws[PBRANCH]; + const ScalarT qbranch = ws[QBRANCH]; + const ScalarT freq = ws[FREQ]; + const ScalarT freqref = ws[FREQREF]; + const ScalarT vref = ws[VREF]; + const ScalarT qref = toComponentBase(ws[QREF]); + const ScalarT pplantref = toComponentBase(ws[PPLANTREF]); + + const ScalarT vldc_r = vr - Rc_ * ibranchr + Xc_ * ibranchi; + const ScalarT vldc_i = vi - Rc_ * ibranchi - Xc_ * ibranchr; + const ScalarT pfreq = Ddn_ * Math::ramp(ef) - Dup_ * Math::ramp(-ef); + + f[VMEAS] = -vmeas_dot + (vctrl - vmeas) / Tfltr_; + f[QMEAS] = -qmeas_dot + (toComponentBase(qbranch) - qmeas) / Tfltr_; + f[XQPI] = -xqpi_dot + sfrz * Math::antiwindup(qpi, Ki_ * erqlim, Qmin_, Qmax_); + f[XQLAG] = -xqlag_dot + (qpi - xqlag) / Tfv_; + f[PMEAS] = -pmeas_dot + (toComponentBase(pbranch) - pmeas) / Tp_; + f[XPPI] = -xppi_dot + Math::antiwindup(ppi, Kig_ * eplim, Pmin_, Pmax_); + f[PREF] = -pref_dot + (ppi - pref) / Tlag_; + + f[V] = -v * v + vr * vr + vi * vi; + f[VLDC] = -vldc * vldc + vldc_r * vldc_r + vldc_i * vldc_i; + f[VDROOP] = -vdroop + v + Kc_ * toComponentBase(qbranch); + f[VCTRL] = -vctrl + vcomp_on_ * vldc + vcomp_off_ * vdroop; + f[SFRZ] = -sfrz + Math::above(v, Vfrz_); + f[ERQ] = -erq + ref_on_ * (vref - vmeas) + ref_off_ * (qref - qmeas); + f[ERQDB] = -erqdb + Math::deadband2(erq, dbdlow_, dbdupper_); + f[ERQLIM] = -erqlim + Math::clamp(erqdb, emin_, emax_); + f[QPI] = -qpi + Math::clamp(Kp_ * erqlim + xqpi, Qmin_, Qmax_); + f[QEXT] = -Tfv_ * (qext - xqlag) + Tft_ * (qpi - xqlag); + + f[EF] = -ef + Math::deadband2(freqref - freq, fdbd1_, fdbd2_); + f[EP] = -ep + pplantref - pmeas + pfreq; + f[EPLIM] = -eplim + Math::clamp(ep, femin_, femax_); + f[PPI] = -ppi + Math::clamp(Kpg_ * eplim + xppi, Pmin_, Pmax_); + f[PEXT] = -pext + freq_on_ * pref; + + return 0; + } + + template + int Repca::evaluateResidual() + { + const auto IBRANCHR = static_cast(RepcaExternalVariables::IBRANCHR); + const auto IBRANCHI = static_cast(RepcaExternalVariables::IBRANCHI); + const auto PBRANCH = static_cast(RepcaExternalVariables::PBRANCH); + const auto QBRANCH = static_cast(RepcaExternalVariables::QBRANCH); + const auto FREQ = static_cast(RepcaExternalVariables::FREQ); + const auto FREQREF = static_cast(RepcaExternalVariables::FREQREF); + const auto VREF = static_cast(RepcaExternalVariables::VREF); + const auto QREF = static_cast(RepcaExternalVariables::QREF); + const auto PPLANTREF = static_cast(RepcaExternalVariables::PPLANTREF); + + std::fill(ws_.begin(), ws_.end(), ZERO); + std::fill(ws_indices_.begin(), ws_indices_.end(), INVALID_INDEX); + + ws_[FREQREF] = freqref_set_; + ws_[VREF] = vref_set_; + ws_[QREF] = qref_set_; + ws_[PPLANTREF] = pplantref_set_; + + if (signals_.template isAttached()) + { + ws_[IBRANCHR] = + signals_.template readExternalVariable(); + ws_indices_[IBRANCHR] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[IBRANCHI] = + signals_.template readExternalVariable(); + ws_indices_[IBRANCHI] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[PBRANCH] = + signals_.template readExternalVariable(); + ws_indices_[PBRANCH] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[QBRANCH] = + signals_.template readExternalVariable(); + ws_indices_[QBRANCH] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[FREQ] = signals_.template readExternalVariable(); + ws_indices_[FREQ] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[FREQREF] = + signals_.template readExternalVariable(); + ws_indices_[FREQREF] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[VREF] = signals_.template readExternalVariable(); + ws_indices_[VREF] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[QREF] = signals_.template readExternalVariable(); + ws_indices_[QREF] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[PPLANTREF] = + signals_.template readExternalVariable(); + ws_indices_[PPLANTREF] = + signals_.template readExternalVariableIndex(); + } + + wb_[0] = Vr(); + wb_[1] = Vi(); + + evaluateInternalResidual(y_.data(), yp_.data(), wb_.data(), ws_.data(), f_.data()); + return 0; + } + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md index 5bff41ad4..a7837c5e8 100644 --- a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md +++ b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md @@ -144,16 +144,17 @@ are specified: `Genrou` | 6th order machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Tdop`, `Tdopp`, `Tqop`, `Tqopp`, `Xd`, `Xdp`, `Xdpp`, `Xq`, `Xqp`, `Xqpp`, `Xl`, `S10`, `S12`, `mva` | `ir`, `ii`, `p`, `q`, `delta`, `omega`, `speed` `Gensal` | 5th order salient-pole machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Tdop`, `Tdopp`, `Tqopp`, `Xd`, `Xdp`, `Xdpp`, `Xq`, `Xl`, `S10`, `S12`, `mva` | `ir`, `ii`, `p`, `q`, `delta`, `omega`, `speed`, `Eqp`, `psidp`, `psiqpp`, `psidpp`, `vd`, `vq`, `te`, `id`, `iq` `GenClassical` | the classical machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Xdp`, `mva` | `ir`, `ii`, `p`, `q`, `delta`, `omega` + `Repca` | the REPCA renewable plant-control model | `bus`, `ibranchr`, `ibranchi`, `qbranch`, `pbranch`, `vref`\*, `qref`\*, `pplantref`\*, `freq`\*, `freqref`\*, `qext`, `pext` | `mva`, `VcompFlag`, `RefFlag`, `Freqflag`, `Tfltr`, `Tft`, `Tfv`, `Tp`, `Tlag`, `Vfrz`, `Rc`, `Xc`, `Kc`, `dbdlow`, `dbdupper`, `emax`, `emin`, `Kp`, `Ki`, `Qmax`, `Qmin`, `fdbd1`, `fdbd2`, `Ddn`, `Dup`, `femax`, `femin`, `Kpg`, `Kig`, `Pmax`, `Pmin` | `qext`, `pext`, `vmeas`, `qmeas`, `pmeas`, `pref`, `vctrl`, `sfrz`, `qpi`, `pfreq`, `ppi` `Tgov1 ` | the TGOV1 governor model | `pmech`, `speed` | `Trate`, `R`, `T1`, `T2`, `T3`, `Pvmax`, `Pvmin`, `Dt` | `none` `Ieeet1` | the IEEET1 exciter model | `bus`, `speed`, `efd`, `vs`\* | `Tr`, `Ka`, `Ta`, `Ke`, `Te`, `Kf`, `Tf`, `Vrmin`, `Vrmax`, `E1`, `E2`, `Se1`, `Se2`, `Ispdlim` | `efd`, `ksat` `SexsPti` | the SEXS-PTI simplified exciter model | `bus`, `efd`, `vs`\* | `Ta`, `Tb`, `Te`, `K`, `Efdmax`, `Efdmin` | `efd` - `Ieeest` | the IEEEST stabilizer model | `input`, `output` | `A1`, `A2`, `A3`, `A4`, `A5`, `A6`, `T1`, `T2`, `T3`, `T4`, `T5`, `T6`, `Ks`, `Lsmin`, `Lsmax`, `Vcl`, `Vcu`, `Tdelay` | `vss` + `Ieeest` | the IEEEST stabilizer model | `input`, `output` | `A1`, `A2`, `A3`, `A4`, `A5`, `A6`, `T1`, `T2`, `T3`, `T4`, `T5`, `T6`, `Ks`, `Lsmin`, `Lsmax`, `Vcl`, `Vcu`, `Tdelay` | `vss` `BusFault` | simple impedance-based fault at a bus | `bus`, `status`\* | `state0`, `R`, `X` | `state`, `ir`, `ii` `BusToSignalAdapter` | signal adapter component for a bus | `bus`, `vr`, `vi`, `ir`, `ii` | | Ports marked with \* are optional and, if missing, will be assumed to be -connected to a constant value. This list is subject to change. - +connected to a constant value. +This list is subject to change. For `Branch`, `tap` and `phase` are optional parameters. If omitted, `tap` defaults to `1.0` and `phase` defaults to `0.0` radians. Bus `bus1` is the tap diff --git a/GridKit/Model/PhasorDynamics/SystemModelData.hpp b/GridKit/Model/PhasorDynamics/SystemModelData.hpp index 8a6e73056..230cc891a 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelData.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelData.hpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -41,6 +42,7 @@ namespace GridKit using BusDataT = BusData; using BusToSignalAdapterDataT = BusToSignalAdapterData; using BusFaultDataT = BusFaultData; + using RepcaDataT = Converter::RepcaData; using Tgov1DataT = Governor::Tgov1Data; using Ieeet1DataT = Exciter::Ieeet1Data; using SexsPtiDataT = Exciter::SexsPtiData; @@ -95,6 +97,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 repca; ///< REPCA plant controllers within the model std::vector genrou; ///< GENROU instances within the model std::vector gensal; ///< GENSAL instances within the model std::vector genclassical; ///< Classical generator instances within the model diff --git a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp index a3682feca..555e87ca7 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 == "Repca") + { + typename SystemModelData::RepcaDataT repca; + raw_component.get_to(repca); + sm.repca.push_back(repca); + } else if (kind == "Tgov1") { typename SystemModelData::Tgov1DataT gov; diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index 287b328dc..c364ad806 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -42,6 +42,7 @@ namespace GridKit using namespace Governor; using namespace Exciter; using namespace Stabilizer; + using namespace Converter; owns_components_ = true; @@ -106,6 +107,87 @@ namespace GridKit addComponent(adapter); } + // Add REPCA plant controllers + for (const auto& repcadata : data.repca) + { + IdxT bus_index = 0; + if (repcadata.ports.contains(RepcaPorts::bus)) + { + bus_index = repcadata.ports.at(RepcaPorts::bus); + } + + auto* repca = new Repca(getBus(bus_index), repcadata); + + if (repcadata.ports.contains(RepcaPorts::ibranchr)) + { + const IdxT signal = repcadata.ports.at(RepcaPorts::ibranchr); + constexpr auto IBRANCHR = RepcaExternalVariables::IBRANCHR; + repca->getSignals().template attachSignalNode(getSignal(signal)); + } + if (repcadata.ports.contains(RepcaPorts::ibranchi)) + { + const IdxT signal = repcadata.ports.at(RepcaPorts::ibranchi); + constexpr auto IBRANCHI = RepcaExternalVariables::IBRANCHI; + repca->getSignals().template attachSignalNode(getSignal(signal)); + } + if (repcadata.ports.contains(RepcaPorts::qbranch)) + { + const IdxT signal = repcadata.ports.at(RepcaPorts::qbranch); + constexpr auto QBRANCH = RepcaExternalVariables::QBRANCH; + repca->getSignals().template attachSignalNode(getSignal(signal)); + } + if (repcadata.ports.contains(RepcaPorts::pbranch)) + { + const IdxT signal = repcadata.ports.at(RepcaPorts::pbranch); + constexpr auto PBRANCH = RepcaExternalVariables::PBRANCH; + repca->getSignals().template attachSignalNode(getSignal(signal)); + } + if (repcadata.ports.contains(RepcaPorts::vref)) + { + const IdxT signal = repcadata.ports.at(RepcaPorts::vref); + constexpr auto VREF = RepcaExternalVariables::VREF; + repca->getSignals().template attachSignalNode(getSignal(signal)); + } + if (repcadata.ports.contains(RepcaPorts::qref)) + { + const IdxT signal = repcadata.ports.at(RepcaPorts::qref); + constexpr auto QREF = RepcaExternalVariables::QREF; + repca->getSignals().template attachSignalNode(getSignal(signal)); + } + if (repcadata.ports.contains(RepcaPorts::pplantref)) + { + const IdxT signal = repcadata.ports.at(RepcaPorts::pplantref); + constexpr auto PPLANTREF = RepcaExternalVariables::PPLANTREF; + repca->getSignals().template attachSignalNode(getSignal(signal)); + } + if (repcadata.ports.contains(RepcaPorts::freq)) + { + const IdxT signal = repcadata.ports.at(RepcaPorts::freq); + constexpr auto FREQ = RepcaExternalVariables::FREQ; + repca->getSignals().template attachSignalNode(getSignal(signal)); + } + if (repcadata.ports.contains(RepcaPorts::freqref)) + { + const IdxT signal = repcadata.ports.at(RepcaPorts::freqref); + constexpr auto FREQREF = RepcaExternalVariables::FREQREF; + repca->getSignals().template attachSignalNode(getSignal(signal)); + } + if (repcadata.ports.contains(RepcaPorts::qext)) + { + const IdxT signal = repcadata.ports.at(RepcaPorts::qext); + constexpr auto QEXT = RepcaInternalVariables::QEXT; + repca->getSignals().template assignSignalNode(getSignal(signal)); + } + if (repcadata.ports.contains(RepcaPorts::pext)) + { + const IdxT signal = repcadata.ports.at(RepcaPorts::pext); + constexpr auto PEXT = RepcaInternalVariables::PEXT; + repca->getSignals().template assignSignalNode(getSignal(signal)); + } + + addComponent(repca); + } + // Add branches for (const auto& branchdata : data.branch) { diff --git a/docs/Figures/PhasorDynamics_REPCA_Diagram.png b/docs/Figures/PhasorDynamics_REPCA_Diagram.png new file mode 100644 index 000000000..08fef3c15 Binary files /dev/null and b/docs/Figures/PhasorDynamics_REPCA_Diagram.png differ diff --git a/docs/GridKit/Model/PhasorDynamics/Converter/README.md b/docs/GridKit/Model/PhasorDynamics/Converter/README.md index fa56e8dd8..857859edf 100644 --- a/docs/GridKit/Model/PhasorDynamics/Converter/README.md +++ b/docs/GridKit/Model/PhasorDynamics/Converter/README.md @@ -8,6 +8,7 @@ REGCA REGCB REECA +REPCA ``` ```{include} ../../../../../GridKit/Model/PhasorDynamics/Converter/README.md diff --git a/docs/GridKit/Model/PhasorDynamics/Converter/REPCA/README.md b/docs/GridKit/Model/PhasorDynamics/Converter/REPCA/README.md new file mode 100644 index 000000000..855a3ea10 --- /dev/null +++ b/docs/GridKit/Model/PhasorDynamics/Converter/REPCA/README.md @@ -0,0 +1,6 @@ +# REPCA + +```{include} ../../../../../../GridKit/Model/PhasorDynamics/Converter/REPCA/README.md +:start-line: 1 +:relative-images: +``` 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 3d2973b62..cc7c11cc8 100644 --- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt +++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt @@ -104,6 +104,14 @@ target_link_libraries( GridKit::phasor_dynamics_systemmodel_dependency_tracking GridKit::testing) +add_executable(test_phasor_converter_repca runConverterRepcaTests.cpp) +target_link_libraries( + test_phasor_converter_repca + 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 @@ -149,6 +157,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 PhasorDynamicsConverterRepcaTest COMMAND test_phasor_converter_repca) add_test(NAME PhasorDynamicsStabilizerIeeestTest COMMAND test_phasor_stabilizer_ieeest) add_test(NAME PhasorDynamicsGenClassicalTest COMMAND test_phasor_gen_classical) add_test(NAME PhasorDynamicsLoadZTest COMMAND test_phasor_loadz) @@ -171,6 +180,7 @@ install( test_phasor_exciter_ieeet1 test_phasor_gensal test_phasor_exciter_sexspti + test_phasor_converter_repca test_phasor_stabilizer_ieeest test_phasor_gen_classical test_phasor_system diff --git a/tests/UnitTests/PhasorDynamics/ConverterRepcaTests.hpp b/tests/UnitTests/PhasorDynamics/ConverterRepcaTests.hpp new file mode 100644 index 000000000..0c90178ee --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/ConverterRepcaTests.hpp @@ -0,0 +1,847 @@ +#pragma once + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + template + class ConverterRepcaTests + { + public: + using ScalarT = scalar_type; + using IdxT = index_type; + using RealT = typename PhasorDynamics::Component::RealT; + + static constexpr ScalarT kTol = static_cast(1.0e-8); + + TestOutcome constructionAndValidation() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + ScalarT ir{0.0}; + ScalarT ii{0.0}; + ScalarT p{0.4}; + ScalarT q{0.1}; + ScalarT freq{1.0}; + IdxT ir_i = 10; + IdxT ii_i = 11; + IdxT p_i = 12; + IdxT q_i = 13; + IdxT f_i = 14; + + PhasorDynamics::SignalNode ir_node; + PhasorDynamics::SignalNode ii_node; + PhasorDynamics::SignalNode p_node; + PhasorDynamics::SignalNode q_node; + PhasorDynamics::SignalNode freq_node; + ir_node.set(&ir, &ir_i); + ii_node.set(&ii, &ii_i); + p_node.set(&p, &p_i); + q_node.set(&q, &q_i); + freq_node.set(&freq, &f_i); + + auto data = makeData(); + PhasorDynamics::Converter::Repca repca(&bus, data); + auto& signals = repca.getSignals(); + signals.template attachSignalNode(&ir_node); + signals.template attachSignalNode(&ii_node); + signals.template attachSignalNode(&p_node); + signals.template attachSignalNode(&q_node); + signals.template attachSignalNode(&freq_node); + + success *= (repca.size() == static_cast(Vars::MAXIMUM)); + success *= (repca.getMonitor() != nullptr); + success *= (repca.verify() == 0); + + auto bad = data; + bad.parameters[Params::mva] = static_cast(0.0); + PhasorDynamics::Converter::Repca bad_mva(&bus, bad); + bad_mva.getSignals().template attachSignalNode(&ir_node); + bad_mva.getSignals().template attachSignalNode(&ii_node); + bad_mva.getSignals().template attachSignalNode(&p_node); + bad_mva.getSignals().template attachSignalNode(&q_node); + bad_mva.getSignals().template attachSignalNode(&freq_node); + success *= (bad_mva.verify() > 0); + + bad = data; + bad.parameters[Params::Tfv] = static_cast(0.0); + PhasorDynamics::Converter::Repca bad_tfv(&bus, bad); + bad_tfv.getSignals().template attachSignalNode(&ir_node); + bad_tfv.getSignals().template attachSignalNode(&ii_node); + bad_tfv.getSignals().template attachSignalNode(&p_node); + bad_tfv.getSignals().template attachSignalNode(&q_node); + bad_tfv.getSignals().template attachSignalNode(&freq_node); + success *= (bad_tfv.verify() > 0); + + bad = data; + bad.parameters[Params::Qmin] = static_cast(2.0); + PhasorDynamics::Converter::Repca bad_q_limits(&bus, bad); + bad_q_limits.getSignals().template attachSignalNode(&ir_node); + bad_q_limits.getSignals().template attachSignalNode(&ii_node); + bad_q_limits.getSignals().template attachSignalNode(&p_node); + bad_q_limits.getSignals().template attachSignalNode(&q_node); + bad_q_limits.getSignals().template attachSignalNode(&freq_node); + success *= (bad_q_limits.verify() > 0); + + bad = data; + bad.parameters[Params::fdbd1] = static_cast(0.1); + PhasorDynamics::Converter::Repca bad_deadband(&bus, bad); + bad_deadband.getSignals().template attachSignalNode(&ir_node); + bad_deadband.getSignals().template attachSignalNode(&ii_node); + bad_deadband.getSignals().template attachSignalNode(&p_node); + bad_deadband.getSignals().template attachSignalNode(&q_node); + bad_deadband.getSignals().template attachSignalNode(&freq_node); + success *= (bad_deadband.verify() > 0); + + return success.report(__func__); + } + + TestOutcome signalVerification() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + PhasorDynamics::Converter::Repca repca(&bus, makeData()); + + PhasorDynamics::SignalNode ir_node; + repca.getSignals().template attachSignalNode(&ir_node); + success *= (repca.verify() > 0); + + ScalarT ir{0.0}; + ScalarT ii{0.0}; + ScalarT p{0.4}; + ScalarT q{0.1}; + ScalarT freq{1.0}; + ScalarT vref{1.0}; + IdxT ir_i = 20; + IdxT ii_i = 21; + IdxT p_i = 22; + IdxT q_i = 23; + IdxT f_i = 24; + IdxT v_i = 25; + + PhasorDynamics::SignalNode ii_node; + PhasorDynamics::SignalNode p_node; + PhasorDynamics::SignalNode q_node; + PhasorDynamics::SignalNode freq_node; + PhasorDynamics::SignalNode vref_node; + ir_node.set(&ir, &ir_i); + ii_node.set(&ii, &ii_i); + p_node.set(&p, &p_i); + q_node.set(&q, &q_i); + freq_node.set(&freq, &f_i); + + repca.getSignals().template attachSignalNode(&ii_node); + repca.getSignals().template attachSignalNode(&p_node); + repca.getSignals().template attachSignalNode(&q_node); + repca.getSignals().template attachSignalNode(&freq_node); + success *= (repca.verify() == 0); + + repca.getSignals().template attachSignalNode(&vref_node); + success *= (repca.verify() > 0); + vref_node.set(&vref, &v_i); + success *= (repca.verify() == 0); + + return success.report(__func__); + } + + TestOutcome initializationAndResidual() + { + TestStatus success = true; + + auto data = makeData(); + data.parameters[Params::mva] = static_cast(50.0); + data.parameters[Params::Freqflag] = true; + + PhasorDynamics::Bus bus(0.8, 0.6); + bus.allocate(); + bus.initialize(); + + ScalarT ir{0.2}; + ScalarT ii{-0.1}; + ScalarT p{0.4}; + ScalarT q{0.1}; + ScalarT freq{1.0}; + ScalarT freqref{-4.0}; + ScalarT vref{-3.0}; + ScalarT qref{-2.0}; + ScalarT pplantref{-1.0}; + IdxT ir_i = 30; + IdxT ii_i = 31; + IdxT p_i = 32; + IdxT q_i = 33; + IdxT f_i = 34; + IdxT freqref_i = 35; + IdxT vref_i = 36; + IdxT qref_i = 37; + IdxT plantref_i = 38; + + PhasorDynamics::SignalNode ir_node; + PhasorDynamics::SignalNode ii_node; + PhasorDynamics::SignalNode p_node; + PhasorDynamics::SignalNode q_node; + PhasorDynamics::SignalNode freq_node; + PhasorDynamics::SignalNode freqref_node; + PhasorDynamics::SignalNode vref_node; + PhasorDynamics::SignalNode qref_node; + PhasorDynamics::SignalNode plantref_node; + PhasorDynamics::SignalNode qext_node; + PhasorDynamics::SignalNode pext_node; + + ir_node.set(&ir, &ir_i); + ii_node.set(&ii, &ii_i); + p_node.set(&p, &p_i); + q_node.set(&q, &q_i); + freq_node.set(&freq, &f_i); + freqref_node.set(&freqref, &freqref_i); + vref_node.set(&vref, &vref_i); + qref_node.set(&qref, &qref_i); + plantref_node.set(&pplantref, &plantref_i); + + PhasorDynamics::Converter::Repca repca(&bus, data); + auto& signals = repca.getSignals(); + signals.template attachSignalNode(&ir_node); + signals.template attachSignalNode(&ii_node); + signals.template attachSignalNode(&p_node); + signals.template attachSignalNode(&q_node); + signals.template attachSignalNode(&freq_node); + signals.template attachSignalNode(&freqref_node); + signals.template attachSignalNode(&vref_node); + signals.template attachSignalNode(&qref_node); + signals.template attachSignalNode(&plantref_node); + signals.template assignSignalNode(&qext_node); + signals.template assignSignalNode(&pext_node); + + success *= (repca.allocate() == 0); + repca.y()[index(Vars::QEXT)] = static_cast(0.25); + repca.y()[index(Vars::PEXT)] = static_cast(0.70); + success *= (repca.initialize() == 0); + success *= (repca.tagDifferentiable() == 0); + success *= (repca.evaluateResidual() == 0); + + const ScalarT base_ratio = static_cast(2.0); + const ScalarT qmeas = base_ratio * q; + const ScalarT pmeas = base_ratio * p; + const ScalarT pfreq0 = + static_cast(20.0) * Math::ramp(static_cast(0.0)); + + success *= isEqual(repca.y()[index(Vars::QEXT)], static_cast(0.25), kTol); + success *= isEqual(repca.y()[index(Vars::PEXT)], static_cast(0.70), kTol); + success *= isEqual(qext_node.read(), static_cast(0.25), kTol); + success *= isEqual(pext_node.read(), static_cast(0.70), kTol); + success *= isEqual(freqref, freq, kTol); + success *= isEqual(vref, repca.y()[index(Vars::VMEAS)], kTol); + success *= isEqual(qref, qmeas / base_ratio, kTol); + success *= isEqual(pplantref, (pmeas - pfreq0) / base_ratio, kTol); + + for (size_t i = 0; i < repca.getResidual().size(); ++i) + { + if (!isEqual(repca.getResidual()[i], static_cast(0.0), kTol)) + { + std::cout << "REPCA residual row " << i << " is " + << repca.getResidual()[i] << "\n"; + success = false; + } + } + + return success.report(__func__); + } + + TestOutcome residualEquations() + { + TestStatus success = true; + + auto data = makeResidualData(); + PhasorDynamics::Bus bus(0.9, 0.4); + bus.allocate(); + bus.initialize(); + + ScalarT ir{0.08}; + ScalarT ii{-0.02}; + ScalarT p{0.41}; + ScalarT q{0.13}; + ScalarT freq{0.99}; + ScalarT freqref{1.0}; + ScalarT vref{1.01}; + ScalarT qref{0.12}; + ScalarT pplantref{0.55}; + IdxT ir_i = 40; + IdxT ii_i = 41; + IdxT p_i = 42; + IdxT q_i = 43; + IdxT f_i = 44; + IdxT freqref_i = 45; + IdxT vref_i = 46; + IdxT qref_i = 47; + IdxT plantref_i = 48; + + PhasorDynamics::SignalNode ir_node; + PhasorDynamics::SignalNode ii_node; + PhasorDynamics::SignalNode p_node; + PhasorDynamics::SignalNode q_node; + PhasorDynamics::SignalNode freq_node; + PhasorDynamics::SignalNode freqref_node; + PhasorDynamics::SignalNode vref_node; + PhasorDynamics::SignalNode qref_node; + PhasorDynamics::SignalNode plantref_node; + + ir_node.set(&ir, &ir_i); + ii_node.set(&ii, &ii_i); + p_node.set(&p, &p_i); + q_node.set(&q, &q_i); + freq_node.set(&freq, &f_i); + freqref_node.set(&freqref, &freqref_i); + vref_node.set(&vref, &vref_i); + qref_node.set(&qref, &qref_i); + plantref_node.set(&pplantref, &plantref_i); + + PhasorDynamics::Converter::Repca repca(&bus, data); + repca.getSignals().template attachSignalNode(&ir_node); + repca.getSignals().template attachSignalNode(&ii_node); + repca.getSignals().template attachSignalNode(&p_node); + repca.getSignals().template attachSignalNode(&q_node); + repca.getSignals().template attachSignalNode(&freq_node); + repca.getSignals().template attachSignalNode(&freqref_node); + repca.getSignals().template attachSignalNode(&vref_node); + repca.getSignals().template attachSignalNode(&qref_node); + repca.getSignals().template attachSignalNode(&plantref_node); + + success *= (repca.allocate() == 0); + + repca.y()[index(Vars::VMEAS)] = static_cast(0.98); + repca.y()[index(Vars::QMEAS)] = static_cast(0.11); + repca.y()[index(Vars::XQPI)] = static_cast(0.07); + repca.y()[index(Vars::XQLAG)] = static_cast(0.14); + repca.y()[index(Vars::PMEAS)] = static_cast(0.44); + repca.y()[index(Vars::XPPI)] = static_cast(0.21); + repca.y()[index(Vars::PREF)] = static_cast(0.60); + repca.y()[index(Vars::V)] = static_cast(1.0); + repca.y()[index(Vars::VLDC)] = static_cast(0.99); + repca.y()[index(Vars::VDROOP)] = static_cast(1.05); + repca.y()[index(Vars::VCTRL)] = static_cast(1.02); + repca.y()[index(Vars::SFRZ)] = static_cast(0.8); + repca.y()[index(Vars::ERQ)] = static_cast(0.03); + repca.y()[index(Vars::ERQDB)] = static_cast(0.02); + repca.y()[index(Vars::ERQLIM)] = static_cast(0.02); + repca.y()[index(Vars::QPI)] = static_cast(0.27); + repca.y()[index(Vars::QEXT)] = static_cast(0.20); + repca.y()[index(Vars::EF)] = static_cast(0.01); + repca.y()[index(Vars::EP)] = static_cast(0.04); + repca.y()[index(Vars::EPLIM)] = static_cast(0.04); + repca.y()[index(Vars::PPI)] = static_cast(0.66); + repca.y()[index(Vars::PEXT)] = static_cast(0.61); + + repca.yp()[index(Vars::VMEAS)] = static_cast(0.01); + repca.yp()[index(Vars::QMEAS)] = static_cast(-0.02); + repca.yp()[index(Vars::XQPI)] = static_cast(0.03); + repca.yp()[index(Vars::XQLAG)] = static_cast(-0.04); + repca.yp()[index(Vars::PMEAS)] = static_cast(0.02); + repca.yp()[index(Vars::XPPI)] = static_cast(-0.01); + repca.yp()[index(Vars::PREF)] = static_cast(0.05); + + success *= (repca.evaluateResidual() == 0); + + const ScalarT vldc_r = bus.Vr() - static_cast(0.02) * ir + + static_cast(0.03) * ii; + const ScalarT vldc_i = bus.Vi() - static_cast(0.02) * ii + - static_cast(0.03) * ir; + const ScalarT pfreq = static_cast(2.0) * Math::ramp(repca.y()[index(Vars::EF)]) + - Math::ramp(-repca.y()[index(Vars::EF)]); + + std::vector expected(static_cast(Vars::MAXIMUM), + static_cast(0.0)); + expected[index(Vars::VMEAS)] = + -repca.yp()[index(Vars::VMEAS)] + + (repca.y()[index(Vars::VCTRL)] - repca.y()[index(Vars::VMEAS)]) + / static_cast(0.2); + expected[index(Vars::QMEAS)] = + -repca.yp()[index(Vars::QMEAS)] + + (q - repca.y()[index(Vars::QMEAS)]) / static_cast(0.2); + expected[index(Vars::XQPI)] = + -repca.yp()[index(Vars::XQPI)] + + repca.y()[index(Vars::SFRZ)] + * Math::antiwindup(repca.y()[index(Vars::QPI)], + static_cast(3.0) + * repca.y()[index(Vars::ERQLIM)], + static_cast(-0.8), + static_cast(0.9)); + expected[index(Vars::XQLAG)] = + -repca.yp()[index(Vars::XQLAG)] + + (repca.y()[index(Vars::QPI)] - repca.y()[index(Vars::XQLAG)]) + / static_cast(1.5); + expected[index(Vars::PMEAS)] = + -repca.yp()[index(Vars::PMEAS)] + + (p - repca.y()[index(Vars::PMEAS)]) / static_cast(0.4); + expected[index(Vars::XPPI)] = + -repca.yp()[index(Vars::XPPI)] + + Math::antiwindup(repca.y()[index(Vars::PPI)], + static_cast(1.8) * repca.y()[index(Vars::EPLIM)], + static_cast(0.0), + static_cast(1.2)); + expected[index(Vars::PREF)] = + -repca.yp()[index(Vars::PREF)] + + (repca.y()[index(Vars::PPI)] - repca.y()[index(Vars::PREF)]) + / static_cast(0.5); + + expected[index(Vars::V)] = + -repca.y()[index(Vars::V)] * repca.y()[index(Vars::V)] + + bus.Vr() * bus.Vr() + bus.Vi() * bus.Vi(); + expected[index(Vars::VLDC)] = + -repca.y()[index(Vars::VLDC)] * repca.y()[index(Vars::VLDC)] + + vldc_r * vldc_r + vldc_i * vldc_i; + expected[index(Vars::VDROOP)] = + -repca.y()[index(Vars::VDROOP)] + repca.y()[index(Vars::V)] + + static_cast(0.4) * q; + expected[index(Vars::VCTRL)] = + -repca.y()[index(Vars::VCTRL)] + repca.y()[index(Vars::VLDC)]; + expected[index(Vars::SFRZ)] = + -repca.y()[index(Vars::SFRZ)] + + Math::above(repca.y()[index(Vars::V)], static_cast(0.7)); + expected[index(Vars::ERQ)] = + -repca.y()[index(Vars::ERQ)] + vref - repca.y()[index(Vars::VMEAS)]; + expected[index(Vars::ERQDB)] = + -repca.y()[index(Vars::ERQDB)] + + Math::deadband2(repca.y()[index(Vars::ERQ)], + static_cast(-0.02), + static_cast(0.03)); + expected[index(Vars::ERQLIM)] = + -repca.y()[index(Vars::ERQLIM)] + + Math::clamp(repca.y()[index(Vars::ERQDB)], + static_cast(-0.7), + static_cast(0.8)); + expected[index(Vars::QPI)] = + -repca.y()[index(Vars::QPI)] + + Math::clamp(static_cast(2.0) * repca.y()[index(Vars::ERQLIM)] + + repca.y()[index(Vars::XQPI)], + static_cast(-0.8), + static_cast(0.9)); + expected[index(Vars::QEXT)] = + -static_cast(1.5) + * (repca.y()[index(Vars::QEXT)] - repca.y()[index(Vars::XQLAG)]) + + static_cast(0.2) + * (repca.y()[index(Vars::QPI)] - repca.y()[index(Vars::XQLAG)]); + + expected[index(Vars::EF)] = + -repca.y()[index(Vars::EF)] + + Math::deadband2(freqref - freq, + static_cast(-0.01), + static_cast(0.015)); + expected[index(Vars::EP)] = + -repca.y()[index(Vars::EP)] + pplantref - repca.y()[index(Vars::PMEAS)] + + pfreq; + expected[index(Vars::EPLIM)] = + -repca.y()[index(Vars::EPLIM)] + + Math::clamp(repca.y()[index(Vars::EP)], + static_cast(-0.5), + static_cast(0.6)); + expected[index(Vars::PPI)] = + -repca.y()[index(Vars::PPI)] + + Math::clamp(static_cast(1.7) * repca.y()[index(Vars::EPLIM)] + + repca.y()[index(Vars::XPPI)], + static_cast(0.0), + static_cast(1.2)); + expected[index(Vars::PEXT)] = + -repca.y()[index(Vars::PEXT)] + repca.y()[index(Vars::PREF)]; + + for (size_t i = 0; i < expected.size(); ++i) + { + if (!isEqual(repca.getResidual()[i], expected[i], kTol)) + { + std::cout << "REPCA residual mismatch at row " << i << ": " + << repca.getResidual()[i] << " != " << expected[i] << "\n"; + success = false; + } + } + + return success.report(__func__); + } + + TestOutcome jsonParseAndSystemAssembly() + { + TestStatus success = true; + + std::istringstream input(R"json( +{ + "header": { + "format_version": 0, + "format_revision": 1, + "case_name": "REPCA parser", + "case_description": "REPCA 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 + } + ], + "signals": [ + { "signal_id": 10, "name": "Ir" }, + { "signal_id": 11, "name": "Ii" }, + { "signal_id": 12, "name": "Pbranch" }, + { "signal_id": 13, "name": "Qbranch" }, + { "signal_id": 14, "name": "Freq" }, + { "signal_id": 15, "name": "Qext" }, + { "signal_id": 16, "name": "Pext" } + ], + "devices": [ + { + "class": "Repca", + "ports": { + "bus": 1, + "ibranchr": 10, + "ibranchi": 11, + "pbranch": 12, + "qbranch": 13, + "freq": 14, + "qext": 15, + "pext": 16 + }, + "id": "RP1", + "params": { "mva": 100.0 }, + "mon": ["qext", "pext", "vmeas", "qmeas", "pmeas"] + } + ] +} +)json"); + + auto data = PhasorDynamics::parseSystemModelData(input); + success *= (data.repca.size() == 1); + success *= (data.repca[0].device_class == "Repca"); + success *= (data.repca[0].ports.at(PhasorDynamics::Converter::RepcaPorts::bus) == 1); + success *= (std::get_if( + &data.repca[0].parameters.at(Params::mva)) + != nullptr); + + PhasorDynamics::SystemModel system(data); + success *= (system.size() == 0); + + 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::RepcaParameters; + using Vars = PhasorDynamics::Converter::RepcaInternalVariables; + using Ext = PhasorDynamics::Converter::RepcaExternalVariables; + + static size_t index(Vars variable) + { + return static_cast(variable); + } + + auto makeData() -> PhasorDynamics::Converter::RepcaData + { + using Mon = PhasorDynamics::Converter::RepcaMonitorableVariables; + + PhasorDynamics::Converter::RepcaData data; + data.device_class = "Repca"; + data.disambiguation_string = "repca_test"; + data.monitored_variables.insert(Mon::qext); + data.monitored_variables.insert(Mon::pext); + data.monitored_variables.insert(Mon::vmeas); + data.monitored_variables.insert(Mon::qmeas); + data.monitored_variables.insert(Mon::pmeas); + return data; + } + + auto makeResidualData() -> PhasorDynamics::Converter::RepcaData + { + auto data = makeData(); + + data.parameters[Params::mva] = static_cast(100.0); + data.parameters[Params::VcompFlag] = true; + data.parameters[Params::RefFlag] = true; + data.parameters[Params::Freqflag] = true; + data.parameters[Params::Tfltr] = static_cast(0.2); + data.parameters[Params::Vfrz] = static_cast(0.7); + data.parameters[Params::Rc] = static_cast(0.02); + data.parameters[Params::Xc] = static_cast(0.03); + data.parameters[Params::Kc] = static_cast(0.4); + data.parameters[Params::dbdlow] = static_cast(-0.02); + data.parameters[Params::dbdupper] = static_cast(0.03); + data.parameters[Params::emax] = static_cast(0.8); + data.parameters[Params::emin] = static_cast(-0.7); + data.parameters[Params::Kp] = static_cast(2.0); + data.parameters[Params::Ki] = static_cast(3.0); + data.parameters[Params::Qmax] = static_cast(0.9); + data.parameters[Params::Qmin] = static_cast(-0.8); + data.parameters[Params::Tft] = static_cast(0.2); + data.parameters[Params::Tfv] = static_cast(1.5); + data.parameters[Params::Tp] = static_cast(0.4); + data.parameters[Params::fdbd1] = static_cast(-0.01); + data.parameters[Params::fdbd2] = static_cast(0.015); + data.parameters[Params::Ddn] = static_cast(2.0); + data.parameters[Params::Dup] = static_cast(1.0); + data.parameters[Params::femax] = static_cast(0.6); + data.parameters[Params::femin] = static_cast(-0.5); + data.parameters[Params::Kpg] = static_cast(1.7); + data.parameters[Params::Kig] = static_cast(1.8); + data.parameters[Params::Pmax] = static_cast(1.2); + data.parameters[Params::Pmin] = static_cast(0.0); + data.parameters[Params::Tlag] = static_cast(0.5); + return data; + } + +#ifdef GRIDKIT_ENABLE_ENZYME + void setValue(double& target, double value) + { + target = value; + } + + void setValue(DependencyTracking::Variable& target, double value) + { + target.setValue(value); + } + + template + void setJacobianState(model_type& repca, bus_type& bus) + { + setValue(bus.Vr(), 0.9); + setValue(bus.Vi(), 0.4); + + setValue(repca.y()[index(Vars::VMEAS)], 0.98); + setValue(repca.y()[index(Vars::QMEAS)], 0.11); + setValue(repca.y()[index(Vars::XQPI)], 0.07); + setValue(repca.y()[index(Vars::XQLAG)], 0.14); + setValue(repca.y()[index(Vars::PMEAS)], 0.44); + setValue(repca.y()[index(Vars::XPPI)], 0.21); + setValue(repca.y()[index(Vars::PREF)], 0.60); + setValue(repca.y()[index(Vars::V)], 1.0); + setValue(repca.y()[index(Vars::VLDC)], 0.99); + setValue(repca.y()[index(Vars::VDROOP)], 1.05); + setValue(repca.y()[index(Vars::VCTRL)], 1.02); + setValue(repca.y()[index(Vars::SFRZ)], 0.8); + setValue(repca.y()[index(Vars::ERQ)], 0.03); + setValue(repca.y()[index(Vars::ERQDB)], 0.02); + setValue(repca.y()[index(Vars::ERQLIM)], 0.02); + setValue(repca.y()[index(Vars::QPI)], 0.27); + setValue(repca.y()[index(Vars::QEXT)], 0.20); + setValue(repca.y()[index(Vars::EF)], 0.01); + setValue(repca.y()[index(Vars::EP)], 0.04); + setValue(repca.y()[index(Vars::EPLIM)], 0.04); + setValue(repca.y()[index(Vars::PPI)], 0.66); + setValue(repca.y()[index(Vars::PEXT)], 0.61); + + setValue(repca.yp()[index(Vars::VMEAS)], 0.01); + setValue(repca.yp()[index(Vars::QMEAS)], -0.02); + setValue(repca.yp()[index(Vars::XQPI)], 0.03); + setValue(repca.yp()[index(Vars::XQLAG)], -0.04); + setValue(repca.yp()[index(Vars::PMEAS)], 0.02); + setValue(repca.yp()[index(Vars::XPPI)], -0.01); + setValue(repca.yp()[index(Vars::PREF)], 0.05); + } + + std::vector dependencyTrackingJacobian() + { + using DepVar = DependencyTracking::Variable; + + auto data = makeResidualData(); + + PhasorDynamics::Bus bus(DepVar{0.9}, DepVar{0.4}); + PhasorDynamics::Converter::Repca repca(&bus, data); + + DepVar ir{0.08}; + DepVar ii{-0.02}; + DepVar p{0.41}; + DepVar q{0.13}; + DepVar freq{0.99}; + DepVar freqref{1.0}; + DepVar vref{1.01}; + DepVar qref{0.12}; + DepVar pplantref{0.55}; + IdxT ir_i = 24; + IdxT ii_i = 25; + IdxT p_i = 26; + IdxT q_i = 27; + IdxT f_i = 28; + IdxT freqref_i = 29; + IdxT vref_i = 30; + IdxT qref_i = 31; + IdxT plantref_i = 32; + + ir.setVariableNumber(ir_i); + ii.setVariableNumber(ii_i); + p.setVariableNumber(p_i); + q.setVariableNumber(q_i); + freq.setVariableNumber(f_i); + freqref.setVariableNumber(freqref_i); + vref.setVariableNumber(vref_i); + qref.setVariableNumber(qref_i); + pplantref.setVariableNumber(plantref_i); + + PhasorDynamics::SignalNode ir_node; + PhasorDynamics::SignalNode ii_node; + PhasorDynamics::SignalNode p_node; + PhasorDynamics::SignalNode q_node; + PhasorDynamics::SignalNode freq_node; + PhasorDynamics::SignalNode freqref_node; + PhasorDynamics::SignalNode vref_node; + PhasorDynamics::SignalNode qref_node; + PhasorDynamics::SignalNode plantref_node; + ir_node.set(&ir, &ir_i); + ii_node.set(&ii, &ii_i); + p_node.set(&p, &p_i); + q_node.set(&q, &q_i); + freq_node.set(&freq, &f_i); + freqref_node.set(&freqref, &freqref_i); + vref_node.set(&vref, &vref_i); + qref_node.set(&qref, &qref_i); + plantref_node.set(&pplantref, &plantref_i); + + repca.getSignals().template attachSignalNode(&ir_node); + repca.getSignals().template attachSignalNode(&ii_node); + repca.getSignals().template attachSignalNode(&p_node); + repca.getSignals().template attachSignalNode(&q_node); + repca.getSignals().template attachSignalNode(&freq_node); + repca.getSignals().template attachSignalNode(&freqref_node); + repca.getSignals().template attachSignalNode(&vref_node); + repca.getSignals().template attachSignalNode(&qref_node); + repca.getSignals().template attachSignalNode(&plantref_node); + + bus.allocate(); + repca.allocate(); + for (IdxT i = 0; i < repca.size(); ++i) + { + repca.y()[static_cast(i)].setVariableNumber(i); + repca.yp()[static_cast(i)].setVariableNumber(i); + } + for (IdxT i = 0; i < bus.size(); ++i) + { + bus.y()[static_cast(i)].setVariableNumber(i + repca.size()); + } + + setJacobianState(repca, bus); + repca.evaluateResidual(); + + std::vector dependencies( + static_cast(repca.size())); + for (IdxT i = 0; i < repca.size(); ++i) + { + dependencies[static_cast(i)] = + repca.getResidual()[static_cast(i)].getDependencies(); + } + return dependencies; + } + + std::vector enzymeJacobian() + { + auto data = makeResidualData(); + + PhasorDynamics::Bus bus(0.9, 0.4); + PhasorDynamics::Converter::Repca repca(&bus, data); + + ScalarT ir{0.08}; + ScalarT ii{-0.02}; + ScalarT p{0.41}; + ScalarT q{0.13}; + ScalarT freq{0.99}; + ScalarT freqref{1.0}; + ScalarT vref{1.01}; + ScalarT qref{0.12}; + ScalarT pplantref{0.55}; + IdxT ir_i = 24; + IdxT ii_i = 25; + IdxT p_i = 26; + IdxT q_i = 27; + IdxT f_i = 28; + IdxT freqref_i = 29; + IdxT vref_i = 30; + IdxT qref_i = 31; + IdxT plantref_i = 32; + + PhasorDynamics::SignalNode ir_node; + PhasorDynamics::SignalNode ii_node; + PhasorDynamics::SignalNode p_node; + PhasorDynamics::SignalNode q_node; + PhasorDynamics::SignalNode freq_node; + PhasorDynamics::SignalNode freqref_node; + PhasorDynamics::SignalNode vref_node; + PhasorDynamics::SignalNode qref_node; + PhasorDynamics::SignalNode plantref_node; + ir_node.set(&ir, &ir_i); + ii_node.set(&ii, &ii_i); + p_node.set(&p, &p_i); + q_node.set(&q, &q_i); + freq_node.set(&freq, &f_i); + freqref_node.set(&freqref, &freqref_i); + vref_node.set(&vref, &vref_i); + qref_node.set(&qref, &qref_i); + plantref_node.set(&pplantref, &plantref_i); + + repca.getSignals().template attachSignalNode(&ir_node); + repca.getSignals().template attachSignalNode(&ii_node); + repca.getSignals().template attachSignalNode(&p_node); + repca.getSignals().template attachSignalNode(&q_node); + repca.getSignals().template attachSignalNode(&freq_node); + repca.getSignals().template attachSignalNode(&freqref_node); + repca.getSignals().template attachSignalNode(&vref_node); + repca.getSignals().template attachSignalNode(&qref_node); + repca.getSignals().template attachSignalNode(&plantref_node); + + bus.allocate(); + repca.allocate(); + for (IdxT i = 0; i < bus.size(); ++i) + { + bus.setVariableIndex(i, i + repca.size()); + bus.setResidualIndex(i, i + repca.size()); + } + + setJacobianState(repca, bus); + repca.updateTime(0.0, 1.0); + repca.evaluateResidual(); + repca.evaluateJacobian(); + repca.constructCsr(); + return MapFromCsr(repca.getCsrJacobian()); + } +#endif + }; + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/runConverterRepcaTests.cpp b/tests/UnitTests/PhasorDynamics/runConverterRepcaTests.cpp new file mode 100644 index 000000000..3103f8106 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/runConverterRepcaTests.cpp @@ -0,0 +1,19 @@ +#include "ConverterRepcaTests.hpp" + +int main() +{ + GridKit::Testing::TestingResults result; + + GridKit::Testing::ConverterRepcaTests test; + + result += test.constructionAndValidation(); + result += test.signalVerification(); + result += test.initializationAndResidual(); + result += test.residualEquations(); + result += test.jsonParseAndSystemAssembly(); +#ifdef GRIDKIT_ENABLE_ENZYME + result += test.jacobian(); +#endif + + return result.summary(); +}