diff --git a/CHANGELOG.md b/CHANGELOG.md index d6076d9f..161b37c3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -38,6 +38,7 @@ - Updated Jacobian value storage from `ScalarT` to `RealT`. - Added a header file defining constants to be used throughout the code. - Added `GridKitDocs` target for Doxygen documentation. +- Added a header file defining common math functions (e.g., sigmoid) to be used throughout the code. ## v0.1 diff --git a/GridKit/CMakeLists.txt b/GridKit/CMakeLists.txt index 2a57de7e..a01f8a6c 100644 --- a/GridKit/CMakeLists.txt +++ b/GridKit/CMakeLists.txt @@ -28,6 +28,7 @@ add_subdirectory(Solver) install( FILES Constants.hpp + CommonMath.hpp ScalarTraits.hpp DESTINATION include/GridKit) diff --git a/GridKit/CommonMath.hpp b/GridKit/CommonMath.hpp new file mode 100644 index 00000000..7866d8f1 --- /dev/null +++ b/GridKit/CommonMath.hpp @@ -0,0 +1,96 @@ +#pragma once + +#include + +#include +#include + +namespace GridKit +{ + namespace Math + { + /** + * @brief Scaled sigmoid activation function + * + * @note The sigmoid constant (mu) value is chosen to balance accuracy + * and finite derivatives. Large values more closely approximate a step + * function, but lead to inf or NaN derivatives. + * The value of 240 corresponds to a time step of 1/4 of a 60Hz cycle. + * From experiments, numerical derivatives break down for mu above 700. + * + * @tparam ScalarT - scalar data type + * + * @param[in] x + * @return value of the sigmoid function + */ + template + __attribute__((always_inline)) inline ScalarT sigmoid(const ScalarT x) + { + using RealT = typename GridKit::ScalarTraits::RealT; + static constexpr RealT MU = 240.0; + return ONE / (ONE + std::exp(-MU * x)); + } + + /** + * @brief Low indicator function for regulator limits + * + * @tparam ScalarT - Scalar data type + * @tparam RealT - Real data type (see GridKit::ScalarTraits::RealT) + * + * @param[in] limit_min - Minimum limit + * @param[in] x - State variable + * @param[in] f - Conditional derivative of state variable + * @return Scalar value indicating limit activation + */ + template + __attribute__((always_inline)) inline ScalarT indicator_low( + const RealT limit_min, + const ScalarT x, + const ScalarT f) + { + return sigmoid(limit_min - x) * sigmoid(-f); + } + + /** + * @brief High indicator function for regulator limits + * + * @tparam ScalarT - Scalar data type + * @tparam RealT - Real data type (see GridKit::ScalarTraits::RealT) + * + * @param[in] limit_max - Maximum limit + * @param[in] x - State variable + * @param[in] f - Conditional derivative of state variable + * @return Scalar value indicating limit activation + */ + template + __attribute__((always_inline)) inline ScalarT indicator_high( + const RealT limit_max, + const ScalarT x, + const ScalarT f) + { + return sigmoid(x - limit_max) * sigmoid(f); + } + + /** + * @brief Net Indicator function for regulator limits + * + * @tparam ScalarT - Scalar data type + * @tparam RealT - Real data type (see GridKit::ScalarTraits::RealT) + * + * @param[in] limit_min - Minimum limit + * @param[in] limit_max - Maximum limit + * @param[in] x - State variable + * @param[in] f - Conditional derivative of state variable + * @return Scalar value indicating limit activation + */ + template + __attribute__((always_inline)) inline ScalarT indicator( + const RealT limit_min, + const RealT limit_max, + const ScalarT x, + const ScalarT f) + { + return (ONE - indicator_low(limit_min, x, f)) * (ONE - indicator_high(limit_max, x, f)); + } + } // namespace Math +} // namespace GridKit diff --git a/GridKit/Model/Evaluator.hpp b/GridKit/Model/Evaluator.hpp index 98789eba..12b760c6 100644 --- a/GridKit/Model/Evaluator.hpp +++ b/GridKit/Model/Evaluator.hpp @@ -2,6 +2,7 @@ #include +#include #include #include #include diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp index 33fa4d98..ff20ca4c 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1.hpp @@ -150,15 +150,6 @@ namespace GridKit /// Component signal extension ComponentSignals signals_; - // Scale of Sigmoid function - // (temporary local implementation) - // This value gave higher precision. - const RealT mu_{400000.0}; - - // Activation function (sigmoid approximation) - ScalarT sigmoid(ScalarT x); - ScalarT indicator(ScalarT x, ScalarT f); - // Parameter initialization function void initModelParams(const model_data_type& data); }; diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp index 1623c17c..17640d4a 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/Ieeet1Impl.hpp @@ -236,48 +236,6 @@ namespace GridKit return 0; } - /** - * @brief Scaled sigmoid activation function - * - * @param[in] x State variable - * @tparam ScalarT Scalar data type - * @tparam IdxT Index data type - * @return Sigmoid approximation evaluated at x - * - * @warning This needs to be abstracted to be used - * across phasor dynamics. Identical pattern is - * being used in TGOV1 model. - * - * Algebraic approximation of transcendental sigmoid. - */ - template - ScalarT Ieeet1::sigmoid(ScalarT x) - { - return ((HALF * mu_ * x) / (ONE + std::abs(mu_ * x))) + HALF; - } - - /** - * @brief Net Indicator function for regulator limits - * - * @param[in] x State variable - * @param[in] f Conditional derivative of state variable - * @tparam ScalarT Scalar data type - * @tparam IdxT Index data type - * @return Scalar value indicating limit activation. - * - * @warning This needs to be abstracted to be used - * across phasor dynamics. Identical pattern is - * being used in TGOV1 model. - */ - template - ScalarT Ieeet1::indicator(ScalarT x, ScalarT f) - { - - ScalarT ind_low = (this->sigmoid(Vrmin_ - x)) * (this->sigmoid(-f)); - ScalarT ind_high = (this->sigmoid(x - Vrmax_)) * (this->sigmoid(f)); - return (ONE - ind_low) * (ONE - ind_high); - } - /** * @brief Residuals of system equations * @@ -326,7 +284,7 @@ namespace GridKit // The 'pre-limit' derivative of Pv ScalarT f = -vr + Ka_ * vtr; - ScalarT vr_ind = this->indicator(vr, f); + ScalarT vr_ind = Math::indicator(Vrmin_, Vrmax_, vr, f); // Internal Differential Equations f_[0] = -vts_dot + (Ec_ - vts) / Tr_; @@ -340,11 +298,8 @@ namespace GridKit f_[6] = -ve + ksat * efdp; f_[7] = -efd + efdp + omega * efdp * Ispdlim_; - // TODO smooth approaximation for Enzyme - // NOTE seems about double PW saturation. - f_[8] = -ksat; - if (efdp > SA_) - f_[8] += SB_ * (efdp - SA_) * (efdp - SA_); + ScalarT efd_sat = (efdp - SA_) * (Math::sigmoid(efdp - SA_)); + f_[8] = -ksat + SB_ * efd_sat * efd_sat; return 0; } diff --git a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md index 5f4422f3..64aa85e3 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md +++ b/GridKit/Model/PhasorDynamics/Exciter/IEEET1/README.md @@ -155,58 +155,11 @@ The indicator function $\phi$ can be defined in terms of a scaled activation fun \end{aligned} ``` -The scale of the sigmoid function ($\alpha$ on the order of $10^3$) should be chosen so that for all practical parameters of the IEEET1 model, the sigmoid acts as a step function. This is further approximated by an algebraic form to obtain a practical function during implementation. +The scale of the sigmoid function ($\alpha$ on the order of $10^3$) should be chosen so that for all practical parameters of the IEEET1 model, the sigmoid acts as a step function. ```math \begin{aligned} \sigma(x) = \dfrac{1}{1+\exp(-\alpha x)} - \approx - \dfrac{1}{2} - \left( - \dfrac{\alpha x}{1+|\alpha x|} + 1 - \right) -\end{aligned} -``` - -Applying these approximations produces a smooth approximation of the explicit differential equation of $V_R$. The approximation written out below approaches an exact solution as $\alpha\to\infty$. - - -```math -\begin{aligned} - \dot{V}_R - &\approx f - \left[ - 1 + - \dfrac{\alpha^2}{4} - \dfrac{V_{rmin}\text{-}V_R} - {1+\alpha |V_{rmin}\text{-}V_R|} - \dfrac{f}{1+\alpha |f|} - \right] - \left[ - 1- - \dfrac{\alpha^2}{4} - \dfrac{V_R\text{-}V_{rmax}} - {1+\alpha| V_R\text{-}V_{rmax}|} - \dfrac{f}{1+\alpha |f|} - \right]\\ \end{aligned} ``` @@ -228,23 +181,20 @@ The algebraic equations of the exciter. \end{cases} \\ \end{aligned} ``` - #### Smooth Piecewise Approximation (Algebraic) For the algebraic piecewise functions (non-flags), this implementation is straightforward when the approximation above is used. ```math \begin{aligned} + E_{fd} + &=(1 + \omega I_{spdlm})E_{fd}' \\ k_{sat} - &=\sigma (E_{fd}' -S_A) \cdot S_B(E_{fd}' -S_A)^2 -\end{aligned} -``` -The approximation written out below approaches an exact solution as $\alpha\to\infty$. -```math -\begin{aligned} - k_{sat}&\approx\dfrac{\alpha S_B}{2} \dfrac{(E_{fd}' -S_A)^3}{1+\alpha|E_{fd}' -S_A|} + &=S_B\left[(E_{fd}' -S_A) \cdot \sigma (E_{fd}' -S_A)\right]^2 \end{aligned} ``` +The approximation approaches an exact solution as $\alpha\to\infty$. + ## Initialization At steady state we assume that $V_R$ is at or within its limits, diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/README.md b/GridKit/Model/PhasorDynamics/Governor/Tgov1/README.md index 27431fbf..b4d695a6 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/README.md +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/README.md @@ -101,10 +101,10 @@ The indicator function $\phi$ can be defined in terms of a scaled activation fun \end{aligned} ``` -The scale of the sigmoid function ($\alpha$ on the order of $10^3$) should be chosen so that for all practical parameters of the TGOV1 model, the sigmoid acts as a step function. This is further approximated by an algebraic form to obtain a practical function during implementation. +The scale of the sigmoid function ($\alpha$ on the order of $10^3$) should be chosen so that for all practical parameters of the TGOV1 model, the sigmoid acts as a step function. ```math \begin{aligned} - \sigma(x) =\dfrac{1}{1+\exp(-\alpha x)}\approx \dfrac{1}{2}\left(\dfrac{\alpha x}{1+|\alpha x|}\right) + \sigma(x) =\dfrac{1}{1+\exp(-\alpha x)} \end{aligned} ``` diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp index c4412d5f..95bb6fda 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1.hpp @@ -121,20 +121,10 @@ namespace GridKit // Input States (which can be parameters) ScalarT pref_{0}; - // Scale of Sigmoid function (temporary local implementation) - const ScalarT mu_{4000.0}; - /// Component signal extension ComponentSignals signals_; - // Activation function (sigmoid approximation) - ScalarT sigmoid(ScalarT x); - - // Indicator of Valve limit states - ScalarT indicator_low(ScalarT x, ScalarT f); - ScalarT indicator_high(ScalarT x, ScalarT f); - ScalarT indicator(ScalarT x, ScalarT f); - + // Parameter initialization function void initializeParameters(const model_data_type& data); /* Local copies of external variables */ diff --git a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp index 50c6709f..f6add0df 100644 --- a/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp +++ b/GridKit/Model/PhasorDynamics/Governor/Tgov1/Tgov1Impl.hpp @@ -229,48 +229,6 @@ namespace GridKit return 0; } - /** - * @brief Scaled sigmoid activation function - * - * Temporary local implementation of smooth approximation - * of a piecewise differential equation. Ideally this is - * a more abstracted capability with GK. - * - * Algebraic approximation of transcendental sigmoid. - */ - template - ScalarT Tgov1::sigmoid(ScalarT x) - { - return ((HALF * mu_ * x) / (ONE + std::abs(mu_ * x))) + HALF; - } - - /** - * @brief Indicator function for lower valve limit violation. - */ - template - ScalarT Tgov1::indicator_low(ScalarT x, ScalarT f) - { - return (this->sigmoid(Pvmin_ - x)) * (this->sigmoid(-f)); - } - - /** - * @brief Indicator function for high valve limit violation. - */ - template - ScalarT Tgov1::indicator_high(ScalarT x, ScalarT f) - { - return (this->sigmoid(x - Pvmax_)) * (this->sigmoid(f)); - } - - /** - * @brief Net Indicator function for valve limits. - */ - template - ScalarT Tgov1::indicator(ScalarT x, ScalarT f) - { - return (ONE - this->indicator_low(x, f)) * (ONE - this->indicator_high(x, f)); - } - /** * @brief Internal residuals * @@ -297,7 +255,7 @@ namespace GridKit // The 'pre-limit' derivative of Pv ScalarT func = (-pv + (pref_ - omega) / R_) / T1_; - ScalarT valv_ind = this->indicator(pv, func); + ScalarT valv_ind = Math::indicator(Pvmin_, Pvmax_, pv, func); // Internal Differential Equations f[0] = -ptx_dot + pv - (ptx + T2_ * pv) / T3_;