diff --git a/CHANGELOG.md b/CHANGELOG.md index ba0f1dde6..b5ec500ff 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 `REECB` converter model for PhasorDynamics. ## v0.1 diff --git a/GridKit/CommonMath.hpp b/GridKit/CommonMath.hpp index e4a7ec87f..f3e10d6d3 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) @@ -315,7 +339,8 @@ namespace GridKit * @brief Smooth anti-windup indicator for a limited state variable * * @tparam ScalarT - Scalar data type - * @tparam RealT - Real data type (see GridKit::ScalarTraits::RealT) + * @tparam LowerT - data type of the lower limit + * @tparam UpperT - data type of the upper limit * * @param[in] x - State variable * @param[in] f - Pre-limit derivative of the state variable @@ -323,14 +348,19 @@ namespace GridKit * @param[in] limit_max - Maximum limit * @return Scalar value in [0, 1]: 1 when dynamics should pass through, * 0 when integration should be blocked. + * + * @note The limit types intentionally may differ from the scalar type so + * that constant Real limits and algebraic-variable limits both work. */ - template + template __attribute__((always_inline)) inline ScalarT indicator( const ScalarT x, const ScalarT f, - const RealT limit_min, - const RealT limit_max) + const LowerT limit_min, + const UpperT limit_max) { + using RealT = typename GridKit::ScalarTraits::RealT; + assert(limit_min <= limit_max); ScalarT above_min = above(x, limit_min); @@ -350,20 +380,24 @@ namespace GridKit * and blocks motion that would push further into saturation. * * @tparam ScalarT - Scalar data type - * @tparam RealT - Real data type (see GridKit::ScalarTraits::RealT) + * @tparam LowerT - data type of the lower limit + * @tparam UpperT - data type of the upper limit * * @param[in] x - Limited state or limited output signal * @param[in] f - Pre-limit derivative * @param[in] limit_min - Minimum limit * @param[in] limit_max - Maximum limit * @return Smooth anti-windup limited derivative + * + * @note The limit types intentionally may differ from the scalar type so + * that constant Real limits and algebraic-variable limits both work. */ - template + template __attribute__((always_inline)) inline ScalarT antiwindup( const ScalarT x, const ScalarT f, - const RealT limit_min, - const RealT limit_max) + const LowerT limit_min, + const UpperT limit_max) { return indicator(x, f, limit_min, limit_max) * f; } 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..e7df7a31f 100644 --- a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp +++ b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include diff --git a/GridKit/Model/PhasorDynamics/Converter/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Converter/CMakeLists.txt new file mode 100644 index 000000000..2002173ee --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/CMakeLists.txt @@ -0,0 +1,6 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +add_subdirectory(REECB) diff --git a/GridKit/Model/PhasorDynamics/Converter/README.md b/GridKit/Model/PhasorDynamics/Converter/README.md index ad38ba19a..ecd71d986 100644 --- a/GridKit/Model/PhasorDynamics/Converter/README.md +++ b/GridKit/Model/PhasorDynamics/Converter/README.md @@ -9,6 +9,6 @@ models and the bus equations, typically through commanded active and reactive cu The GridKit converter documentation includes: -- Renewable Energy Generator/Converter Model REGCA (See [REGCA](REGCA/README.md)) - Renewable Energy Generator/Converter Model REGCB (See [REGCB](REGCB/README.md)) - Renewable Energy Electrical Control Model REECA (See [REECA](REECA/README.md)) +- Renewable Energy Electrical Control Model REECB (See [REECB](REECB/README.md)) diff --git a/GridKit/Model/PhasorDynamics/Converter/REECB/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Converter/REECB/CMakeLists.txt new file mode 100644 index 000000000..03492a503 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REECB/CMakeLists.txt @@ -0,0 +1,54 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +set(_install_headers Reecb.hpp ReecbData.hpp) + +if(GRIDKIT_ENABLE_ENZYME) + gridkit_add_library( + phasor_dynamics_converter_reecb + SOURCES ReecbEnzyme.cpp + HEADERS ${_install_headers} + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal + PRIVATE + ClangEnzymeFlags + COMPILE_OPTIONS + PRIVATE + -mllvm + -enzyme-auto-sparsity=1 + -fno-math-errno) +else() + gridkit_add_library( + phasor_dynamics_converter_reecb + SOURCES Reecb.cpp + HEADERS ${_install_headers} + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal) +endif() + +gridkit_add_library( + phasor_dynamics_converter_reecb_dependency_tracking + SOURCES ReecbDependencyTracking.cpp + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal_dependency_tracking) + +target_link_libraries( + phasor_dynamics_components + INTERFACE GridKit::phasor_dynamics_converter_reecb) +target_link_libraries( + phasor_dynamics_components_dependency_tracking + INTERFACE GridKit::phasor_dynamics_converter_reecb_dependency_tracking) diff --git a/GridKit/Model/PhasorDynamics/Converter/REECB/README.md b/GridKit/Model/PhasorDynamics/Converter/REECB/README.md new file mode 100644 index 000000000..fb4b86525 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REECB/README.md @@ -0,0 +1,456 @@ +# **Renewable Energy Electrical Control Model (REECB)** + +REECB is a WECC renewable electrical-control model for inverter-coupled +resources. + +## Notes + +- When used with REPCA active-power control, connect REPCA `pext` to REECB `pref`. + +## Block Diagram + +Standard REECB block diagram. + +![](../../../../../docs/Figures/PhasorDynamics/REECB_diagram.png) + +Figure 1: REECB block diagram. Figure courtesy of [PowerWorld](https://www.powerworld.com/WebHelp/) + +## Model Parameters + +Symbol | Units | JSON | Description | Typical Value | Note +------------------------------------|----------|----------|---------------------------------------------------------|---------------|------ +$S^\mathrm{base}$ | [MVA] | `mva` | REECB component power base | 100.0 | Required positive value; block name: `MVABase` +$s_{\mathrm{pf}}$ | [binary] | `PfFlag` | Power-factor control flag | 0 | Block name: `PfFlag`; 1 = power-factor control, 0 = Q control +$s_V$ | [binary] | `VFlag` | Voltage-control mode flag | 0 | Block name: `VFlag`; 1 = Q control, 0 = voltage control +$s_Q$ | [binary] | `QFlag` | Reactive-power control flag | 0 | Block name: `QFlag`; 1 = voltage/Q control, 0 = constant pf or Q control +$s_{PQ}$ | [binary] | `Pqflag` | P/Q priority flag for converter current limit | 0 | Block name: `Pqflag`; 0 = Q priority, 1 = P priority +$T_{\mathrm{rv}}$ | [sec] | `Trv` | Voltage-measurement filter time constant | 0.02 | State 1 +$T_{\mathrm{p}}$ | [sec] | `Tp` | Electrical-power measurement filter time constant | 0.0 | State 2 +$V_0^\mathrm{ref}$ | [p.u.] | `Vref0` | Outer-loop voltage reference | $V_{T,0}$ | Initialized from terminal voltage if omitted +$V_{\mathrm{dip}}$ | [p.u.] | `Vdip` | Low-voltage threshold for the voltage-band gate | 0.85 | +$V_{\mathrm{up}}$ | [p.u.] | `Vup` | High-voltage threshold for the voltage-band gate | 1.15 | +$D_1^\mathrm{db}$ | [p.u.] | `dbd1` | Lower deadband threshold for voltage-error response | 0.0 | +$D_2^\mathrm{db}$ | [p.u.] | `dbd2` | Upper deadband threshold for voltage-error response | 0.0 | +$K_{\mathrm{qv}}$ | [p.u.] | `kqv` | Reactive-current injection gain outside the voltage band | 5.0 | +$I_{q,\mathrm{inj}}^{\min}$ | [p.u.] | `Iql1` | Minimum reactive-current injection limit | -1.1 | +$I_{q,\mathrm{inj}}^{\max}$ | [p.u.] | `Iqh1` | Maximum reactive-current injection limit | 1.1 | +$Q^{\max}$ | [p.u.] | `Qmax` | Maximum reactive-power control limit | 0.436 | +$Q^{\min}$ | [p.u.] | `Qmin` | Minimum reactive-power control limit | -0.436 | +$K_{\mathrm{qp}}$ | [p.u.] | `Kqp` | Reactive-power control proportional gain | 0.0 | +$K_{\mathrm{qi}}$ | [p.u./s] | `Kqi` | Reactive-power control integral gain | 0.1 | +$V^{\max}$ | [p.u.] | `Vmax` | Maximum voltage-control limit | 1.1 | +$V^{\min}$ | [p.u.] | `Vmin` | Minimum voltage-control limit | 0.9 | +$K_{\mathrm{vp}}$ | [p.u.] | `Kvp` | Voltage-control proportional gain | 18.0 | +$K_{\mathrm{vi}}$ | [p.u./s] | `Kvi` | Voltage-control integral gain | 5.0 | +$T_{\mathrm{iq}}$ | [sec] | `Tiq` | Reactive-current command lag time constant | 0.02 | State 5 +$T_{\mathrm{pord}}$ | [sec] | `Tpord` | Active-power order filter time constant | 0.02 | State 6 +$R_P^{\max}$ | [p.u./s] | `dPmax` | Positive active-power order ramp-rate limit | 99.0 | +$R_P^{\min}$ | [p.u./s] | `dPmin` | Negative active-power order ramp-rate limit | -99.0 | +$P^{\max}$ | [p.u.] | `Pmax` | Maximum active-power order limit | 1.0 | +$P^{\min}$ | [p.u.] | `Pmin` | Minimum active-power order limit | 0.0 | +$I^{\max}$ | [p.u.] | `Imax` | Maximum total converter current | 1.3 | + +### Parameter Validation + +Invalid REECB parameter sets are rejected by the following checks. The displayed +equations use effective time constants with $\epsilon_T=10^{-3}$. + +```math +\begin{aligned} + T &\leftarrow \max\!\left(T, \epsilon_T\right) + \quad T\in\{T_{\mathrm{rv}},T_{\mathrm{p}}\} \\ + S^\mathrm{base} &> 0 \\ + s_{\mathrm{pf}}, s_V, s_Q, s_{PQ} + &\in \{0,1\} \\ + T_{\mathrm{rv}}, T_{\mathrm{p}} + &\ge 0 \\ + V_{\mathrm{dip}} + &< V_{\mathrm{up}} \\ + D_1^\mathrm{db} + &\le 0 \le D_2^\mathrm{db} \\ + I_{q,\mathrm{inj}}^{\min} + &\le I_{q,\mathrm{inj}}^{\max} \\ + Q^{\min} + &\le Q^{\max} \\ + V^{\min} + &\le V^{\max} \\ + T_{\mathrm{iq}}, T_{\mathrm{pord}} + &> 0 \\ + R_P^{\min} + &< 0 < R_P^{\max} \\ + P^{\min} + &\le P^{\max} \\ + I^{\max} + &\ge 0 +\end{aligned} +``` + +### Model Derived Parameters + +```math +\begin{aligned} + s_{\mathrm{pf}}^\mathrm{off} + &= 1 - s_{\mathrm{pf}} \\ + s_V^\mathrm{off} + &= 1 - s_V \\ + s_Q^\mathrm{off} + &= 1 - s_Q \\ + 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 terminal voltage | State 1 in Fig. 1 +$P^\mathrm{meas}$ | [p.u.] | Filtered electrical power | State 2 in Fig. 1 +$x_Q^\mathrm{PI}$ | [p.u.] | Reactive-power PI controller state | State 3 in Fig. 1 +$x_V^\mathrm{PI}$ | [p.u.] | Voltage PI controller state | State 4 in Fig. 1 +$Q_V$ | [p.u.] | Reactive-current command lag state | State 5 in Fig. 1 +$P^\mathrm{ord}$ | [p.u.] | Filtered active-power order | State 6 in Fig. 1 + +#### Algebraic + +Symbol | Units | Description | Note +------------------------------------|----------|-------------------------------------|------ +$V_T$ | [p.u.] | Terminal voltage magnitude | +$V_{\mathrm{safe}}^\mathrm{meas}$ | [p.u.] | Safe filtered terminal voltage for divider blocks | Lower bounded by 0.01 +$s_{\mathrm{dip}}$ | [-] | Voltage inside-band control gate | +$e_V^\mathrm{db}$ | [p.u.] | Deadbanded voltage error | +$I_q^\mathrm{inj}$ | [p.u.] | Reactive-current injection candidate | Component base +$Q^\mathrm{ref}$ | [p.u.] | Selected reactive-power reference | +$e_Q$ | [p.u.] | Reactive-power control error | +$V_Q^\mathrm{PI}$ | [p.u.] | Reactive-power control PI output | +$e_V^\mathrm{PI}$ | [p.u.] | Voltage-control PI error | +$f_P^\mathrm{ord}$ | [p.u./s] | Active-power order derivative target | Before ramp-rate limit +$r_P^\mathrm{ord}$ | [p.u./s] | Ramp-rate-limited active-power order derivative target | +$I_q^\mathrm{circ}$ | [p.u.] | Reactive-current limit from converter current circle | +$I_p^\mathrm{circ}$ | [p.u.] | Active-current limit from converter current circle | +$I_q^{\max}$ | [p.u.] | Final reactive-current upper limit | Component base +$I_p^{\max}$ | [p.u.] | Final active-current upper limit | Component base +$I_q^\mathrm{base}$ | [p.u.] | Base reactive-current command | Component base +$I_q^\mathrm{raw}$ | [p.u.] | Raw reactive-current command before final limit | Component base +$I_q^\mathrm{cmd}$ | [p.u.] | Reactive-current command output | System base +$I_p^\mathrm{cmd}$ | [p.u.] | Active-current command output | System base + +### External Variables + +#### Differential +None. + +#### Algebraic + +Symbol | Units | Type | Description | Note +-------------------------------------|--------|---------|-----------------------------------|------ +$V_{\mathrm{r}}$ | [p.u.] | Known | Terminal voltage, real component | Bus input +$V_{\mathrm{i}}$ | [p.u.] | Known | Terminal voltage, imaginary component | Bus input +$P_e$ | [p.u.] | Unknown | Electrical active-power feedback | Signal port `pe`; system base +$Q^\mathrm{gen}$ | [p.u.] | Unknown | Reactive-power feedback | Signal port `qgen`; system base +$Q^\mathrm{ext}$ | [p.u.] | Unknown | External reactive-power command | Optional signal port `qext`; system base +$\phi^\mathrm{ref}$ | [rad] | Unknown | Power-factor angle reference | Optional signal port `pfaref` +$P^\mathrm{ref}$ | [p.u.] | Unknown | External active-power reference | Optional signal port `pref`; system base + +## Model Equations + +### Differential Equations + +```math +\begin{aligned} + 0 &= + -\dot{V}^\mathrm{meas} + + \dfrac{1}{T_{\mathrm{rv}}} + \left(V_T - V^\mathrm{meas}\right) \\ + 0 &= + -\dot{P}^\mathrm{meas} + + \dfrac{1}{T_{\mathrm{p}}} + \left(k_{\mathrm{base}}P_e - P^\mathrm{meas}\right) \\ + 0 &= + -\dot{x}_Q^\mathrm{PI} + + s_{\mathrm{dip}}\, + \text{antiwindup} + \left( + K_{\mathrm{qp}}e_Q + x_Q^\mathrm{PI},\, + K_{\mathrm{qi}}e_Q;\, + V^{\min}, V^{\max} + \right) \\ + 0 &= + -\dot{x}_V^\mathrm{PI} + + s_{\mathrm{dip}}\, + \text{antiwindup} + \left( + K_{\mathrm{vp}}e_V^\mathrm{PI} + x_V^\mathrm{PI},\, + K_{\mathrm{vi}}e_V^\mathrm{PI};\, + -I_q^{\max}, I_q^{\max} + \right) \\ + 0 &= + -\dot{Q}_V + + \dfrac{s_{\mathrm{dip}}}{T_{\mathrm{iq}}} + \left( + \dfrac{Q^\mathrm{ref}}{V_{\mathrm{safe}}^\mathrm{meas}} + - Q_V + \right) \\ + 0 &= + -\dot{P}^\mathrm{ord} + + s_{\mathrm{dip}}\, + \text{antiwindup} + \left(P^\mathrm{ord}, r_P^\mathrm{ord};\, P^{\min}, P^{\max}\right) +\end{aligned} +``` + +CommonMath defines the [Anti-Windup](../../../../CommonMath.md#anti-windup-indicator) +target and smooth approximation. + +### Algebraic Equations + +```math +\begin{aligned} + 0 &= + -V_T^2 + + V_{\mathrm{r}}^2 + + V_{\mathrm{i}}^2 \\ + 0 &= + -V_{\mathrm{safe}}^\mathrm{meas} + + \text{max} + \left(V^\mathrm{meas}, 0.01\right) \\ + 0 &= + -s_{\mathrm{dip}} + + \text{inside} + \left(V_T;\, V_{\mathrm{dip}}, V_{\mathrm{up}}\right) \\ + 0 &= + -e_V^\mathrm{db} + + \text{deadband2} + \left(V_0^\mathrm{ref} - V^\mathrm{meas};\, + D_1^\mathrm{db}, D_2^\mathrm{db}\right) \\ + 0 &= + -I_q^\mathrm{inj} + + \text{clamp} + \left(K_{\mathrm{qv}}e_V^\mathrm{db};\, + I_{q,\mathrm{inj}}^{\min}, I_{q,\mathrm{inj}}^{\max}\right) \\ + 0 &= + -Q^\mathrm{ref} + + s_{\mathrm{pf}}P^\mathrm{meas}\tan\!\left(\phi^\mathrm{ref}\right) + + s_{\mathrm{pf}}^\mathrm{off}k_{\mathrm{base}}Q^\mathrm{ext} \\ + 0 &= + -e_Q + + \text{clamp} + \left(Q^\mathrm{ref};\, Q^{\min}, Q^{\max}\right) + - k_{\mathrm{base}}Q^\mathrm{gen} \\ + 0 &= + -V_Q^\mathrm{PI} + + \text{clamp} + \left(K_{\mathrm{qp}}e_Q + x_Q^\mathrm{PI};\, + V^{\min}, V^{\max}\right) \\ + 0 &= + -e_V^\mathrm{PI} + + s_V V_Q^\mathrm{PI} + + s_V^\mathrm{off}Q^\mathrm{ref} + - V^\mathrm{meas} \\ + 0 &= + -f_P^\mathrm{ord} + + \dfrac{1}{T_{\mathrm{pord}}} + \left(k_{\mathrm{base}}P^\mathrm{ref} - P^\mathrm{ord}\right) \\ + 0 &= + -r_P^\mathrm{ord} + + \text{clamp} + \left(f_P^\mathrm{ord};\, R_P^{\min}, R_P^{\max}\right) \\ + 0 &= + -\left(I_q^\mathrm{circ}\right)^2 + + \left(I^{\max}\right)^2 + - s_{PQ}\left(k_{\mathrm{base}}I_p^\mathrm{cmd}\right)^2 \\ + 0 &= + -\left(I_p^\mathrm{circ}\right)^2 + + \left(I^{\max}\right)^2 + - \left(1-s_{PQ}\right)\left(k_{\mathrm{base}}I_q^\mathrm{cmd}\right)^2 \\ + 0 &= + -I_q^{\max} + + \left(1-s_{PQ}\right)I^{\max} + + s_{PQ}I_q^\mathrm{circ} \\ + 0 &= + -I_p^{\max} + + s_{PQ}I^{\max} + + \left(1-s_{PQ}\right)I_p^\mathrm{circ} \\ + 0 &= + -I_q^\mathrm{base} + + \text{clamp} + \left(K_{\mathrm{vp}}e_V^\mathrm{PI} + x_V^\mathrm{PI};\, + -I_q^{\max}, I_q^{\max}\right) \\ + 0 &= + -I_q^\mathrm{raw} + + s_Q I_q^\mathrm{base} + + s_Q^\mathrm{off}Q_V + + \left(1-s_{\mathrm{dip}}\right)I_q^\mathrm{inj} \\ + 0 &= + -k_{\mathrm{base}}I_q^\mathrm{cmd} + + \text{clamp} + \left(I_q^\mathrm{raw};\, + -I_q^{\max}, I_q^{\max}\right) \\ + 0 &= + -k_{\mathrm{base}}I_p^\mathrm{cmd} + + \text{clamp} + \left( + \dfrac{P^\mathrm{ord}}{V_{\mathrm{safe}}^\mathrm{meas}};\, + 0,\, + I_p^{\max} + \right) +\end{aligned} +``` + +CommonMath defines helper targets and smooth approximations for +[max, clamp, deadband2, and inside](../../../../CommonMath.md#derived-functions). + +## Initialization + +### External Priors + +```math +\begin{aligned} + V_{\mathrm{r},0}, V_{\mathrm{i},0} + &\leftarrow \text{reference voltage measurement} \\ + I_{q,0}^\mathrm{cmd}, I_{p,0}^\mathrm{cmd} + &\leftarrow \text{current-command start} +\end{aligned} +``` + +### Internal Initialization + +Define + +```math +\begin{aligned} + \text{awinit}(x^\star,f;\ell,u) + &= + \begin{cases} + x^\star & f = 0 \\ + u + \epsilon_{\mathrm{sat}} & f > 0 \\ + \ell - \epsilon_{\mathrm{sat}} & f < 0 + \end{cases} +\end{aligned} +``` + +with $\epsilon_{\mathrm{sat}}>0$. + +Initialization is performed by evaluating the steady-state residuals in +dependency order. Let subscript $0$ denote initial values and set all internal +derivatives to zero: + +```math +\begin{aligned} + V_{T,0} + &= \sqrt{V_{\mathrm{r},0}^2 + V_{\mathrm{i},0}^2} \\ + V_0^\mathrm{ref} + &= V_{T,0}\quad\text{if omitted} \\ + V_0^\mathrm{meas} + &= V_{T,0} \\ + V_{\mathrm{safe},0}^\mathrm{meas} + &= \text{max}\left(V_0^\mathrm{meas}, 0.01\right) \\ + P_0^\mathrm{meas} + &= k_{\mathrm{base}}P_{e,0} \\ + s_{\mathrm{dip},0} + &= \text{inside} + \left(V_{T,0};\, V_{\mathrm{dip}}, V_{\mathrm{up}}\right) \\ + e_{V,0}^\mathrm{db} + &= + \text{deadband2} + \left(V_0^\mathrm{ref} - V_0^\mathrm{meas};\, + D_1^\mathrm{db}, D_2^\mathrm{db}\right) \\ + I_{q,0}^\mathrm{inj} + &= + \text{clamp} + \left(K_{\mathrm{qv}}e_{V,0}^\mathrm{db};\, + I_{q,\mathrm{inj}}^{\min}, I_{q,\mathrm{inj}}^{\max}\right) \\ + Q_0^\mathrm{ref} + &= + s_{\mathrm{pf}}P_0^\mathrm{meas} + \tan\!\left(\phi_0^\mathrm{ref}\right) + + s_{\mathrm{pf}}^\mathrm{off}k_{\mathrm{base}}Q_0^\mathrm{ext} \\ + e_{Q,0} + &= + \text{clamp} + \left(Q_0^\mathrm{ref};\, Q^{\min}, Q^{\max}\right) + - k_{\mathrm{base}}Q_0^\mathrm{gen} \\ + Q_{V,0} + &= \dfrac{Q_0^\mathrm{ref}}{V_{\mathrm{safe},0}^\mathrm{meas}} \\ + P_0^\mathrm{ord} + &= k_{\mathrm{base}}P_0^\mathrm{ref} \\ + f_{P,0}^\mathrm{ord} + &= 0 \\ + r_{P,0}^\mathrm{ord} + &= 0 \\ + u_{Q,0}^\mathrm{PI} + &= + \text{awinit} + \left( + s_V V_0^\mathrm{meas} + + s_V^\mathrm{off}Q_0^\mathrm{ref},\, + K_{\mathrm{qi}}e_{Q,0};\, + V^{\min}, V^{\max} + \right) \\ + V_{Q,0}^\mathrm{PI} + &= + \text{clamp} + \left(u_{Q,0}^\mathrm{PI};\, V^{\min}, V^{\max}\right) \\ + e_{V,0}^\mathrm{PI} + &= + s_V V_{Q,0}^\mathrm{PI} + + s_V^\mathrm{off}Q_0^\mathrm{ref} + - V_0^\mathrm{meas} \\ + x_{Q,0}^\mathrm{PI} + &= u_{Q,0}^\mathrm{PI} - K_{\mathrm{qp}}e_{Q,0} \\ + u_{V,0}^\mathrm{PI} + &= + \text{awinit} + \left( + k_{\mathrm{base}}Q_0^\mathrm{gen}/V_{\mathrm{safe},0}^\mathrm{meas},\, + K_{\mathrm{vi}}e_{V,0}^\mathrm{PI};\, + -I_{q,0}^{\max}, I_{q,0}^{\max} + \right) \\ + I_{q,0}^\mathrm{base} + &= + \text{clamp} + \left(u_{V,0}^\mathrm{PI};\, + -I_{q,0}^{\max}, + I_{q,0}^{\max}\right) \\ + x_{V,0}^\mathrm{PI} + &= u_{V,0}^\mathrm{PI} - K_{\mathrm{vp}}e_{V,0}^\mathrm{PI} +\end{aligned} +``` + +Initialization rejects negative current-circle radicands. + +### External Solved + +```math +\begin{aligned} + P_{e,0} + &\leftarrow V_{\mathrm{safe},0}^\mathrm{meas} I_{p,0}^\mathrm{cmd} \\ + Q_0^\mathrm{gen} + &\leftarrow V_{\mathrm{safe},0}^\mathrm{meas} I_{q,0}^\mathrm{cmd} \\ + Q_0^\mathrm{ext} + &\leftarrow Q_0^\mathrm{gen} \\ + \phi_0^\mathrm{ref} + &\leftarrow + \begin{cases} + \tan^{-1}\!\left(Q_0^\mathrm{gen}/P_{e,0}\right) & P_{e,0} \ne 0 \\ + 0 & P_{e,0} = 0 + \end{cases} \\ + P_0^\mathrm{ref} + &\leftarrow + \dfrac{1}{k_{\mathrm{base}}}\text{clamp} + \left(k_{\mathrm{base}}P_{e,0};\, P^{\min}, P^{\max}\right) +\end{aligned} +``` + +## Monitorable Outputs + +Output | Units | Description | Note +----------------|--------|-------------------------------------|------ +`iqcmd` | [p.u.] | Reactive-current command output | $I_q^\mathrm{cmd}$ (system base) +`ipcmd` | [p.u.] | Active-current command output | $I_p^\mathrm{cmd}$ (system base) +`vmeas` | [p.u.] | Filtered terminal voltage | $V^\mathrm{meas}$ +`pmeas` | [p.u.] | Filtered electrical power | $P^\mathrm{meas}$ (component base) diff --git a/GridKit/Model/PhasorDynamics/Converter/REECB/Reecb.cpp b/GridKit/Model/PhasorDynamics/Converter/REECB/Reecb.cpp new file mode 100644 index 000000000..3b6cacabf --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REECB/Reecb.cpp @@ -0,0 +1,27 @@ +/** + * @file Reecb.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Non-Enzyme instantiation for the REECB electrical-control model. + */ + +#include "ReecbImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + template + int Reecb::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Reecb..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Reecb; + template class Reecb; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REECB/Reecb.hpp b/GridKit/Model/PhasorDynamics/Converter/REECB/Reecb.hpp new file mode 100644 index 000000000..7172ef229 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REECB/Reecb.hpp @@ -0,0 +1,200 @@ +/** + * @file Reecb.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Declaration of the REECB electrical-control model. + */ + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + template + class BusBase; + + template + class SignalNode; + + namespace Converter + { + /// Internal variables of a `Reecb`. + enum class ReecbInternalVariables : size_t + { + VMEAS, ///< Filtered terminal voltage + PMEAS, ///< Filtered electrical power + XPIQ, ///< Reactive-power PI state + XPIV, ///< Voltage PI state + QV, ///< Reactive-current command lag state + PORD, ///< Filtered active-power order + VT, ///< Terminal voltage magnitude + VMEASSAFE, ///< Safe filtered terminal voltage + SDIP, ///< Voltage inside-band control gate + VERR, ///< Deadbanded voltage error + IQV, ///< Reactive-current injection candidate + QREF, ///< Selected reactive-power reference + EQ, ///< Reactive-power control error + VPIQ, ///< Reactive-power PI output + EPIV, ///< Voltage-control PI error + FPORD, ///< Active-power order derivative before ramp limiting + RPORD, ///< Ramp-rate-limited active-power derivative + IQCIRC, ///< Reactive-current limit from current circle + IPCIRC, ///< Active-current limit from current circle + IQMAX, ///< Final reactive-current upper limit + IPMAX, ///< Final active-current upper limit + IQBASE, ///< Base reactive-current command + IQRAW, ///< Raw reactive-current command + IQCMD, ///< Reactive-current command output + IPCMD, ///< Active-current command output + MAXIMUM, + }; + + /// External variables of a `Reecb`. + enum class ReecbExternalVariables : size_t + { + PE, ///< Electrical active-power signal + QGEN, ///< Reactive-power signal + QEXT, ///< External reactive-power command + PFAREF, ///< Power-factor angle reference + PREF, ///< External active-power reference + MAXIMUM, + }; + + template + class Reecb : public Component + { + using Component::alpha_; + using Component::abs_tol_; + 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 = ReecbData; + using MonitorT = Model::VariableMonitor; + + Reecb(BusT* bus); + Reecb(BusT* bus, const ModelDataT& data); + ~Reecb(); + + int setGridKitComponentID(IdxT) override final; + int allocate() override final; + int verify() const override final; + int initialize() override final; + int tagDifferentiable() override final; + int 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 initModelParams(const ModelDataT& data); + void initializeMonitor(); + void setDerivedParameters(); + + ScalarT toComponentBase(ScalarT value) const; + ScalarT toSystemBase(ScalarT value) const; + + ScalarT& Vr(); + ScalarT& Vi(); + + static constexpr RealT TIME_CONSTANT_MINIMUM = static_cast(1.0e-3); + static constexpr RealT VMEAS_MINIMUM = static_cast(0.01); + static constexpr RealT INIT_TOL = static_cast(1.0e-10); + static constexpr RealT SAT_MARGIN = static_cast(0.1); + + BusT* bus_{nullptr}; + + RealT mva_base_{static_cast(100.0)}; + RealT PfFlag_{ZERO}; + RealT VFlag_{ZERO}; + RealT QFlag_{ZERO}; + RealT Pqflag_{ZERO}; + RealT Trv_{static_cast(0.02)}; + RealT Tp_{ZERO}; + RealT Vref0_{ZERO}; + RealT Vdip_{static_cast(0.85)}; + RealT Vup_{static_cast(1.15)}; + RealT dbd1_{ZERO}; + RealT dbd2_{ZERO}; + RealT kqv_{static_cast(5.0)}; + RealT Iql1_{static_cast(-1.1)}; + RealT Iqh1_{static_cast(1.1)}; + RealT Qmax_{static_cast(0.436)}; + RealT Qmin_{static_cast(-0.436)}; + RealT Kqp_{ZERO}; + RealT Kqi_{static_cast(0.1)}; + RealT Vmax_{static_cast(1.1)}; + RealT Vmin_{static_cast(0.9)}; + RealT Kvp_{static_cast(18.0)}; + RealT Kvi_{static_cast(5.0)}; + RealT Tiq_{static_cast(0.02)}; + RealT Tpord_{static_cast(0.02)}; + RealT dPmax_{static_cast(99.0)}; + RealT dPmin_{static_cast(-99.0)}; + RealT Pmax_{ONE}; + RealT Pmin_{ZERO}; + RealT Imax_{static_cast(1.3)}; + RealT va_converter_base_{0}; + RealT Trv_eff_{TIME_CONSTANT_MINIMUM}; + RealT Tp_eff_{TIME_CONSTANT_MINIMUM}; + RealT pf_off_{1}; + RealT v_off_{1}; + RealT q_off_{1}; + + bool Vref0_given_{false}; + IdxT parameter_error_count_{0}; + + ScalarT qext_set_{0}; + ScalarT pe_set_{0}; + ScalarT qgen_set_{0}; + ScalarT pfaref_set_{0}; + ScalarT pref_set_{0}; + + ComponentSignals signals_; + std::unique_ptr monitor_; + + std::vector ws_; + std::vector ws_indices_; + }; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbData.hpp b/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbData.hpp new file mode 100644 index 000000000..11f7cab41 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbData.hpp @@ -0,0 +1,89 @@ +/** + * @file ReecbData.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Modeling data for the REECB electrical-control model. + */ + +#pragma once + +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + /// Parameter keys for the REECB electrical-control model. + enum class ReecbParameters + { + mva, ///< REECB component power base + PfFlag, ///< Power-factor control flag + VFlag, ///< Voltage-control mode flag + QFlag, ///< Reactive-power control flag + Pqflag, ///< P/Q current-priority flag + Trv, ///< Voltage-measurement filter time constant + Tp, ///< Electrical-power measurement filter time constant + Vref0, ///< Outer-loop voltage reference + Vdip, ///< Low-voltage dip threshold + Vup, ///< High-voltage threshold + dbd1, ///< Lower voltage-error deadband threshold + dbd2, ///< Upper voltage-error deadband threshold + kqv, ///< Reactive-current injection gain + Iql1, ///< Minimum reactive-current injection limit + Iqh1, ///< Maximum reactive-current injection limit + Qmax, ///< Maximum reactive-power control limit + Qmin, ///< Minimum reactive-power control limit + Kqp, ///< Reactive-power proportional gain + Kqi, ///< Reactive-power integral gain + Vmax, ///< Maximum voltage-control limit + Vmin, ///< Minimum voltage-control limit + Kvp, ///< Voltage-control proportional gain + Kvi, ///< Voltage-control integral gain + Tiq, ///< Reactive-current command lag time constant + Tpord, ///< Active-power order filter time constant + dPmax, ///< Positive active-power ramp-rate limit + dPmin, ///< Negative active-power ramp-rate limit + Pmax, ///< Maximum active-power order limit + Pmin, ///< Minimum active-power order limit + Imax ///< Maximum converter current + }; + + /// Ports for the REECB electrical-control model. + enum class ReecbPorts + { + bus, ///< Terminal bus ID + pe, ///< Electrical active-power signal ID + qgen, ///< Reactive-power signal ID + qext, ///< Optional reactive-power command signal ID + pfaref, ///< Optional power-factor angle reference signal ID + pref, ///< Optional active-power reference signal ID + iqcmd, ///< Reactive-current command output signal ID + ipcmd ///< Active-current command output signal ID + }; + + /// Variables available through the monitor interface. + enum class ReecbMonitorableVariables + { + iqcmd, ///< Reactive-current command output + ipcmd, ///< Active-current command output + vmeas, ///< Filtered terminal voltage + pmeas ///< Filtered electrical power + }; + + template + struct ReecbData : public ComponentData + { + ReecbData() = default; + + using Parameters = ReecbParameters; + using Ports = ReecbPorts; + using MonitorableVariables = ReecbMonitorableVariables; + }; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbDependencyTracking.cpp new file mode 100644 index 000000000..2f81e8274 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbDependencyTracking.cpp @@ -0,0 +1,27 @@ +/** + * @file ReecbDependencyTracking.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Dependency-tracking instantiations for the REECB electrical-control model. + */ + +#include "ReecbImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + template + int Reecb::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Reecb..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Reecb; + template class Reecb; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbEnzyme.cpp b/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbEnzyme.cpp new file mode 100644 index 000000000..471bd0014 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbEnzyme.cpp @@ -0,0 +1,104 @@ +/** + * @file ReecbEnzyme.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Enzyme sparse Jacobian for the REECB electrical-control model. + */ + +#include + +#include "ReecbImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + template + int Reecb::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Reecb..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + + if (J_rows_buffer_ == nullptr) + { + auto size = static_cast(size_); + auto bus_size = static_cast(bus_->size()); + auto signal_size = 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 ModelT = GridKit::PhasorDynamics::Converter::Reecb; + 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 Reecb; + template class Reecb; + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbImpl.hpp b/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbImpl.hpp new file mode 100644 index 000000000..92b25f828 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Converter/REECB/ReecbImpl.hpp @@ -0,0 +1,662 @@ +/** + * @file ReecbImpl.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Definition of the REECB electrical-control model. + */ + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Converter + { + using Log = ::GridKit::Utilities::Logger; + + template + Reecb::Reecb(BusT* bus) + : bus_(bus) + { + size_ = static_cast(ReecbInternalVariables::MAXIMUM); + setDerivedParameters(); + } + + template + Reecb::Reecb(BusT* bus, const ModelDataT& data) + : bus_(bus), + monitor_(std::make_unique(data)) + { + initModelParams(data); + initializeMonitor(); + size_ = static_cast(ReecbInternalVariables::MAXIMUM); + } + + template + Reecb::~Reecb() + { + } + + template + scalar_type& Reecb::Vr() + { + return bus_->Vr(); + } + + template + scalar_type& Reecb::Vi() + { + return bus_->Vi(); + } + + template + void Reecb::setDerivedParameters() + { + va_converter_base_ = mva_base_ * static_cast(1.0e6); + Trv_eff_ = std::max(Trv_, TIME_CONSTANT_MINIMUM); + Tp_eff_ = std::max(Tp_, TIME_CONSTANT_MINIMUM); + pf_off_ = ONE - PfFlag_; + v_off_ = ONE - VFlag_; + q_off_ = ONE - QFlag_; + } + + template + scalar_type Reecb::toComponentBase( + scalar_type value) const + { + return value * va_system_base_ / va_converter_base_; + } + + template + scalar_type Reecb::toSystemBase( + scalar_type value) const + { + return value / toComponentBase(static_cast(ONE)); + } + + template + void Reecb::initModelParams(const ModelDataT& data) + { + using Params = typename ModelDataT::Parameters; + + parameter_error_count_ = 0; + Vref0_given_ = false; + + 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() << "Reecb: parameter '" << name << "' must be numeric\n"; + ++parameter_error_count_; + } + }; + + auto load_switch = [&](auto key, RealT& 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 = ZERO; + if (*bool_value) + { + target = ONE; + } + } + else if (const auto* index_value = std::get_if(&value); + index_value && (*index_value == 0 || *index_value == 1)) + { + target = static_cast(*index_value); + } + else if (const auto* real_value = std::get_if(&value); + real_value && (*real_value == ZERO || *real_value == ONE) ) + { + target = *real_value; + } + else + { + Log::error() << "Reecb: parameter '" << name << "' must be bool or 0/1\n"; + ++parameter_error_count_; + } + }; + + if (!data.parameters.contains(Params::mva)) + { + Log::error() << "Reecb: missing required parameter 'mva'\n"; + ++parameter_error_count_; + } + load_real(Params::mva, mva_base_, "mva"); + load_switch(Params::PfFlag, PfFlag_, "PfFlag"); + load_switch(Params::VFlag, VFlag_, "VFlag"); + load_switch(Params::QFlag, QFlag_, "QFlag"); + load_switch(Params::Pqflag, Pqflag_, "Pqflag"); + load_real(Params::Trv, Trv_, "Trv"); + load_real(Params::Tp, Tp_, "Tp"); + if (data.parameters.contains(Params::Vref0)) + { + load_real(Params::Vref0, Vref0_, "Vref0"); + Vref0_given_ = true; + } + load_real(Params::Vdip, Vdip_, "Vdip"); + load_real(Params::Vup, Vup_, "Vup"); + load_real(Params::dbd1, dbd1_, "dbd1"); + load_real(Params::dbd2, dbd2_, "dbd2"); + load_real(Params::kqv, kqv_, "kqv"); + load_real(Params::Iql1, Iql1_, "Iql1"); + load_real(Params::Iqh1, Iqh1_, "Iqh1"); + load_real(Params::Qmax, Qmax_, "Qmax"); + load_real(Params::Qmin, Qmin_, "Qmin"); + load_real(Params::Kqp, Kqp_, "Kqp"); + load_real(Params::Kqi, Kqi_, "Kqi"); + load_real(Params::Vmax, Vmax_, "Vmax"); + load_real(Params::Vmin, Vmin_, "Vmin"); + load_real(Params::Kvp, Kvp_, "Kvp"); + load_real(Params::Kvi, Kvi_, "Kvi"); + load_real(Params::Tiq, Tiq_, "Tiq"); + load_real(Params::Tpord, Tpord_, "Tpord"); + load_real(Params::dPmax, dPmax_, "dPmax"); + load_real(Params::dPmin, dPmin_, "dPmin"); + load_real(Params::Pmax, Pmax_, "Pmax"); + load_real(Params::Pmin, Pmin_, "Pmin"); + load_real(Params::Imax, Imax_, "Imax"); + setDerivedParameters(); + } + + template + const Model::VariableMonitorBase* Reecb::getMonitor() const + { + return monitor_.get(); + } + + template + void Reecb::initializeMonitor() + { + using Variable = typename ModelDataT::MonitorableVariables; + auto index = [](ReecbInternalVariables variable) + { + return static_cast(variable); + }; + + monitor_->set(Variable::iqcmd, [this, index] + { return y_[index(ReecbInternalVariables::IQCMD)]; }); + monitor_->set(Variable::ipcmd, [this, index] + { return y_[index(ReecbInternalVariables::IPCMD)]; }); + monitor_->set(Variable::vmeas, [this, index] + { return y_[index(ReecbInternalVariables::VMEAS)]; }); + monitor_->set(Variable::pmeas, [this, index] + { return y_[index(ReecbInternalVariables::PMEAS)]; }); + } + + template + int Reecb::setGridKitComponentID(IdxT component_id) + { + gridkit_component_id_ = component_id; + return 0; + } + + template + int Reecb::allocate() + { + size_ = static_cast(ReecbInternalVariables::MAXIMUM); + auto size = static_cast(size_); + + f_.assign(size, ScalarT{0}); + y_.assign(size, ScalarT{0}); + yp_.assign(size, ScalarT{0}); + tag_.assign(size, false); + abs_tol_.resize(size); + variable_indices_.resize(size); + residual_indices_.resize(size); + + wb_.assign(2, ScalarT{0}); + + auto signal_size = static_cast(ReecbExternalVariables::MAXIMUM); + ws_.assign(signal_size, ScalarT{0}); + ws_indices_.assign(signal_size, INVALID_INDEX); + + for (IdxT j = 0; j < size_; ++j) + { + this->setVariableIndex(j, j); + this->setResidualIndex(j, j); + } + + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y_[static_cast(ReecbInternalVariables::IQCMD)], + &(this->getVariableIndex(static_cast(ReecbInternalVariables::IQCMD)))); + } + + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y_[static_cast(ReecbInternalVariables::IPCMD)], + &(this->getVariableIndex(static_cast(ReecbInternalVariables::IPCMD)))); + } + + return 0; + } + + template + int Reecb::verify() const + { + int ret = static_cast(parameter_error_count_); + + auto check = [&](bool condition, const char* message) + { + if (!condition) + { + Log::error() << "Reecb: " << message << '\n'; + ret += 1; + } + }; + + if (bus_ == nullptr) + { + Log::error() << "Reecb: bus pointer is null\n"; + ret += 1; + } + + check(mva_base_ > ZERO, "mva must be positive"); + check(va_converter_base_ > ZERO, "converter VA base must be positive"); + check(PfFlag_ == ZERO || PfFlag_ == ONE, "PfFlag must be 0 or 1"); + check(VFlag_ == ZERO || VFlag_ == ONE, "VFlag must be 0 or 1"); + check(QFlag_ == ZERO || QFlag_ == ONE, "QFlag must be 0 or 1"); + check(Pqflag_ == ZERO || Pqflag_ == ONE, "Pqflag must be 0 or 1"); + check(Trv_ >= ZERO, "Trv must be non-negative"); + check(Tp_ >= ZERO, "Tp must be non-negative"); + check(Vdip_ < Vup_, "Vdip must be less than Vup"); + check(dbd1_ <= ZERO && ZERO <= dbd2_, "dbd1 <= 0 <= dbd2 is required"); + check(Iql1_ <= Iqh1_, "Iql1 must be less than or equal to Iqh1"); + check(Qmin_ <= Qmax_, "Qmin must be less than or equal to Qmax"); + check(Vmin_ <= Vmax_, "Vmin must be less than or equal to Vmax"); + check(Tiq_ > ZERO, "Tiq must be positive"); + check(Tpord_ > ZERO, "Tpord must be positive"); + check(dPmin_ < ZERO && ZERO < dPmax_, "dPmin < 0 < dPmax is required"); + check(Pmin_ <= Pmax_, "Pmin must be less than or equal to Pmax"); + check(Imax_ >= ZERO, "Imax must be non-negative"); + + auto check_attached_signal = [&](bool attached, bool linked, const char* name) + { + if (attached && !linked) + { + Log::error() << "Reecb: " << name << " signal attached with no linked variable\n"; + ret += 1; + } + }; + + check_attached_signal( + signals_.template isAttached(), + signals_.template isAttached() + && signals_.template isLinked(), + "pe"); + check_attached_signal( + signals_.template isAttached(), + signals_.template isAttached() + && signals_.template isLinked(), + "qgen"); + check_attached_signal( + signals_.template isAttached(), + signals_.template isAttached() + && signals_.template isLinked(), + "qext"); + check_attached_signal( + signals_.template isAttached(), + signals_.template isAttached() + && signals_.template isLinked(), + "pfaref"); + check_attached_signal( + signals_.template isAttached(), + signals_.template isAttached() + && signals_.template isLinked(), + "pref"); + + return ret; + } + + template + int Reecb::initialize() + { + if (parameter_error_count_ > 0 || verify() > 0) + { + Log::error() << "Reecb: cannot initialize with invalid configuration\n"; + return 1; + } + + const auto VMEAS = static_cast(ReecbInternalVariables::VMEAS); + const auto PMEAS = static_cast(ReecbInternalVariables::PMEAS); + const auto XPIQ = static_cast(ReecbInternalVariables::XPIQ); + const auto XPIV = static_cast(ReecbInternalVariables::XPIV); + const auto QV = static_cast(ReecbInternalVariables::QV); + const auto PORD = static_cast(ReecbInternalVariables::PORD); + const auto VT = static_cast(ReecbInternalVariables::VT); + const auto VMEASSAFE = static_cast(ReecbInternalVariables::VMEASSAFE); + const auto SDIP = static_cast(ReecbInternalVariables::SDIP); + const auto VERR = static_cast(ReecbInternalVariables::VERR); + const auto IQV = static_cast(ReecbInternalVariables::IQV); + const auto QREF = static_cast(ReecbInternalVariables::QREF); + const auto EQ = static_cast(ReecbInternalVariables::EQ); + const auto VPIQ = static_cast(ReecbInternalVariables::VPIQ); + const auto EPIV = static_cast(ReecbInternalVariables::EPIV); + const auto FPORD = static_cast(ReecbInternalVariables::FPORD); + const auto RPORD = static_cast(ReecbInternalVariables::RPORD); + const auto IQCIRC = static_cast(ReecbInternalVariables::IQCIRC); + const auto IPCIRC = static_cast(ReecbInternalVariables::IPCIRC); + const auto IQMAX = static_cast(ReecbInternalVariables::IQMAX); + const auto IPMAX = static_cast(ReecbInternalVariables::IPMAX); + const auto IQBASE = static_cast(ReecbInternalVariables::IQBASE); + const auto IQRAW = static_cast(ReecbInternalVariables::IQRAW); + const auto IQCMD = static_cast(ReecbInternalVariables::IQCMD); + const auto IPCMD = static_cast(ReecbInternalVariables::IPCMD); + + const ScalarT vr = Vr(); + const ScalarT vi = Vi(); + + y_[VT] = std::sqrt(vr * vr + vi * vi); + if (!Vref0_given_) + { + Vref0_ = static_cast(y_[VT]); + } + + y_[VMEAS] = y_[VT]; + y_[VMEASSAFE] = Math::max(y_[VMEAS], VMEAS_MINIMUM); + + const ScalarT ipcmd0 = toComponentBase(y_[IPCMD]); + const ScalarT iqcmd0 = toComponentBase(y_[IQCMD]); + ScalarT pe0 = ipcmd0 * y_[VMEASSAFE]; + ScalarT qgen0 = iqcmd0 * y_[VMEASSAFE]; + + const ScalarT qext0 = qgen0; + const ScalarT pref0 = Math::clamp(pe0, Pmin_, Pmax_); + + pe_set_ = toSystemBase(pe0); + qgen_set_ = toSystemBase(qgen0); + qext_set_ = toSystemBase(qext0); + pfaref_set_ = std::abs(static_cast(pe0)) > INIT_TOL ? static_cast(std::atan(static_cast(qgen0 / pe0))) : static_cast(ZERO); + pref_set_ = toSystemBase(pref0); + + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(pe_set_); + } + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(qgen_set_); + } + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(qext_set_); + } + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(pfaref_set_); + } + if (signals_.template isAttached()) + { + signals_.template writeExternalVariable(pref_set_); + } + + y_[PMEAS] = pe0; + y_[SDIP] = Math::inside(y_[VT], Vdip_, Vup_); + y_[VERR] = Math::deadband2(Vref0_ - y_[VMEAS], dbd1_, dbd2_); + y_[IQV] = Math::clamp(kqv_ * y_[VERR], Iql1_, Iqh1_); + y_[QREF] = PfFlag_ * y_[PMEAS] * std::tan(pfaref_set_) + pf_off_ * qext0; + y_[EQ] = Math::clamp(y_[QREF], Qmin_, Qmax_) - qgen0; + y_[QV] = y_[QREF] / y_[VMEASSAFE]; + y_[PORD] = pref0; + y_[FPORD] = ZERO; + y_[RPORD] = ZERO; + + auto awinit = [](const ScalarT target, const ScalarT rate, const ScalarT lower, const ScalarT upper) -> ScalarT + { + if (std::abs(static_cast(rate)) <= INIT_TOL) + { + return target; + } + return rate > ZERO ? upper + static_cast(SAT_MARGIN) : lower - static_cast(SAT_MARGIN); + }; + + const ScalarT vpiq_arg = awinit(VFlag_ * y_[VMEAS] + v_off_ * y_[QREF], Kqi_ * y_[EQ], static_cast(Vmin_), static_cast(Vmax_)); + y_[VPIQ] = Math::clamp(vpiq_arg, Vmin_, Vmax_); + y_[EPIV] = VFlag_ * y_[VPIQ] + v_off_ * y_[QREF] - y_[VMEAS]; + y_[XPIQ] = vpiq_arg - Kqp_ * y_[EQ]; + + const ScalarT iqbase_target = qgen0 / y_[VMEASSAFE]; + const ScalarT ip_star = y_[PORD] / y_[VMEASSAFE]; + + auto initializeReactiveBase = [&]() + { + const ScalarT piv_arg = awinit(iqbase_target, Kvi_ * y_[EPIV], -y_[IQMAX], y_[IQMAX]); + y_[IQBASE] = Math::clamp(piv_arg, -y_[IQMAX], y_[IQMAX]); + y_[XPIV] = piv_arg - Kvp_ * y_[EPIV]; + }; + + auto initializeReactiveCommand = [&]() + { + y_[IQRAW] = QFlag_ * y_[IQBASE] + q_off_ * y_[QV] + (ONE - y_[SDIP]) * y_[IQV]; + y_[IQCMD] = toSystemBase(Math::clamp(y_[IQRAW], -y_[IQMAX], y_[IQMAX])); + }; + + if (Pqflag_ == ZERO) + { + y_[IQCIRC] = Imax_; + y_[IQMAX] = Imax_; + initializeReactiveBase(); + initializeReactiveCommand(); + + const ScalarT iqcmd = toComponentBase(y_[IQCMD]); + const ScalarT ip_radicand = Imax_ * Imax_ - iqcmd * iqcmd; + if (static_cast(ip_radicand) < ZERO) + { + Log::error() << "Reecb: initial active-current circle radicand is negative\n"; + return 1; + } + y_[IPCIRC] = std::sqrt(ip_radicand); + y_[IPMAX] = y_[IPCIRC]; + y_[IPCMD] = toSystemBase(Math::clamp(ip_star, ZERO, y_[IPMAX])); + } + else + { + y_[IPCIRC] = Imax_; + y_[IPMAX] = Imax_; + y_[IPCMD] = toSystemBase(Math::clamp(ip_star, ZERO, y_[IPMAX])); + + const ScalarT ipcmd = toComponentBase(y_[IPCMD]); + const ScalarT iq_radicand = Imax_ * Imax_ - ipcmd * ipcmd; + if (static_cast(iq_radicand) < ZERO) + { + Log::error() << "Reecb: initial reactive-current circle radicand is negative\n"; + return 1; + } + y_[IQCIRC] = std::sqrt(iq_radicand); + y_[IQMAX] = y_[IQCIRC]; + initializeReactiveBase(); + initializeReactiveCommand(); + } + + std::fill(yp_.begin(), yp_.end(), ZERO); + + return 0; + } + + template + int Reecb::tagDifferentiable() + { + std::fill(tag_.begin(), tag_.end(), false); + tag_[static_cast(ReecbInternalVariables::VMEAS)] = true; + tag_[static_cast(ReecbInternalVariables::PMEAS)] = true; + tag_[static_cast(ReecbInternalVariables::XPIQ)] = true; + tag_[static_cast(ReecbInternalVariables::XPIV)] = true; + tag_[static_cast(ReecbInternalVariables::QV)] = true; + tag_[static_cast(ReecbInternalVariables::PORD)] = true; + return 0; + } + + template + int Reecb::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + + template + __attribute__((always_inline)) inline int + Reecb::evaluateInternalResidual( + const ScalarT* y, + const ScalarT* yp, + const ScalarT* wb, + const ScalarT* ws, + ScalarT* f) + { + const auto VMEAS = static_cast(ReecbInternalVariables::VMEAS); + const auto PMEAS = static_cast(ReecbInternalVariables::PMEAS); + const auto XPIQ = static_cast(ReecbInternalVariables::XPIQ); + const auto XPIV = static_cast(ReecbInternalVariables::XPIV); + const auto QV = static_cast(ReecbInternalVariables::QV); + const auto PORD = static_cast(ReecbInternalVariables::PORD); + const auto VT = static_cast(ReecbInternalVariables::VT); + const auto VMEASSAFE = static_cast(ReecbInternalVariables::VMEASSAFE); + const auto SDIP = static_cast(ReecbInternalVariables::SDIP); + const auto VERR = static_cast(ReecbInternalVariables::VERR); + const auto IQV = static_cast(ReecbInternalVariables::IQV); + const auto QREF = static_cast(ReecbInternalVariables::QREF); + const auto EQ = static_cast(ReecbInternalVariables::EQ); + const auto VPIQ = static_cast(ReecbInternalVariables::VPIQ); + const auto EPIV = static_cast(ReecbInternalVariables::EPIV); + const auto FPORD = static_cast(ReecbInternalVariables::FPORD); + const auto RPORD = static_cast(ReecbInternalVariables::RPORD); + const auto IQCIRC = static_cast(ReecbInternalVariables::IQCIRC); + const auto IPCIRC = static_cast(ReecbInternalVariables::IPCIRC); + const auto IQMAX = static_cast(ReecbInternalVariables::IQMAX); + const auto IPMAX = static_cast(ReecbInternalVariables::IPMAX); + const auto IQBASE = static_cast(ReecbInternalVariables::IQBASE); + const auto IQRAW = static_cast(ReecbInternalVariables::IQRAW); + const auto IQCMD = static_cast(ReecbInternalVariables::IQCMD); + const auto IPCMD = static_cast(ReecbInternalVariables::IPCMD); + + const auto PE = static_cast(ReecbExternalVariables::PE); + const auto QGEN = static_cast(ReecbExternalVariables::QGEN); + const auto QEXT = static_cast(ReecbExternalVariables::QEXT); + const auto PFAREF = static_cast(ReecbExternalVariables::PFAREF); + const auto PREF = static_cast(ReecbExternalVariables::PREF); + + const ScalarT vr = wb[0]; + const ScalarT vi = wb[1]; + + const ScalarT pe = toComponentBase(ws[PE]); + const ScalarT qgen = toComponentBase(ws[QGEN]); + const ScalarT qext = toComponentBase(ws[QEXT]); + const ScalarT pfaref = ws[PFAREF]; + const ScalarT pref = toComponentBase(ws[PREF]); + const ScalarT iqcmd = toComponentBase(y[IQCMD]); + const ScalarT ipcmd = toComponentBase(y[IPCMD]); + + f[VMEAS] = -yp[VMEAS] + (y[VT] - y[VMEAS]) / Trv_eff_; + f[PMEAS] = -yp[PMEAS] + (pe - y[PMEAS]) / Tp_eff_; + f[XPIQ] = -yp[XPIQ] + y[SDIP] * Math::antiwindup(Kqp_ * y[EQ] + y[XPIQ], Kqi_ * y[EQ], Vmin_, Vmax_); + f[XPIV] = -yp[XPIV] + y[SDIP] * Math::antiwindup(Kvp_ * y[EPIV] + y[XPIV], Kvi_ * y[EPIV], -y[IQMAX], y[IQMAX]); + f[QV] = -yp[QV] + y[SDIP] * (y[QREF] / y[VMEASSAFE] - y[QV]) / Tiq_; + f[PORD] = -yp[PORD] + y[SDIP] * Math::antiwindup(y[PORD], y[RPORD], Pmin_, Pmax_); + f[VT] = -y[VT] * y[VT] + vr * vr + vi * vi; + f[VMEASSAFE] = -y[VMEASSAFE] + Math::max(y[VMEAS], VMEAS_MINIMUM); + f[SDIP] = -y[SDIP] + Math::inside(y[VT], Vdip_, Vup_); + f[VERR] = -y[VERR] + Math::deadband2(Vref0_ - y[VMEAS], dbd1_, dbd2_); + f[IQV] = -y[IQV] + Math::clamp(kqv_ * y[VERR], Iql1_, Iqh1_); + f[QREF] = -y[QREF] + PfFlag_ * y[PMEAS] * std::tan(pfaref) + pf_off_ * qext; + f[EQ] = -y[EQ] + Math::clamp(y[QREF], Qmin_, Qmax_) - qgen; + f[VPIQ] = -y[VPIQ] + Math::clamp(Kqp_ * y[EQ] + y[XPIQ], Vmin_, Vmax_); + f[EPIV] = -y[EPIV] + VFlag_ * y[VPIQ] + v_off_ * y[QREF] - y[VMEAS]; + f[FPORD] = -y[FPORD] + (pref - y[PORD]) / Tpord_; + f[RPORD] = -y[RPORD] + Math::clamp(y[FPORD], dPmin_, dPmax_); + f[IQCIRC] = -y[IQCIRC] * y[IQCIRC] + Imax_ * Imax_ - Pqflag_ * ipcmd * ipcmd; + f[IPCIRC] = -y[IPCIRC] * y[IPCIRC] + Imax_ * Imax_ - (ONE - Pqflag_) * iqcmd * iqcmd; + f[IQMAX] = -y[IQMAX] + (ONE - Pqflag_) * Imax_ + Pqflag_ * y[IQCIRC]; + f[IPMAX] = -y[IPMAX] + Pqflag_ * Imax_ + (ONE - Pqflag_) * y[IPCIRC]; + f[IQBASE] = -y[IQBASE] + Math::clamp(Kvp_ * y[EPIV] + y[XPIV], -y[IQMAX], y[IQMAX]); + f[IQRAW] = -y[IQRAW] + QFlag_ * y[IQBASE] + q_off_ * y[QV] + (ONE - y[SDIP]) * y[IQV]; + f[IQCMD] = -y[IQCMD] + toSystemBase(Math::clamp(y[IQRAW], -y[IQMAX], y[IQMAX])); + f[IPCMD] = -y[IPCMD] + toSystemBase(Math::clamp(y[PORD] / y[VMEASSAFE], ZERO, y[IPMAX])); + + return 0; + } + + template + int Reecb::evaluateResidual() + { + const auto PE = static_cast(ReecbExternalVariables::PE); + const auto QGEN = static_cast(ReecbExternalVariables::QGEN); + const auto QEXT = static_cast(ReecbExternalVariables::QEXT); + const auto PFAREF = static_cast(ReecbExternalVariables::PFAREF); + const auto PREF = static_cast(ReecbExternalVariables::PREF); + + ws_[PE] = pe_set_; + ws_[QGEN] = qgen_set_; + ws_[QEXT] = qext_set_; + ws_[PFAREF] = pfaref_set_; + ws_[PREF] = pref_set_; + std::fill(ws_indices_.begin(), ws_indices_.end(), INVALID_INDEX); + + if (signals_.template isAttached()) + { + ws_[PE] = signals_.template readExternalVariable(); + ws_indices_[PE] = signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[QGEN] = signals_.template readExternalVariable(); + ws_indices_[QGEN] = signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[QEXT] = signals_.template readExternalVariable(); + ws_indices_[QEXT] = signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[PFAREF] = signals_.template readExternalVariable(); + ws_indices_[PFAREF] = signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[PREF] = signals_.template readExternalVariable(); + ws_indices_[PREF] = signals_.template readExternalVariableIndex(); + } + + wb_[0] = Vr(); + wb_[1] = Vi(); + + evaluateInternalResidual(y_.data(), yp_.data(), wb_.data(), ws_.data(), f_.data()); + return 0; + } + } // namespace Converter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md index 5bff41ad4..246473672 100644 --- a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md +++ b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md @@ -144,16 +144,17 @@ are specified: `Genrou` | 6th order machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Tdop`, `Tdopp`, `Tqop`, `Tqopp`, `Xd`, `Xdp`, `Xdpp`, `Xq`, `Xqp`, `Xqpp`, `Xl`, `S10`, `S12`, `mva` | `ir`, `ii`, `p`, `q`, `delta`, `omega`, `speed` `Gensal` | 5th order salient-pole machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Tdop`, `Tdopp`, `Tqopp`, `Xd`, `Xdp`, `Xdpp`, `Xq`, `Xl`, `S10`, `S12`, `mva` | `ir`, `ii`, `p`, `q`, `delta`, `omega`, `speed`, `Eqp`, `psidp`, `psiqpp`, `psidpp`, `vd`, `vq`, `te`, `id`, `iq` `GenClassical` | the classical machine model | `bus`, `pmech`\*, `speed`\*, `efd`\* | `p0`, `q0`, `H`, `D`, `Ra`, `Xdp`, `mva` | `ir`, `ii`, `p`, `q`, `delta`, `omega` + `Reecb` | the REECB renewable electrical-control model | `bus`, `pe`\*, `qgen`\*, `qext`\*, `pfaref`\*, `pref`\*, `iqcmd`\*, `ipcmd`\* | `mva`, `PfFlag`, `VFlag`, `QFlag`, `Pqflag`, `Trv`, `Tp`, `Vref0`, `Vdip`, `Vup`, `dbd1`, `dbd2`, `kqv`, `Iql1`, `Iqh1`, `Qmax`, `Qmin`, `Kqp`, `Kqi`, `Vmax`, `Vmin`, `Kvp`, `Kvi`, `Tiq`, `Tpord`, `dPmax`, `dPmin`, `Pmax`, `Pmin`, `Imax` | `iqcmd`, `ipcmd`, `vmeas`, `pmeas` `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..054303e18 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 ReecbDataT = Converter::ReecbData; 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 reecb; ///< REECB electrical controllers within the model std::vector genrou; ///< GENROU instances within the model std::vector gensal; ///< GENSAL instances within the model std::vector genclassical; ///< Classical generator instances within the model diff --git a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp index a3682feca..d263e92d5 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp @@ -136,6 +136,12 @@ namespace GridKit raw_component.get_to(loadzip); sm.loadzip.push_back(loadzip); } + else if (kind == "Reecb") + { + typename SystemModelData::ReecbDataT reecb; + raw_component.get_to(reecb); + sm.reecb.push_back(reecb); + } else if (kind == "Tgov1") { typename SystemModelData::Tgov1DataT gov; diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index 287b328dc..bc82dcfd0 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,63 @@ namespace GridKit addComponent(adapter); } + // Add REECB electrical controllers + for (const auto& reecbdata : data.reecb) + { + IdxT bus_index = 0; + if (reecbdata.ports.contains(ReecbPorts::bus)) + { + bus_index = reecbdata.ports.at(ReecbPorts::bus); + } + + auto* reecb = new Reecb(getBus(bus_index), reecbdata); + + if (reecbdata.ports.contains(ReecbPorts::pe)) + { + const IdxT signal = reecbdata.ports.at(ReecbPorts::pe); + constexpr auto PE = ReecbExternalVariables::PE; + reecb->getSignals().template attachSignalNode(getSignal(signal)); + } + if (reecbdata.ports.contains(ReecbPorts::qgen)) + { + const IdxT signal = reecbdata.ports.at(ReecbPorts::qgen); + constexpr auto QGEN = ReecbExternalVariables::QGEN; + reecb->getSignals().template attachSignalNode(getSignal(signal)); + } + if (reecbdata.ports.contains(ReecbPorts::qext)) + { + const IdxT signal = reecbdata.ports.at(ReecbPorts::qext); + constexpr auto QEXT = ReecbExternalVariables::QEXT; + reecb->getSignals().template attachSignalNode(getSignal(signal)); + } + if (reecbdata.ports.contains(ReecbPorts::pfaref)) + { + const IdxT signal = reecbdata.ports.at(ReecbPorts::pfaref); + constexpr auto PFAREF = ReecbExternalVariables::PFAREF; + reecb->getSignals().template attachSignalNode(getSignal(signal)); + } + if (reecbdata.ports.contains(ReecbPorts::pref)) + { + const IdxT signal = reecbdata.ports.at(ReecbPorts::pref); + constexpr auto PREF = ReecbExternalVariables::PREF; + reecb->getSignals().template attachSignalNode(getSignal(signal)); + } + if (reecbdata.ports.contains(ReecbPorts::iqcmd)) + { + const IdxT signal = reecbdata.ports.at(ReecbPorts::iqcmd); + constexpr auto IQCMD = ReecbInternalVariables::IQCMD; + reecb->getSignals().template assignSignalNode(getSignal(signal)); + } + if (reecbdata.ports.contains(ReecbPorts::ipcmd)) + { + const IdxT signal = reecbdata.ports.at(ReecbPorts::ipcmd); + constexpr auto IPCMD = ReecbInternalVariables::IPCMD; + reecb->getSignals().template assignSignalNode(getSignal(signal)); + } + + addComponent(reecb); + } + // Add branches for (const auto& branchdata : data.branch) { diff --git a/docs/Figures/PhasorDynamics/REECB_diagram.png b/docs/Figures/PhasorDynamics/REECB_diagram.png new file mode 100644 index 000000000..c3aaa1478 Binary files /dev/null and b/docs/Figures/PhasorDynamics/REECB_diagram.png differ diff --git a/docs/GridKit/Model/PhasorDynamics/Converter/README.md b/docs/GridKit/Model/PhasorDynamics/Converter/README.md index fa56e8dd8..e868967e5 100644 --- a/docs/GridKit/Model/PhasorDynamics/Converter/README.md +++ b/docs/GridKit/Model/PhasorDynamics/Converter/README.md @@ -8,6 +8,7 @@ REGCA REGCB REECA +REECB ``` ```{include} ../../../../../GridKit/Model/PhasorDynamics/Converter/README.md diff --git a/docs/GridKit/Model/PhasorDynamics/Converter/REECB/README.md b/docs/GridKit/Model/PhasorDynamics/Converter/REECB/README.md new file mode 100644 index 000000000..ff972cda5 --- /dev/null +++ b/docs/GridKit/Model/PhasorDynamics/Converter/REECB/README.md @@ -0,0 +1,6 @@ +# REECB + +```{include} ../../../../../../GridKit/Model/PhasorDynamics/Converter/REECB/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..13e3b4d03 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_reecb runConverterReecbTests.cpp) +target_link_libraries( + test_phasor_converter_reecb + 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 PhasorDynamicsConverterReecbTest COMMAND test_phasor_converter_reecb) add_test(NAME PhasorDynamicsStabilizerIeeestTest COMMAND test_phasor_stabilizer_ieeest) add_test(NAME PhasorDynamicsGenClassicalTest COMMAND test_phasor_gen_classical) add_test(NAME PhasorDynamicsLoadZTest COMMAND test_phasor_loadz) @@ -171,6 +180,7 @@ install( test_phasor_exciter_ieeet1 test_phasor_gensal test_phasor_exciter_sexspti + test_phasor_converter_reecb test_phasor_stabilizer_ieeest test_phasor_gen_classical test_phasor_system diff --git a/tests/UnitTests/PhasorDynamics/ConverterReecbTests.hpp b/tests/UnitTests/PhasorDynamics/ConverterReecbTests.hpp new file mode 100644 index 000000000..faf2fd13d --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/ConverterReecbTests.hpp @@ -0,0 +1,641 @@ +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + template + class ConverterReecbTests + { + public: + using ScalarT = scalar_type; + using IdxT = index_type; + using RealT = typename PhasorDynamics::Component::RealT; + using ReecbT = PhasorDynamics::Converter::Reecb; + using Var = PhasorDynamics::Converter::ReecbInternalVariables; + using Ext = PhasorDynamics::Converter::ReecbExternalVariables; + using Params = PhasorDynamics::Converter::ReecbParameters; + + static constexpr ScalarT kTol = static_cast(1.0e-8); + + TestOutcome validation() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + + ReecbT reecb(&bus, makeData()); + success *= (reecb.size() == static_cast(Var::MAXIMUM)); + success *= (reecb.getMonitor() != nullptr); + success *= (reecb.verify() == 0); + + auto minimal = makeMinimalData(); + minimal.parameters[Params::mva] = static_cast(100.0); + ReecbT minimal_model(&bus, minimal); + success *= (minimal_model.verify() == 0); + + ReecbT missing_mva(&bus, makeMinimalData()); + success *= (missing_mva.verify() > 0); + + auto bad_band = makeData(); + bad_band.parameters[Params::Vdip] = static_cast(1.2); + bad_band.parameters[Params::Vup] = static_cast(1.2); + ReecbT bad_band_model(&bus, bad_band); + success *= (bad_band_model.verify() > 0); + + auto bad_imax = makeData(); + bad_imax.parameters[Params::Imax] = static_cast(-1.0); + ReecbT bad_imax_model(&bus, bad_imax); + success *= (bad_imax_model.verify() > 0); + + ScalarT pe_value{0.5}; + IdxT pe_index = 20; + PhasorDynamics::SignalNode pe_node; + pe_node.set(&pe_value, &pe_index); + + ReecbT half_connected(&bus, makeData()); + half_connected.getSignals().template attachSignalNode(&pe_node); + success *= (half_connected.verify() == 0); + + PhasorDynamics::SignalNode unlinked_pe_node; + ReecbT unlinked(&bus, makeData()); + unlinked.getSignals().template attachSignalNode(&unlinked_pe_node); + success *= (unlinked.verify() > 0); + + return success.report(__func__); + } + + TestOutcome signals() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + ScalarT iqcmd_value{0.2}; + ScalarT ipcmd_value{0.6}; + IdxT iqcmd_index = 21; + IdxT ipcmd_index = 22; + + PhasorDynamics::SignalNode iqcmd_node; + PhasorDynamics::SignalNode ipcmd_node; + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + + ReecbT reecb(&bus, makeData()); + reecb.getSignals().template assignSignalNode(&iqcmd_node); + reecb.getSignals().template assignSignalNode(&ipcmd_node); + + success *= (reecb.allocate() == 0); + iqcmd_node.init(static_cast(0.2)); + ipcmd_node.init(static_cast(0.6)); + success *= (reecb.verify() == 0); + success *= (reecb.initialize() == 0); + success *= (reecb.tagDifferentiable() == 0); + success *= (reecb.evaluateResidual() == 0); + + success *= isEqual(reecb.y()[index(Var::VMEAS)], static_cast(1.0), kTol); + success *= isEqual(reecb.y()[index(Var::PMEAS)], static_cast(0.6), kTol); + success *= isEqual(reecb.y()[index(Var::QREF)], static_cast(0.2), kTol); + success *= isEqual(reecb.y()[index(Var::PORD)], static_cast(0.6), kTol); + success *= isEqual(iqcmd_node.read(), reecb.y()[index(Var::IQCMD)], kTol); + success *= isEqual(ipcmd_node.read(), reecb.y()[index(Var::IPCMD)], kTol); + success *= (reecb.tag()[index(Var::VMEAS)] == true); + success *= (reecb.tag()[index(Var::PMEAS)] == true); + success *= allZero(reecb); + + return success.report(__func__); + } + + TestOutcome publishRefs() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(0.8, 0.6); + bus.allocate(); + bus.initialize(); + + auto data = makeData(); + data.parameters[Params::Qmin] = static_cast(-2.0); + data.parameters[Params::Qmax] = static_cast(2.0); + data.parameters[Params::Pmax] = static_cast(2.0); + data.parameters[Params::Imax] = static_cast(2.0); + + ScalarT iqcmd_value{-0.2}; + ScalarT ipcmd_value{1.6}; + ScalarT qext_value{99.0}; + ScalarT pfaref_value{99.0}; + ScalarT pref_value{99.0}; + IdxT iqcmd_index = 30; + IdxT ipcmd_index = 31; + IdxT qext_index = 32; + IdxT pfref_index = 33; + IdxT pref_index = 34; + + PhasorDynamics::SignalNode iqcmd_node; + PhasorDynamics::SignalNode ipcmd_node; + PhasorDynamics::SignalNode qext_node; + PhasorDynamics::SignalNode pfaref_node; + PhasorDynamics::SignalNode pref_node; + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + qext_node.set(&qext_value, &qext_index); + pfaref_node.set(&pfaref_value, &pfref_index); + pref_node.set(&pref_value, &pref_index); + + ReecbT reecb(&bus, data); + reecb.getSignals().template assignSignalNode(&iqcmd_node); + reecb.getSignals().template assignSignalNode(&ipcmd_node); + reecb.getSignals().template attachSignalNode(&qext_node); + reecb.getSignals().template attachSignalNode(&pfaref_node); + reecb.getSignals().template attachSignalNode(&pref_node); + + success *= (reecb.allocate() == 0); + iqcmd_node.init(static_cast(-0.2)); + ipcmd_node.init(static_cast(1.6)); + success *= (reecb.verify() == 0); + success *= (reecb.initialize() == 0); + success *= (reecb.evaluateResidual() == 0); + + const ScalarT expected_pfaref = static_cast(std::atan(static_cast(-0.2 / 1.6))); + success *= isEqual(qext_node.read(), static_cast(-0.2), kTol); + success *= isEqual(pfaref_node.read(), expected_pfaref, kTol); + success *= isEqual(pref_node.read(), static_cast(1.6), kTol); + success *= isEqual(reecb.y()[index(Var::QREF)], static_cast(-0.2), kTol); + success *= allZero(reecb); + + return success.report(__func__); + } + + TestOutcome baseSignals() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + auto data = makeData(); + data.parameters[Params::mva] = static_cast(50.0); + + ScalarT pe_value{99.0}; + ScalarT qgen_value{99.0}; + ScalarT qext_value{99.0}; + ScalarT pfaref_value{99.0}; + ScalarT pref_value{99.0}; + ScalarT iqcmd_value{0.0}; + ScalarT ipcmd_value{0.0}; + IdxT pe_index = 40; + IdxT qgen_index = 41; + IdxT qext_index = 42; + IdxT pfaref_index = 43; + IdxT pref_index = 44; + IdxT iqcmd_index = 45; + IdxT ipcmd_index = 46; + + PhasorDynamics::SignalNode pe_node; + PhasorDynamics::SignalNode qgen_node; + PhasorDynamics::SignalNode qext_node; + PhasorDynamics::SignalNode pfaref_node; + PhasorDynamics::SignalNode pref_node; + PhasorDynamics::SignalNode iqcmd_node; + PhasorDynamics::SignalNode ipcmd_node; + pe_node.set(&pe_value, &pe_index); + qgen_node.set(&qgen_value, &qgen_index); + qext_node.set(&qext_value, &qext_index); + pfaref_node.set(&pfaref_value, &pfaref_index); + pref_node.set(&pref_value, &pref_index); + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + + ReecbT reecb(&bus, data); + reecb.getSignals().template attachSignalNode(&pe_node); + reecb.getSignals().template attachSignalNode(&qgen_node); + reecb.getSignals().template attachSignalNode(&qext_node); + reecb.getSignals().template attachSignalNode(&pfaref_node); + reecb.getSignals().template attachSignalNode(&pref_node); + reecb.getSignals().template assignSignalNode(&iqcmd_node); + reecb.getSignals().template assignSignalNode(&ipcmd_node); + + success *= (reecb.allocate() == 0); + iqcmd_node.init(static_cast(0.05)); + ipcmd_node.init(static_cast(0.25)); + success *= (reecb.verify() == 0); + success *= (reecb.initialize() == 0); + success *= (reecb.evaluateResidual() == 0); + + const ScalarT expected_pfaref = static_cast(std::atan(static_cast(0.1 / 0.5))); + success *= isEqual(reecb.y()[index(Var::PMEAS)], static_cast(0.5), kTol); + success *= isEqual(reecb.y()[index(Var::QREF)], static_cast(0.1), kTol); + success *= isEqual(reecb.y()[index(Var::PORD)], static_cast(0.5), kTol); + success *= isEqual(pe_node.read(), static_cast(0.25), kTol); + success *= isEqual(qgen_node.read(), static_cast(0.05), kTol); + success *= isEqual(qext_node.read(), static_cast(0.05), kTol); + success *= isEqual(pfaref_node.read(), expected_pfaref, kTol); + success *= isEqual(pref_node.read(), static_cast(0.25), kTol); + success *= isEqual(iqcmd_node.read(), static_cast(0.05), kTol); + success *= isEqual(ipcmd_node.read(), static_cast(0.25), kTol); + success *= isEqual(reecb.y()[index(Var::IQCMD)], static_cast(0.05), kTol); + success *= isEqual(reecb.y()[index(Var::IPCMD)], static_cast(0.25), kTol); + success *= allZero(reecb); + + return success.report(__func__); + } + + TestOutcome feedbackBase() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + auto data = makeData(); + data.parameters[Params::mva] = static_cast(50.0); + + ScalarT pe_value{0.25}; + ScalarT qgen_value{0.05}; + ScalarT iqcmd_value{0.0}; + ScalarT ipcmd_value{0.0}; + IdxT pe_index = 40; + IdxT qgen_index = 41; + IdxT iqcmd_index = 42; + IdxT ipcmd_index = 43; + + PhasorDynamics::SignalNode pe_node; + PhasorDynamics::SignalNode qgen_node; + PhasorDynamics::SignalNode iqcmd_node; + PhasorDynamics::SignalNode ipcmd_node; + pe_node.set(&pe_value, &pe_index); + qgen_node.set(&qgen_value, &qgen_index); + iqcmd_node.set(&iqcmd_value, &iqcmd_index); + ipcmd_node.set(&ipcmd_value, &ipcmd_index); + + ReecbT reecb(&bus, data); + reecb.getSignals().template attachSignalNode(&pe_node); + reecb.getSignals().template attachSignalNode(&qgen_node); + reecb.getSignals().template assignSignalNode(&iqcmd_node); + reecb.getSignals().template assignSignalNode(&ipcmd_node); + + success *= (reecb.allocate() == 0); + iqcmd_node.init(static_cast(0.05)); + ipcmd_node.init(static_cast(0.25)); + success *= (reecb.verify() == 0); + success *= (reecb.initialize() == 0); + success *= (reecb.evaluateResidual() == 0); + + success *= isEqual(reecb.y()[index(Var::PMEAS)], static_cast(0.5), kTol); + success *= isEqual(reecb.y()[index(Var::QREF)], static_cast(0.1), kTol); + success *= isEqual(reecb.y()[index(Var::PORD)], static_cast(0.5), kTol); + success *= isEqual(pe_node.read(), static_cast(0.25), kTol); + success *= isEqual(qgen_node.read(), static_cast(0.05), kTol); + success *= isEqual(reecb.y()[index(Var::IQCMD)], static_cast(0.05), kTol); + success *= isEqual(reecb.y()[index(Var::IPCMD)], static_cast(0.25), kTol); + success *= allZero(reecb); + + return success.report(__func__); + } + + TestOutcome zeroTime() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + auto data = makeData(); + data.parameters[Params::Trv] = static_cast(0.0); + data.parameters[Params::Tp] = static_cast(0.0); + + ReecbT reecb(&bus, data); + success *= (reecb.allocate() == 0); + success *= (reecb.verify() == 0); + reecb.y()[index(Var::IQCMD)] = static_cast(0.2); + reecb.y()[index(Var::IPCMD)] = static_cast(0.6); + success *= (reecb.initialize() == 0); + success *= (reecb.tagDifferentiable() == 0); + success *= (reecb.tag()[index(Var::VMEAS)] == true); + success *= (reecb.tag()[index(Var::PMEAS)] == true); + + reecb.yp()[index(Var::VMEAS)] = static_cast(1.0); + reecb.yp()[index(Var::PMEAS)] = static_cast(2.0); + success *= (reecb.evaluateResidual() == 0); + success *= isEqual(reecb.getResidual()[index(Var::VMEAS)], static_cast(-1.0), kTol); + success *= isEqual(reecb.getResidual()[index(Var::PMEAS)], static_cast(-2.0), kTol); + + reecb.yp()[index(Var::VMEAS)] = ZERO; + reecb.y()[index(Var::VMEAS)] = static_cast(0.99); + success *= (reecb.evaluateResidual() == 0); + success *= isEqual(reecb.getResidual()[index(Var::VMEAS)], static_cast(10.0), kTol); + + return success.report(__func__); + } + + TestOutcome qPriority() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + auto data = makeData(); + data.parameters[Params::Pqflag] = static_cast(0); + data.parameters[Params::Imax] = static_cast(1.1); + + ReecbT reecb(&bus, data); + success *= (reecb.allocate() == 0); + success *= (reecb.verify() == 0); + reecb.y()[index(Var::IQCMD)] = static_cast(0.3); + reecb.y()[index(Var::IPCMD)] = static_cast(0.9); + success *= (reecb.initialize() == 0); + success *= (reecb.evaluateResidual() == 0); + + const ScalarT ipmax = std::sqrt(static_cast(1.1 * 1.1 - 0.3 * 0.3)); + success *= isEqual(reecb.y()[index(Var::IQMAX)], static_cast(1.1), kTol); + success *= isEqual(reecb.y()[index(Var::IPMAX)], ipmax, kTol); + success *= isEqual(reecb.y()[index(Var::IQCMD)], static_cast(0.3), kTol); + success *= isEqual(reecb.y()[index(Var::IPCMD)], static_cast(0.9), kTol); + success *= allZero(reecb); + + return success.report(__func__); + } + + TestOutcome pPriority() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + auto data = makeData(); + data.parameters[Params::Pqflag] = static_cast(1); + data.parameters[Params::Imax] = static_cast(1.1); + + ReecbT reecb(&bus, data); + success *= (reecb.allocate() == 0); + success *= (reecb.verify() == 0); + reecb.y()[index(Var::IQCMD)] = static_cast(0.3); + reecb.y()[index(Var::IPCMD)] = static_cast(0.9); + success *= (reecb.initialize() == 0); + success *= (reecb.evaluateResidual() == 0); + + const ScalarT iqmax = std::sqrt(static_cast(1.1 * 1.1 - 0.9 * 0.9)); + success *= isEqual(reecb.y()[index(Var::IPMAX)], static_cast(1.1), kTol); + success *= isEqual(reecb.y()[index(Var::IQMAX)], iqmax, kTol); + success *= isEqual(reecb.y()[index(Var::IQCMD)], static_cast(0.3), kTol); + success *= isEqual(reecb.y()[index(Var::IPCMD)], static_cast(0.9), kTol); + success *= allZero(reecb); + + return success.report(__func__); + } + + TestOutcome voltageBand() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(0.8, 0.0); + bus.allocate(); + bus.initialize(); + + auto data = makeData(); + data.parameters[Params::QFlag] = static_cast(1); + data.parameters[Params::Vref0] = static_cast(1.0); + data.parameters[Params::Vdip] = static_cast(0.9); + data.parameters[Params::Vup] = static_cast(1.1); + data.parameters[Params::dbd1] = static_cast(0.0); + data.parameters[Params::dbd2] = static_cast(0.0); + data.parameters[Params::kqv] = static_cast(1.0); + data.parameters[Params::Imax] = static_cast(2.0); + + ReecbT reecb(&bus, data); + success *= (reecb.allocate() == 0); + success *= (reecb.verify() == 0); + reecb.y()[index(Var::IQCMD)] = ZERO; + reecb.y()[index(Var::IPCMD)] = static_cast(0.5); + success *= (reecb.initialize() == 0); + success *= (reecb.evaluateResidual() == 0); + + const ScalarT expected_sdip = Math::inside(reecb.y()[index(Var::VT)], static_cast(0.9), static_cast(1.1)); + const ScalarT expected_iqraw = reecb.y()[index(Var::IQBASE)] + (ONE - reecb.y()[index(Var::SDIP)]) * reecb.y()[index(Var::IQV)]; + success *= isEqual(reecb.y()[index(Var::SDIP)], expected_sdip, kTol); + success *= isEqual(reecb.y()[index(Var::IQRAW)], expected_iqraw, kTol); + success *= (reecb.y()[index(Var::SDIP)] < static_cast(0.01)); + success *= (reecb.y()[index(Var::IQV)] > static_cast(0.1)); + success *= allZero(reecb); + + return success.report(__func__); + } + + TestOutcome piSaturation() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.13, 0.0); + bus.allocate(); + bus.initialize(); + + auto data = makeData(); + data.parameters[Params::QFlag] = static_cast(0); + data.parameters[Params::Pqflag] = static_cast(0); + data.parameters[Params::Vmin] = static_cast(0.9); + data.parameters[Params::Vmax] = static_cast(1.05); + data.parameters[Params::Kvp] = static_cast(10.0); + data.parameters[Params::Kvi] = static_cast(60.0); + data.parameters[Params::Vup] = static_cast(99.0); + data.parameters[Params::Vdip] = static_cast(-99.0); + data.parameters[Params::Imax] = static_cast(1.1); + + ReecbT reecb(&bus, data); + success *= (reecb.allocate() == 0); + success *= (reecb.verify() == 0); + reecb.y()[index(Var::IQCMD)] = static_cast(0.15 / 1.13); + reecb.y()[index(Var::IPCMD)] = static_cast(0.5 / 1.13); + success *= (reecb.initialize() == 0); + success *= (reecb.evaluateResidual() == 0); + + const ScalarT piv_arg = static_cast(10.0) * reecb.y()[index(Var::EPIV)] + reecb.y()[index(Var::XPIV)]; + success *= (piv_arg < -reecb.y()[index(Var::IQMAX)]); + success *= allZero(reecb); + + return success.report(__func__); + } + +#ifdef GRIDKIT_ENABLE_ENZYME + TestOutcome jacobian() + { + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + bus.allocate(); + bus.initialize(); + + ReecbT reecb(&bus, makeData()); + success *= (reecb.allocate() == 0); + success *= (reecb.verify() == 0); + reecb.y()[index(Var::IQCMD)] = static_cast(0.2); + reecb.y()[index(Var::IPCMD)] = static_cast(0.6); + success *= (reecb.initialize() == 0); + success *= (reecb.evaluateResidual() == 0); + success *= (reecb.evaluateJacobian() == 0); + + auto* jac = reecb.getCooJacobian(); + success *= (jac != nullptr); + if (jac != nullptr) + { + success *= (jac->getNnz() > 0); + const auto* values = jac->getValues(); + for (IdxT i = 0; i < jac->getNnz(); ++i) + { + success *= std::isfinite(values[i]); + } + } + + return success.report(__func__); + } +#endif + + TestOutcome json() + { + TestStatus success = true; + + std::istringstream input(R"json( +{ + "header": { + "format_version": 0, + "format_revision": 1, + "case_name": "renewable electrical control", + "case_description": "REECB parser test", + "case_comments": "", + "freq_base": 60.0, + "va_base": 100000000.0 + }, + "buses": [ + { "number": 1, "class": "bus", "name": "Bus 1", "init": { "Vr": 1.0, "Vi": 0.0 }, "v_base": 1.0 } + ], + "signals": [ + { "signal_id": 12, "name": "Iqcmd" }, + { "signal_id": 13, "name": "Ipcmd" } + ], + "devices": [ + { + "class": "Reecb", + "ports": { "bus": 1, "iqcmd": 12, "ipcmd": 13 }, + "id": "REE1", + "params": { + "mva": 100.0, "PfFlag": 0, "VFlag": 1, "QFlag": 0, "Pqflag": 1, + "Trv": 0.0, "Tp": 0.02, "Vdip": 0.7, "Vup": 1.2, + "dbd1": -0.01, "dbd2": 0.01, "kqv": 0.0, "Iql1": -1.0, "Iqh1": 1.0, + "Qmax": 1.0, "Qmin": -1.0, "Kqp": 1.0, "Kqi": 0.0, + "Vmax": 1.2, "Vmin": 0.8, "Kvp": 1.0, "Kvi": 0.0, + "Tiq": 0.02, "Tpord": 0.02, "dPmax": 1.0, "dPmin": -1.0, + "Pmax": 1.0, "Pmin": 0.0, "Imax": 2.0 + } + } + ] +} +)json"); + + auto data = PhasorDynamics::parseSystemModelData(input); + success *= (data.reecb.size() == 1); + success *= (std::get(data.reecb[0].parameters.at(Params::Pqflag)) == static_cast(1)); + success *= (std::get(data.reecb[0].parameters.at(Params::mva)) == static_cast(100.0)); + + PhasorDynamics::SystemModel system(data); + success *= (system.allocate() == 0); + success *= (system.initialize() == 0); + success *= (system.evaluateResidual() == 0); + success *= (system.size() == 27); + + return success.report(__func__); + } + + private: + static size_t index(Var variable) + { + return static_cast(variable); + } + + TestStatus allZero(const ReecbT& reecb) const + { + TestStatus success = true; + + for (size_t i = 0; i < reecb.getResidual().size(); ++i) + { + success *= isEqual(reecb.getResidual()[i], static_cast(0.0), kTol); + success *= isEqual(reecb.yp()[i], static_cast(0.0), kTol); + } + + return success; + } + + auto makeMinimalData() -> PhasorDynamics::Converter::ReecbData + { + using Mon = PhasorDynamics::Converter::ReecbMonitorableVariables; + + PhasorDynamics::Converter::ReecbData data; + data.device_class = "Reecb"; + data.disambiguation_string = "reecb_test"; + data.monitored_variables.insert(Mon::iqcmd); + data.monitored_variables.insert(Mon::ipcmd); + data.monitored_variables.insert(Mon::vmeas); + data.monitored_variables.insert(Mon::pmeas); + return data; + } + + auto makeData() -> PhasorDynamics::Converter::ReecbData + { + auto data = makeMinimalData(); + + data.parameters[Params::mva] = static_cast(100.0); + data.parameters[Params::PfFlag] = static_cast(0); + data.parameters[Params::VFlag] = static_cast(1); + data.parameters[Params::QFlag] = static_cast(0); + data.parameters[Params::Pqflag] = static_cast(1); + data.parameters[Params::Trv] = static_cast(0.0); + data.parameters[Params::Tp] = static_cast(0.02); + data.parameters[Params::Vdip] = static_cast(0.7); + data.parameters[Params::Vup] = static_cast(1.2); + data.parameters[Params::dbd1] = static_cast(-0.01); + data.parameters[Params::dbd2] = static_cast(0.01); + data.parameters[Params::kqv] = static_cast(0.0); + data.parameters[Params::Iql1] = static_cast(-1.0); + data.parameters[Params::Iqh1] = static_cast(1.0); + data.parameters[Params::Qmax] = static_cast(1.0); + data.parameters[Params::Qmin] = static_cast(-1.0); + data.parameters[Params::Kqp] = static_cast(1.0); + data.parameters[Params::Kqi] = static_cast(0.0); + data.parameters[Params::Vmax] = static_cast(1.2); + data.parameters[Params::Vmin] = static_cast(0.8); + data.parameters[Params::Kvp] = static_cast(1.0); + data.parameters[Params::Kvi] = static_cast(0.0); + data.parameters[Params::Tiq] = static_cast(0.02); + data.parameters[Params::Tpord] = static_cast(0.02); + data.parameters[Params::dPmax] = static_cast(1.0); + data.parameters[Params::dPmin] = static_cast(-1.0); + data.parameters[Params::Pmax] = static_cast(1.0); + data.parameters[Params::Pmin] = static_cast(0.0); + data.parameters[Params::Imax] = static_cast(2.0); + + return data; + } + }; + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/runConverterReecbTests.cpp b/tests/UnitTests/PhasorDynamics/runConverterReecbTests.cpp new file mode 100644 index 000000000..90c935e34 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/runConverterReecbTests.cpp @@ -0,0 +1,25 @@ +#include "ConverterReecbTests.hpp" + +int main() +{ + GridKit::Testing::TestingResults result; + + GridKit::Testing::ConverterReecbTests test; + + result += test.validation(); + result += test.signals(); + result += test.publishRefs(); + result += test.baseSignals(); + result += test.feedbackBase(); + result += test.zeroTime(); + result += test.qPriority(); + result += test.pPriority(); + result += test.voltageBand(); + result += test.piSaturation(); +#ifdef GRIDKIT_ENABLE_ENZYME + result += test.jacobian(); +#endif + result += test.json(); + + return result.summary(); +}