diff --git a/CHANGELOG.md b/CHANGELOG.md index 066786e09..2eaacb7b3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -62,6 +62,7 @@ - Added cmake-format hooks, including in pre-commit. - Added off-nominal tap ratio and phase shift support to the PhasorDynamics `Branch` model. - Added portable Vector class to GridKit +- Added `HYGOV` governor model implementation for PhasorDynamics. ## v0.1 diff --git a/GridKit/CommonMath.hpp b/GridKit/CommonMath.hpp index 9bb3a5c3c..daaa11e94 100644 --- a/GridKit/CommonMath.hpp +++ b/GridKit/CommonMath.hpp @@ -153,7 +153,31 @@ namespace GridKit } /** - * @brief Smooth two-sided deadband function + * @brief Smooth Type 1 no-offset two-sided deadband function + * + * Smooth approximation to a deadband that returns zero inside the band and + * passes the input through unchanged outside the band. + * + * @tparam ScalarT - scalar data type + * @tparam RealT - Real data type (see GridKit::ScalarTraits::RealT) + * + * @param[in] x - Input signal + * @param[in] lower - Lower breakpoint + * @param[in] upper - Upper breakpoint + * @return Smooth no-offset deadbanded value + */ + template + __attribute__((always_inline)) inline ScalarT deadband1( + const ScalarT x, + const RealT lower, + const RealT upper) + { + assert(lower <= upper); + return x * (sigmoid(lower - x) + sigmoid(x - upper)); + } + + /** + * @brief Smooth Type 2 offset two-sided deadband function * * Smooth approximation to x - min(max(x, lower), upper), composed from the * smooth ramp function. @@ -164,10 +188,10 @@ namespace GridKit * @param[in] x - Input signal * @param[in] lower - Lower breakpoint * @param[in] upper - Upper breakpoint - * @return Smooth deadbanded value + * @return Smooth offset deadbanded value */ template - __attribute__((always_inline)) inline ScalarT deadband( + __attribute__((always_inline)) inline ScalarT deadband2( const ScalarT x, const RealT lower, const RealT upper) diff --git a/GridKit/CommonMath.md b/GridKit/CommonMath.md index a654487a2..29a9a7901 100644 --- a/GridKit/CommonMath.md +++ b/GridKit/CommonMath.md @@ -49,7 +49,8 @@ The scale $\mu=4\cdot f_{\text{sync}}=240$ is chosen so $\sigma$ behaves like a | `max` | Smooth binary maximum | `REECA` | | `min` | Smooth binary minimum | `REECA` | | `clamp` | Bounded saturation | `IEEEST` | -| `deadband` | Signed two-sided deadband | `REECA` | +| `deadband1` | Type 1 no-offset signed two-sided deadband | - | +| `deadband2` | Type 2 offset signed two-sided deadband | `REECA` | | `slew` | Symmetric slew-rate limiter | - | | `linseg` | Saturated linear segment contribution | `REGCA`, `REECA` | | `above` | Above-lower-limit indicator | - | @@ -84,7 +85,15 @@ x & \ell \le x \le u \\ u & x > u \end{cases} \\ -\text{deadband}(x;\ell,u) +\text{deadband1}(x;\ell,u) +&= +\begin{cases} +x & x < \ell \\ +0 & \ell \le x \le u \\ +x & x > u +\end{cases} +\\ +\text{deadband2}(x;\ell,u) &= \begin{cases} x-\ell & x < \ell \\ @@ -156,7 +165,8 @@ f & x \ge u \land f < 0 \\ \text{max}(x,y) &= y+\rho(x-y) \\ \text{min}(x,y) &= x-\rho(x-y) \\ \text{clamp}(x;\ell,u) &= \ell+\rho(x-\ell)-\rho(x-u) \\ -\text{deadband}(x;\ell,u) &= \rho(x-u)-\rho(\ell-x) \\ +\text{deadband1}(x;\ell,u) &= x\left[\sigma(\ell-x)+\sigma(x-u)\right] \\ +\text{deadband2}(x;\ell,u) &= \rho(x-u)-\rho(\ell-x) \\ \text{slew}(f;r) &= -r+\rho(f+r)-\rho(f-r) \\ \text{linseg}(x;a,b,h) &= \dfrac{h}{b-a}\left[\rho(x-a)-\rho(x-b)\right] \\ \text{above}(x;\ell) &= \sigma(x-\ell) \\ diff --git a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp index e48a8a9c5..0a033850a 100644 --- a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp +++ b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp @@ -6,6 +6,7 @@ #include #include #include +#include #include #include #include diff --git a/GridKit/Model/PhasorDynamics/Governor/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Governor/CMakeLists.txt index 7c7269784..fbe71c740 100644 --- a/GridKit/Model/PhasorDynamics/Governor/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/Governor/CMakeLists.txt @@ -4,3 +4,4 @@ # ]] add_subdirectory(Tgov1) +add_subdirectory(HYGOV) diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Governor/HYGOV/CMakeLists.txt new file mode 100644 index 000000000..1719101b0 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/CMakeLists.txt @@ -0,0 +1,54 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +set(_install_headers Hygov.hpp HygovData.hpp) + +if(GRIDKIT_ENABLE_ENZYME) + gridkit_add_library( + phasor_dynamics_governor_hygov + SOURCES HygovEnzyme.cpp + HEADERS ${_install_headers} + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal + PRIVATE + ClangEnzymeFlags + COMPILE_OPTIONS + PRIVATE + -mllvm + -enzyme-auto-sparsity=1 + -fno-math-errno) +else() + gridkit_add_library( + phasor_dynamics_governor_hygov + SOURCES Hygov.cpp + HEADERS ${_install_headers} + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal) +endif() + +gridkit_add_library( + phasor_dynamics_governor_hygov_dependency_tracking + SOURCES HygovDependencyTracking.cpp + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES + PUBLIC + GridKit::phasor_dynamics_core + PUBLIC + GridKit::phasor_dynamics_signal_dependency_tracking) + +target_link_libraries( + phasor_dynamics_components + INTERFACE GridKit::phasor_dynamics_governor_hygov) +target_link_libraries( + phasor_dynamics_components_dependency_tracking + INTERFACE GridKit::phasor_dynamics_governor_hygov_dependency_tracking) diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.cpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.cpp new file mode 100644 index 000000000..bdb7e15c0 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.cpp @@ -0,0 +1,27 @@ +/** + * @file Hygov.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Non-Enzyme instantiation for the HYGOV governor model. + */ + +#include "HygovImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + template + int Hygov::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Hygov..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Hygov; + template class Hygov; + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp new file mode 100644 index 000000000..19adabede --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp @@ -0,0 +1,154 @@ +/** + * @file Hygov.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Declaration of the HYGOV governor model. + */ + +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + template + class SignalNode; + + namespace Governor + { + /// Internal variables of a `Hygov`. + enum class HygovInternalVariables : size_t + { + XN, ///< Speed lead-lag denominator state + XF, ///< Governor error filter output + C, ///< Desired-gate position + G, ///< Gate position + Q, ///< Turbine flow + OMEGADB, ///< Deadbanded speed deviation + YOMEGA, ///< Lead-lag-conditioned speed signal + EF, ///< Governor error into the filter + FC, ///< Desired-gate derivative target + RC, ///< Rate-limited desired-gate derivative target + PGV, ///< Nonlinear gate-to-power curve output + PGVSAFE, ///< Safe gate-to-power value + H, ///< Turbine head + PMECH, ///< Mechanical-power output + MAXIMUM, + }; + + /// External variables of a `Hygov`. + enum class HygovExternalVariables : size_t + { + OMEGA, ///< Machine speed deviation + PREF, ///< Active-power/load reference + PAUX, ///< Auxiliary power input + MAXIMUM, + }; + + template + class Hygov : public Component + { + using Component::alpha_; + using Component::f_; + using Component::gridkit_component_id_; + using Component::J_; + using Component::J_cols_buffer_; + using Component::J_rows_buffer_; + using Component::J_vals_buffer_; + using Component::residual_indices_; + using Component::size_; + using Component::tag_; + using Component::variable_indices_; + using Component::wb_; + using Component::y_; + using Component::yp_; + + public: + using ScalarT = scalar_type; + using IdxT = index_type; + using RealT = typename Component::RealT; + using SignalT = SignalNode; + using ModelDataT = HygovData; + using MonitorT = Model::VariableMonitor; + + Hygov(); + Hygov(const ModelDataT& data); + ~Hygov(); + + int setGridKitComponentID(IdxT) override final; + int allocate() override final; + int verify() const override final; + int initialize() override final; + int tagDifferentiable() override final; + int evaluateResidual() override final; + int evaluateJacobian() override final; + + auto getSignals() + -> ComponentSignals& + { + return signals_; + } + + const Model::VariableMonitorBase* getMonitor() const override; + + __attribute__((always_inline)) inline int evaluateInternalResidual( + ScalarT*, ScalarT*, ScalarT*, ScalarT*, ScalarT*); + + private: + void initModelParams(const ModelDataT& data); + void initializeMonitor(); + + ScalarT gatePower(ScalarT gate) const; + RealT invertGatePower(RealT pgv) const; + ScalarT toComponentBase(ScalarT value) const; + ScalarT toPmechBase(ScalarT value) const; + + RealT Trate_{0}; + RealT pmech_mva_base_{0}; + RealT Rperm_{0}; + RealT Rtemp_{0}; + RealT Tr_{0}; + RealT Tf_{0}; + RealT Tg_{0}; + RealT Velm_{0}; + RealT Gmax_{0}; + RealT Gmin_{0}; + RealT Tw_{0}; + RealT At_{0}; + RealT Dturb_{0}; + RealT Qnl_{0}; + RealT Tn_{0}; + RealT Tnp_{0}; + RealT leadlag_gain_{0}; + RealT db1_{0}; + RealT db2_{0}; + RealT Hdam_{1}; + std::array Gv_{}; + std::array Pgv_{}; + + IdxT parameter_error_count_{0}; + + ScalarT pref_set_{0}; + ScalarT paux_set_{0}; + + ComponentSignals signals_; + std::unique_ptr monitor_; + + std::vector ws_; + std::vector ws_indices_; + }; + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovData.hpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovData.hpp new file mode 100644 index 000000000..9e2ebfda9 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovData.hpp @@ -0,0 +1,91 @@ +/** + * @file HygovData.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Modeling data for the HYGOV governor model. + */ + +#pragma once + +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + /// Parameter keys for the HYGOV governor model. + enum class HygovParameters + { + Trate, ///< Turbine-rating power base + mva, ///< MVA base of the pmech signal + Rperm, ///< Permanent droop + Rtemp, ///< Temporary droop + Tr, ///< Temporary-droop reset time constant + Tf, ///< Governor error filter time constant + Tg, ///< Gate servo time constant + Velm, ///< Maximum desired-gate velocity magnitude + Gmax, ///< Maximum desired-gate position + Gmin, ///< Minimum desired-gate position + Tw, ///< Water inertia time constant + At, ///< Turbine gain + Dturb, ///< Turbine damping coefficient + Qnl, ///< No-load flow at nominal head + Tn, ///< Speed lead-lag numerator time constant + Tnp, ///< Speed lead-lag denominator time constant + db1, ///< Type 1 speed deadband threshold + db2, ///< Unsupported mechanical backlash deadband + Hdam, ///< Head available at dam + Gv0, ///< Gate point 0 + Gv1, ///< Gate point 1 + Gv2, ///< Gate point 2 + Gv3, ///< Gate point 3 + Gv4, ///< Gate point 4 + Gv5, ///< Gate point 5 + Pgv0, ///< Power point 0 + Pgv1, ///< Power point 1 + Pgv2, ///< Power point 2 + Pgv3, ///< Power point 3 + Pgv4, ///< Power point 4 + Pgv5 ///< Power point 5 + }; + + /// Ports for the HYGOV governor model. + enum class HygovPorts + { + bus, ///< Optional terminal bus ID for case-format compatibility + speed, ///< Machine speed-deviation signal ID + pmech, ///< Mechanical-power output signal ID + pref, ///< Optional active-power/load reference signal ID + paux ///< Optional auxiliary power input signal ID + }; + + /// Variables available through the monitor interface. + enum class HygovMonitorableVariables + { + pmech, ///< Mechanical power output + leadlag, ///< Speed lead-lag denominator state + filter, ///< Governor error filter output + desiredgate, ///< Desired-gate position + gate, ///< Gate position + flow, ///< Turbine flow + head, ///< Turbine head + pgv ///< Nonlinear gate-to-power curve output + }; + + template + struct HygovData : public ComponentData + { + HygovData() = default; + + using Parameters = HygovParameters; + using Ports = HygovPorts; + using MonitorableVariables = HygovMonitorableVariables; + }; + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovDependencyTracking.cpp new file mode 100644 index 000000000..760bac957 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovDependencyTracking.cpp @@ -0,0 +1,27 @@ +/** + * @file HygovDependencyTracking.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Dependency-tracking instantiations for the HYGOV governor model. + */ + +#include "HygovImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + template + int Hygov::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Hygov..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class Hygov; + template class Hygov; + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovEnzyme.cpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovEnzyme.cpp new file mode 100644 index 000000000..7c5ff0aac --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovEnzyme.cpp @@ -0,0 +1,111 @@ +/** + * @file HygovEnzyme.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Enzyme sparse Jacobian for the HYGOV governor model. + */ + +#include + +#include "HygovImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + namespace + { + template + real_type deadband1Derivative(real_type x, real_type lower, real_type upper) + { + using RealT = real_type; + static constexpr RealT MU = static_cast(240.0); + + const RealT lower_sig = Math::sigmoid(lower - x); + const RealT upper_sig = Math::sigmoid(x - upper); + return lower_sig + upper_sig + + x * MU + * (upper_sig * (ONE - upper_sig) + - lower_sig * (ONE - lower_sig)); + } + } // namespace + + template + int Hygov::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Hygov..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + + J_.zeroMatrix(); + if (J_rows_buffer_ == nullptr) + { + J_rows_buffer_ = new IdxT[static_cast(size_) * static_cast(size_)]; + J_cols_buffer_ = new IdxT[static_cast(size_) * static_cast(size_)]; + J_vals_buffer_ = new RealT[static_cast(size_) * static_cast(size_)]; + } + + using ModelT = GridKit::PhasorDynamics::Governor::Hygov; + using Fn = GridKit::Enzyme::Sparse::MemberFunctions; + + GridKit::Enzyme::Sparse::DfDy::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + J_); + + const auto OMEGADB = static_cast(HygovInternalVariables::OMEGADB); + const auto EF = static_cast(HygovInternalVariables::EF); + const auto PMECH = static_cast(HygovInternalVariables::PMECH); + const auto G = static_cast(HygovInternalVariables::G); + + const auto OMEGA = static_cast(HygovExternalVariables::OMEGA); + const auto PREF = static_cast(HygovExternalVariables::PREF); + const auto PAUX = static_cast(HygovExternalVariables::PAUX); + + IdxT nnz = 0; + auto storeSignalDerivative = [&](size_t row, size_t col, RealT value) + { + if (ws_indices_[col] == INVALID_INDEX) + { + return; + } + + J_rows_buffer_[nnz] = residual_indices_.at(row); + J_cols_buffer_[nnz] = ws_indices_[col]; + J_vals_buffer_[nnz] = value; + ++nnz; + }; + + storeSignalDerivative(OMEGADB, + OMEGA, + deadband1Derivative(static_cast(ws_[OMEGA]), + -db1_, + db1_)); + storeSignalDerivative(EF, PREF, ONE); + storeSignalDerivative(EF, PAUX, ONE); + storeSignalDerivative(PMECH, + OMEGA, + -Dturb_ * static_cast(y_[G]) * Trate_ / pmech_mva_base_); + + J_.setValues(1.0, J_rows_buffer_, J_cols_buffer_, J_vals_buffer_, nnz); + return 0; + } + + template class Hygov; + template class Hygov; + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp new file mode 100644 index 000000000..a24158e3d --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp @@ -0,0 +1,553 @@ +/** + * @file HygovImpl.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Definition of the HYGOV governor model. + */ + +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + using Log = ::GridKit::Utilities::Logger; + + template + Hygov::Hygov() + { + size_ = static_cast(HygovInternalVariables::MAXIMUM); + } + + template + Hygov::Hygov(const ModelDataT& data) + : monitor_(std::make_unique(data)) + { + initModelParams(data); + initializeMonitor(); + size_ = static_cast(HygovInternalVariables::MAXIMUM); + } + + template + Hygov::~Hygov() + { + } + + template + scalar_type Hygov::toComponentBase(scalar_type value) const + { + return value * pmech_mva_base_ / Trate_; + } + + template + scalar_type Hygov::toPmechBase(scalar_type value) const + { + return value * Trate_ / pmech_mva_base_; + } + + template + void Hygov::initModelParams(const ModelDataT& data) + { + using Params = typename ModelDataT::Parameters; + + parameter_error_count_ = 0; + + auto load_required_real = [&](auto key, RealT& target, const char* name) + { + if (!data.parameters.contains(key)) + { + Log::error() << "Hygov: missing required parameter '" << name << "'\n"; + ++parameter_error_count_; + return; + } + + const auto& value = data.parameters.at(key); + if (const auto* real_value = std::get_if(&value)) + { + target = *real_value; + } + else if (const auto* index_value = std::get_if(&value)) + { + target = static_cast(*index_value); + } + else + { + Log::error() << "Hygov: parameter '" << name << "' must be numeric\n"; + ++parameter_error_count_; + } + }; + + auto load_optional_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() << "Hygov: parameter '" << name << "' must be numeric\n"; + ++parameter_error_count_; + } + }; + + load_required_real(Params::Trate, Trate_, "Trate"); + load_required_real(Params::mva, pmech_mva_base_, "mva"); + load_required_real(Params::Rperm, Rperm_, "Rperm"); + load_required_real(Params::Rtemp, Rtemp_, "Rtemp"); + load_required_real(Params::Tr, Tr_, "Tr"); + load_required_real(Params::Tf, Tf_, "Tf"); + load_required_real(Params::Tg, Tg_, "Tg"); + load_required_real(Params::Velm, Velm_, "Velm"); + load_required_real(Params::Gmax, Gmax_, "Gmax"); + load_required_real(Params::Gmin, Gmin_, "Gmin"); + load_required_real(Params::Tw, Tw_, "Tw"); + load_required_real(Params::At, At_, "At"); + load_required_real(Params::Dturb, Dturb_, "Dturb"); + load_required_real(Params::Qnl, Qnl_, "Qnl"); + load_required_real(Params::Tn, Tn_, "Tn"); + load_required_real(Params::Tnp, Tnp_, "Tnp"); + load_required_real(Params::db1, db1_, "db1"); + load_optional_real(Params::db2, db2_, "db2"); + load_required_real(Params::Hdam, Hdam_, "Hdam"); + + load_required_real(Params::Gv0, Gv_[0], "Gv0"); + load_required_real(Params::Gv1, Gv_[1], "Gv1"); + load_required_real(Params::Gv2, Gv_[2], "Gv2"); + load_required_real(Params::Gv3, Gv_[3], "Gv3"); + load_required_real(Params::Gv4, Gv_[4], "Gv4"); + load_required_real(Params::Gv5, Gv_[5], "Gv5"); + load_required_real(Params::Pgv0, Pgv_[0], "Pgv0"); + load_required_real(Params::Pgv1, Pgv_[1], "Pgv1"); + load_required_real(Params::Pgv2, Pgv_[2], "Pgv2"); + load_required_real(Params::Pgv3, Pgv_[3], "Pgv3"); + load_required_real(Params::Pgv4, Pgv_[4], "Pgv4"); + load_required_real(Params::Pgv5, Pgv_[5], "Pgv5"); + + const bool source_default_curve = + std::all_of(Gv_.begin(), Gv_.end(), [](RealT value) + { return value == ZERO; }) + && std::all_of(Pgv_.begin(), Pgv_.end(), [](RealT value) + { return value == ZERO; }); + if (source_default_curve) + { + Gv_ = {ZERO, + static_cast(0.2), + static_cast(0.4), + static_cast(0.6), + static_cast(0.8), + ONE}; + Pgv_ = Gv_; + } + + leadlag_gain_ = ZERO; + if (Tnp_ > ZERO) + { + leadlag_gain_ = Tn_ / Tnp_; + } + } + + template + const Model::VariableMonitorBase* Hygov::getMonitor() const + { + return monitor_.get(); + } + + template + void Hygov::initializeMonitor() + { + using Variable = typename ModelDataT::MonitorableVariables; + auto index = [](HygovInternalVariables variable) + { + return static_cast(variable); + }; + + monitor_->set(Variable::pmech, [this, index] + { return y_[index(HygovInternalVariables::PMECH)]; }); + monitor_->set(Variable::leadlag, [this, index] + { return y_[index(HygovInternalVariables::XN)]; }); + monitor_->set(Variable::filter, [this, index] + { return y_[index(HygovInternalVariables::XF)]; }); + monitor_->set(Variable::desiredgate, [this, index] + { return y_[index(HygovInternalVariables::C)]; }); + monitor_->set(Variable::gate, [this, index] + { return y_[index(HygovInternalVariables::G)]; }); + monitor_->set(Variable::flow, [this, index] + { return y_[index(HygovInternalVariables::Q)]; }); + monitor_->set(Variable::head, [this, index] + { return y_[index(HygovInternalVariables::H)]; }); + monitor_->set(Variable::pgv, [this, index] + { return y_[index(HygovInternalVariables::PGV)]; }); + } + + template + int Hygov::setGridKitComponentID(IdxT component_id) + { + gridkit_component_id_ = component_id; + return 0; + } + + template + int Hygov::allocate() + { + size_ = static_cast(HygovInternalVariables::MAXIMUM); + auto size = static_cast(size_); + + f_.assign(size, ScalarT{0}); + y_.assign(size, ScalarT{0}); + yp_.assign(size, ScalarT{0}); + tag_.assign(size, false); + variable_indices_.resize(size); + residual_indices_.resize(size); + + wb_.clear(); + + auto signal_size = static_cast(HygovExternalVariables::MAXIMUM); + ws_.assign(signal_size, ScalarT{0}); + ws_indices_.assign(signal_size, INVALID_INDEX); + + for (IdxT j = 0; j < size_; ++j) + { + this->setVariableIndex(j, j); + this->setResidualIndex(j, j); + } + + if (signals_.template isAssigned()) + { + signals_.template getSignalNode()->set( + &y_[static_cast(HygovInternalVariables::PMECH)], + &(this->getVariableIndex(static_cast(HygovInternalVariables::PMECH)))); + } + + return 0; + } + + template + int Hygov::verify() const + { + int ret = static_cast(parameter_error_count_); + + auto check = [&](bool condition, const char* message) + { + if (!condition) + { + Log::error() << "Hygov: " << message << '\n'; + ret += 1; + } + }; + + check(Trate_ > ZERO, "Trate must be positive"); + check(pmech_mva_base_ > ZERO, "mva must be positive"); + check(Rperm_ > ZERO, "Rperm must be positive"); + check(Rtemp_ > ZERO, "Rtemp must be positive"); + check(Tr_ > ZERO, "Tr must be positive"); + check(Tf_ > ZERO, "Tf must be positive"); + check(Tg_ >= ZERO, "Tg must be non-negative"); + check(Tw_ > ZERO, "Tw must be positive"); + check(Tn_ >= ZERO, "Tn must be non-negative"); + check(Tnp_ >= ZERO, "Tnp must be non-negative"); + check(Tnp_ > ZERO || Tn_ == ZERO, "Tn must be zero when Tnp is zero"); + check(Velm_ > ZERO, "Velm must be positive"); + check(Gmin_ <= Gmax_, "Gmin must be less than or equal to Gmax"); + check(At_ > ZERO, "At must be positive"); + check(Dturb_ >= ZERO, "Dturb must be non-negative"); + check(Qnl_ >= ZERO, "Qnl must be non-negative"); + check(db1_ >= ZERO, "db1 must be non-negative"); + check(db2_ == ZERO, "nonzero db2 is not supported"); + check(Hdam_ > ZERO, "Hdam must be positive"); + + for (size_t i = 1; i < Gv_.size(); ++i) + { + check(Gv_[i - 1] < Gv_[i], "Gv points must be strictly increasing"); + check(Pgv_[i - 1] <= Pgv_[i], "Pgv points must be non-decreasing"); + } + check(Pgv_[0] >= ZERO, "Pgv0 must be non-negative"); + + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "Hygov: omega signal attached with no linked source\n"; + ret += 1; + } + + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "Hygov: pref signal attached with no linked source\n"; + ret += 1; + } + + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "Hygov: paux signal attached with no linked source\n"; + ret += 1; + } + + return ret; + } + + template + scalar_type Hygov::gatePower(scalar_type gate) const + { + return ScalarT{Pgv_[0]} + + Math::linseg(gate, Gv_[0], Gv_[1], Pgv_[1] - Pgv_[0]) + + Math::linseg(gate, Gv_[1], Gv_[2], Pgv_[2] - Pgv_[1]) + + Math::linseg(gate, Gv_[2], Gv_[3], Pgv_[3] - Pgv_[2]) + + Math::linseg(gate, Gv_[3], Gv_[4], Pgv_[4] - Pgv_[3]) + + Math::linseg(gate, Gv_[4], Gv_[5], Pgv_[5] - Pgv_[4]); + } + + template + typename Hygov::RealT + Hygov::invertGatePower( + typename Hygov::RealT pgv) const + { + static constexpr RealT tol = static_cast(1.0e-10); + + if (std::abs(pgv - Pgv_[0]) <= tol) + { + return Gv_[0]; + } + + for (size_t i = 0; i < 5; ++i) + { + if (Pgv_[i + 1] <= Pgv_[i]) + { + continue; + } + + if (Pgv_[i] - tol <= pgv && pgv <= Pgv_[i + 1] + tol) + { + const RealT fraction = (pgv - Pgv_[i]) / (Pgv_[i + 1] - Pgv_[i]); + return Gv_[i] + fraction * (Gv_[i + 1] - Gv_[i]); + } + } + + return std::numeric_limits::quiet_NaN(); + } + + template + int Hygov::initialize() + { + if (parameter_error_count_ > 0 || verify() > 0) + { + Log::error() << "Hygov: cannot initialize with invalid configuration\n"; + return 1; + } + + const auto XN = static_cast(HygovInternalVariables::XN); + const auto XF = static_cast(HygovInternalVariables::XF); + const auto C = static_cast(HygovInternalVariables::C); + const auto G = static_cast(HygovInternalVariables::G); + const auto Q = static_cast(HygovInternalVariables::Q); + const auto OMEGADB = static_cast(HygovInternalVariables::OMEGADB); + const auto YOMEGA = static_cast(HygovInternalVariables::YOMEGA); + const auto EF = static_cast(HygovInternalVariables::EF); + const auto FC = static_cast(HygovInternalVariables::FC); + const auto RC = static_cast(HygovInternalVariables::RC); + const auto PGV = static_cast(HygovInternalVariables::PGV); + const auto PGVSAFE = static_cast(HygovInternalVariables::PGVSAFE); + const auto H = static_cast(HygovInternalVariables::H); + const auto PMECH = static_cast(HygovInternalVariables::PMECH); + + ScalarT omega0{ZERO}; + if (signals_.template isAttached()) + { + omega0 = signals_.template readExternalVariable(); + } + + paux_set_ = ScalarT{ZERO}; + if (signals_.template isAttached()) + { + paux_set_ = signals_.template readExternalVariable(); + } + + const ScalarT pmech0 = toComponentBase(y_[PMECH]); + y_[H] = Hdam_; + y_[Q] = Qnl_ + pmech0 / (At_ * y_[H]); + y_[PGV] = y_[Q] / std::sqrt(y_[H]); + y_[PGVSAFE] = Math::max(y_[PGV], static_cast(0.01)); + + const RealT gate0 = invertGatePower(static_cast(y_[PGV])); + if (std::isnan(gate0)) + { + Log::error() << "Hygov: initial Pgv is outside the invertible gate curve\n"; + return 1; + } + + y_[G] = gate0; + y_[C] = y_[G]; + + if (y_[C] < Gmin_ || y_[C] > Gmax_) + { + Log::error() << "Hygov: initialized gate is outside Gmin/Gmax\n"; + return 1; + } + + y_[OMEGADB] = Math::deadband1(omega0, -db1_, db1_); + y_[XN] = y_[OMEGADB]; + y_[YOMEGA] = y_[OMEGADB]; + y_[XF] = ZERO; + y_[EF] = ZERO; + y_[FC] = ZERO; + y_[RC] = ZERO; + y_[PMECH] = toPmechBase(pmech0); + + pref_set_ = y_[EF] - paux_set_ + y_[YOMEGA] + Rperm_ * y_[C]; + if (signals_.template isAttached()) + { + const ScalarT pref0 = + signals_.template readExternalVariable(); + const RealT pref_err = static_cast(pref0 - pref_set_); + if (std::abs(pref_err) > static_cast(1.0e-10)) + { + Log::error() << "Hygov: pref initial condition is not steady state\n"; + return 1; + } + pref_set_ = pref0; + } + + std::fill(yp_.begin(), yp_.end(), ZERO); + return 0; + } + + template + int Hygov::tagDifferentiable() + { + std::fill(tag_.begin(), tag_.end(), false); + tag_[static_cast(HygovInternalVariables::XN)] = (Tnp_ > ZERO); + tag_[static_cast(HygovInternalVariables::XF)] = true; + tag_[static_cast(HygovInternalVariables::C)] = true; + tag_[static_cast(HygovInternalVariables::G)] = (Tg_ > ZERO); + tag_[static_cast(HygovInternalVariables::Q)] = true; + return 0; + } + + template + __attribute__((always_inline)) inline int + Hygov::evaluateInternalResidual( + ScalarT* y, + ScalarT* yp, + [[maybe_unused]] ScalarT* wb, + ScalarT* ws, + ScalarT* f) + { + const auto XN = static_cast(HygovInternalVariables::XN); + const auto XF = static_cast(HygovInternalVariables::XF); + const auto C = static_cast(HygovInternalVariables::C); + const auto G = static_cast(HygovInternalVariables::G); + const auto Q = static_cast(HygovInternalVariables::Q); + const auto OMEGADB = static_cast(HygovInternalVariables::OMEGADB); + const auto YOMEGA = static_cast(HygovInternalVariables::YOMEGA); + const auto EF = static_cast(HygovInternalVariables::EF); + const auto FC = static_cast(HygovInternalVariables::FC); + const auto RC = static_cast(HygovInternalVariables::RC); + const auto PGV = static_cast(HygovInternalVariables::PGV); + const auto PGVSAFE = static_cast(HygovInternalVariables::PGVSAFE); + const auto H = static_cast(HygovInternalVariables::H); + const auto PMECH = static_cast(HygovInternalVariables::PMECH); + + const auto OMEGA = static_cast(HygovExternalVariables::OMEGA); + const auto PREF = static_cast(HygovExternalVariables::PREF); + const auto PAUX = static_cast(HygovExternalVariables::PAUX); + + const ScalarT xn = y[XN]; + const ScalarT xf = y[XF]; + const ScalarT c = y[C]; + const ScalarT g = y[G]; + const ScalarT q = y[Q]; + const ScalarT omegadb = y[OMEGADB]; + const ScalarT yomega = y[YOMEGA]; + const ScalarT ef = y[EF]; + const ScalarT fc = y[FC]; + const ScalarT rc = y[RC]; + const ScalarT pgv = y[PGV]; + const ScalarT pgv_safe = y[PGVSAFE]; + const ScalarT head = y[H]; + const ScalarT pmech = y[PMECH]; + + const ScalarT omega = ws[OMEGA]; + const ScalarT pref = ws[PREF]; + const ScalarT paux = ws[PAUX]; + + f[XN] = -Tnp_ * yp[XN] - xn + omegadb; + f[XF] = -Tf_ * yp[XF] - xf + ef; + f[FC] = -Rtemp_ * Tr_ * fc + xf + Tr_ * yp[XF]; + f[C] = -yp[C] + Math::antiwindup(c, rc, Gmin_, Gmax_); + f[G] = -Tg_ * yp[G] - g + c; + f[Q] = -Tw_ * yp[Q] + Hdam_ - head; + + f[OMEGADB] = -omegadb + Math::deadband1(omega, -db1_, db1_); + f[YOMEGA] = -yomega + xn + leadlag_gain_ * (omegadb - xn); + f[EF] = -ef + pref + paux - yomega - Rperm_ * c; + f[RC] = -rc + Math::clamp(fc, -Velm_, Velm_); + f[PGV] = -pgv + gatePower(g); + f[PGVSAFE] = -pgv_safe + Math::max(pgv, static_cast(0.01)); + f[H] = -head + (q / pgv_safe) * (q / pgv_safe); + const ScalarT pmech_output = At_ * head * (q - Qnl_) - Dturb_ * omega * g; + f[PMECH] = -pmech + toPmechBase(pmech_output); + + return 0; + } + + template + int Hygov::evaluateResidual() + { + const auto OMEGA = static_cast(HygovExternalVariables::OMEGA); + const auto PREF = static_cast(HygovExternalVariables::PREF); + const auto PAUX = static_cast(HygovExternalVariables::PAUX); + + ws_[OMEGA] = ZERO; + ws_[PREF] = pref_set_; + ws_[PAUX] = paux_set_; + std::fill(ws_indices_.begin(), ws_indices_.end(), INVALID_INDEX); + + if (signals_.template isAttached()) + { + ws_[OMEGA] = signals_.template readExternalVariable(); + ws_indices_[OMEGA] = + signals_.template readExternalVariableIndex(); + } + + if (signals_.template isAttached()) + { + ws_[PREF] = signals_.template readExternalVariable(); + ws_indices_[PREF] = + signals_.template readExternalVariableIndex(); + } + + if (signals_.template isAttached()) + { + ws_[PAUX] = signals_.template readExternalVariable(); + ws_indices_[PAUX] = + signals_.template readExternalVariableIndex(); + } + + evaluateInternalResidual(y_.data(), yp_.data(), wb_.data(), ws_.data(), f_.data()); + return 0; + } + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/HYGOV/README.md b/GridKit/Model/PhasorDynamics/Governor/HYGOV/README.md new file mode 100644 index 000000000..d0957ecdd --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/README.md @@ -0,0 +1,291 @@ +# **Hydro Turbine-Governor Model (HYGOV)** + +HYGOV is a hydro turbine-governor model with a temporary-droop governor, gate +servo, and nonlinear single-penstock water-column turbine. In GridKit it is +represented as a governor model that reads machine speed deviation and supplies +mechanical power to the machine. + +Notes: +- HYGOV equations use the `Trate` base; the `pmech` signal uses `mva` and is converted at the signal boundary. +- GridKit requires `Trate > 0`; for PowerWorld `Trate = 0` records, set `Trate` to the connected machine `mva` before running. +- Aside from `db2`, model parameters are required input values; the table values are typical values, not parser defaults. +- HYGOV uses `db1` with CommonMath `deadband1`; HYGOVD-only `dbL`/`dbH` and nonzero `db2` are not modeled. +- PowerWorld fields `Ttur`, `Eps`, Kaplan blade-servo points, `Bmax`, and `Tblade` are not part of this implementation. +- Permanent-droop feedback uses desired gate $c$ (State 2), not gate $g$. + +## Block Diagram + +Standard HYGOV block diagram. + +
+ + + Figure 1: Governor HYGOV model. Figure courtesy of [PowerWorld](https://www.powerworld.com/WebHelp/) +
+ +## Model Parameters + +Symbol | Units | JSON | Description | Typical Value | Note +--------------------------------|----------|----------|----------------------------------------------|---------------|------ +$P^{\mathrm{rate}}$ | [MW] | `Trate` | Component-base turbine rating | 100.0 | Required positive value +$S^{\mathrm{pmech}}$ | [MVA] | `mva` | Mechanical-power signal base | 100.0 | Use connected machine `mva` +$R_{\mathrm{perm}}$ | [p.u.] | `Rperm` | Permanent droop | 0.04 | Source diagram label: `R` +$R_{\mathrm{temp}}$ | [p.u.] | `Rtemp` | Temporary droop | 0.3 | Source diagram label: `r` +$T_r$ | [sec] | `Tr` | Temporary-droop reset time constant | 5.0 | +$T_f$ | [sec] | `Tf` | Governor error filter time constant | 0.05 | State 1 +$T_g$ | [sec] | `Tg` | Gate servo time constant | 0.5 | State 3; if zero, $g$ is algebraic +$V_{\mathrm{elm}}$ | [p.u./s] | `Velm` | Maximum desired-gate velocity magnitude | 0.2 | Symmetric rate limit on State 2 +$G^{\max}$ | [p.u.] | `Gmax` | Maximum desired-gate position | 1.0 | +$G^{\min}$ | [p.u.] | `Gmin` | Minimum desired-gate position | 0.0 | +$T_w$ | [sec] | `Tw` | Water inertia time constant | 1.0 | State 4 +$A_t$ | [p.u.] | `At` | Turbine gain | 1.2 | +$D_{\mathrm{turb}}$ | [p.u.] | `Dturb` | Turbine damping coefficient | 0.5 | Multiplied by speed deviation and gate +$q_{\mathrm{NL}}$ | [p.u.] | `Qnl` | No-load flow at nominal head | 0.05 | +$T_n$ | [sec] | `Tn` | Speed lead-lag numerator time constant | 0.0 | Source text labels this "lag" +$T_{np}$ | [sec] | `Tnp` | Speed lead-lag denominator time constant | 0.0 | Source text labels this "lead" +$D_{\omega}$ | [p.u.] | `db1` | Type 1 speed deadband threshold | 0.0 | Source data may provide this in Hz +$H_{\mathrm{dam}}$ | [p.u.] | `Hdam` | Head available at dam | 1.0 | + +The nonlinear gate-to-power curve $N_{\mathrm{GV}}$ is represented by six source points: + +Symbol | Units | JSON | Description | Typical Value | Note +--------------------------------|--------|----------|----------------------------------------------|---------------|------ +$G_V^{(k)}$ | [p.u.] | `Gv0`-`Gv5` | Gate point $k$ of the gain curve | 0.0 | $k=0,\ldots,5$ +$P_{\mathrm{GV}}^{(k)}$ | [p.u.] | `Pgv0`-`Pgv5` | Power point $k$ of the gain curve | 0.0 | $k=0,\ldots,5$ + +If all six `Gv` and all six `Pgv` source values are zero, the source data is treated as +omitting the nonlinear curve. The effective curve used by the equations is the identity +curve with points `(0,0), (0.2,0.2), ..., (1,1)`. + +### Parameter Validation + +Invalid HYGOV parameter sets are rejected by the following checks. If source data preprocessing swaps gate limits, expands gate limits to include the initialized gate, adjusts small positive time constants, converts `db1` from Hz to p.u. speed, or replaces an all-zero source gate curve with the identity curve, apply these checks to the effective values used by the equations. Source data with nonzero `db2` is outside the supported equations below. + +The required checks are: + +```math +\begin{aligned} + &P^{\mathrm{rate}} > 0 \\ + &S^{\mathrm{pmech}} > 0 \\ + &R_{\mathrm{perm}} > 0 \\ + &R_{\mathrm{temp}} > 0 \\ + &T_r > 0 \\ + &T_f > 0,\quad T_g \ge 0 \\ + &T_w > 0 \\ + &T_n \ge 0,\quad T_{np} \ge 0,\quad T_{np}=0 \Rightarrow T_n=0 \\ + &V_{\mathrm{elm}} > 0 \\ + &G^{\min} \le G^{\max} \\ + &A_t > 0 \\ + &D_{\mathrm{turb}} \ge 0 \\ + &q_{\mathrm{NL}} \ge 0 \\ + &D_{\omega} \ge 0 \\ + &H_{\mathrm{dam}} > 0 \\ + &G_V^{(0)} < G_V^{(1)} < \cdots < G_V^{(5)} \\ + &0 \le P_{\mathrm{GV}}^{(0)} \le P_{\mathrm{GV}}^{(1)} \le \cdots \le P_{\mathrm{GV}}^{(5)} +\end{aligned} +``` + +Initialization also requires $N_{\mathrm{GV}}^{-1}$ to be single-valued at the solved initial operating point. + +### Model Derived Parameters + +```math +\begin{aligned} + k_n &= + \begin{cases} + T_n/T_{np} & T_{np} > 0 \\ + 0 & T_{np} = 0 + \end{cases} \\ + S_{\mathrm{hygov}}^{\mathrm{base}} + &= P^{\mathrm{rate}} +\end{aligned} +``` + +Here $S^{\mathrm{pmech}}$ is the `mva` parameter. The `pmech` signal uses +that base, while HYGOV turbine equations use the HYGOV component base: + +```math +\begin{aligned} + P^{\mathrm{hygov}} + &= P^{\mathrm{pmech}}\dfrac{S^{\mathrm{pmech}}}{S_{\mathrm{hygov}}^{\mathrm{base}}} \\ + P^{\mathrm{pmech}} + &= P^{\mathrm{hygov}}\dfrac{S_{\mathrm{hygov}}^{\mathrm{base}}}{S^{\mathrm{pmech}}} +\end{aligned} +``` + +The nonlinear gate-to-power function uses GridKit's smooth [Linear Segment](../../../../CommonMath.md#derived-functions) helper: + +```math +\begin{aligned} + N_{\mathrm{GV}}(x) + &= + P_{\mathrm{GV}}^{(0)} + + \sum_{k\in\{0,\ldots,4\}} + \text{linseg}\!\left( + x;\, + G_V^{(k)},\, + G_V^{(k+1)},\, + P_{\mathrm{GV}}^{(k+1)} - P_{\mathrm{GV}}^{(k)} + \right) +\end{aligned} +``` + +## Model Variables + +### Internal Variables + +#### Differential + +Symbol | Units | Description | Note +------------------------|--------|-------------------------------------|------ +$x_n$ | [p.u.] | Speed lead-lag denominator state | Not circled in Fig. 1; realizes the `Tn`/`Tnp` block +$x_f$ | [p.u.] | Governor error filter output | State 1 in Fig. 1 +$c$ | [p.u.] | Desired-gate position | State 2 in Fig. 1 +$g$ | [p.u.] | Gate position | State 3 in Fig. 1; algebraic when $T_g = 0$ +$q$ | [p.u.] | Turbine flow | State 4 in Fig. 1 + +#### Algebraic + +Symbol | Units | Description | Note +----------------------------------|----------|-------------------------------------|------ +$\omega_{\mathrm{db}}$ | [p.u.] | Type 1 deadbanded speed deviation | Defined by CommonMath `deadband1` +$y_{\omega}$ | [p.u.] | Lead-lag-conditioned speed signal | Output of the speed lead-lag +$e_f$ | [p.u.] | Governor error into the filter | Reference path less conditioned speed and permanent-droop feedback +$f_c$ | [p.u./s] | Desired-gate derivative target | Before rate and position limits +$r_c$ | [p.u./s] | Rate-limited desired-gate derivative target | Limited by $\pm V_{\mathrm{elm}}$ +$P_{\mathrm{GV}}$ | [p.u.] | Nonlinear gate-to-power curve output | $N_{\mathrm{GV}}(g)$ +$P_{\mathrm{GV}}^{\mathrm{safe}}$ | [p.u.] | Safe gate-to-power value | Lower bounded by 0.01 +$H$ | [p.u.] | Turbine head | $(q/P_{\mathrm{GV}}^{\mathrm{safe}})^2$ +$P_{\mathrm{mech}}$ | [p.u.] | Mechanical power to generator | `mva` signal; converted at the HYGOV boundary + +### External Variables + +#### Differential +None. + +#### Algebraic + +Symbol | Units | Description | Note +--------------------------------|--------|--------------------------------|------ +$\omega$ | [p.u.] | Machine speed deviation | Read from a machine model +$P_{\mathrm{ref}}$ | [p.u.] | Active-power/load reference | HYGOV component base; external setpoint or constant parameter; source label: `Pref` +$P_{\mathrm{aux}}$ | [p.u.] | Auxiliary power input | HYGOV component base; optional, defaults to zero; source label: `Paux` + +## Model Equations + +The desired-gate transfer block before limits is: + +```math +\begin{aligned} + \dfrac{c(s)}{x_f(s)} + &= \dfrac{1+sT_r}{R_{\mathrm{temp}}T_rs} +\end{aligned} +``` + +The realization below forms the derivative target $f_c$, limits it to $r_c$, and applies anti-windup at the desired-gate limits. + +### Differential Equations + +```math +\begin{aligned} + 0 &= -T_{np}\dot x_n - x_n + \omega_{\mathrm{db}} \\ + 0 &= -T_f\dot x_f - x_f + e_f \\ + 0 &= -R_{\mathrm{temp}}T_r f_c + x_f + T_r\dot x_f \\ + 0 &= + -\dot c + + \text{antiwindup}\!\left( + c, + r_c, + G^{\min}, + G^{\max} + \right) \\ + 0 &= -T_g\dot g - g + c \\ + 0 &= -T_w\dot q + H_{\mathrm{dam}} - H +\end{aligned} +``` + +CommonMath defines the [Anti-Windup](../../../../CommonMath.md#anti-windup-indicator) target and smooth approximation. + +### Algebraic Equations + +```math +\begin{aligned} + 0 &= -\omega_{\mathrm{db}} + + \text{deadband1}\!\left(\omega, -D_{\omega}, D_{\omega}\right) \\ + 0 &= -y_{\omega} + + x_n + + k_n\left(\omega_{\mathrm{db}} - x_n\right) \\ + 0 &= -e_f + P_{\mathrm{ref}} + P_{\mathrm{aux}} - y_{\omega} - R_{\mathrm{perm}}c \\ + 0 &= -r_c + \text{clamp}\!\left(f_c, -V_{\mathrm{elm}}, V_{\mathrm{elm}}\right) \\[0.5ex] + 0 &= -P_{\mathrm{GV}} + N_{\mathrm{GV}}(g) \\ + 0 &= -P_{\mathrm{GV}}^{\mathrm{safe}} + \text{max}\!\left(P_{\mathrm{GV}}, 0.01\right) \\ + 0 &= -H + \left(\dfrac{q}{P_{\mathrm{GV}}^{\mathrm{safe}}}\right)^2 \\ + 0 &= -P_{\mathrm{mech}}^{\mathrm{pmech}} + + \left(A_tH(q - q_{\mathrm{NL}}) - D_{\mathrm{turb}}\omega g\right) + \dfrac{S_{\mathrm{hygov}}^{\mathrm{base}}}{S^{\mathrm{pmech}}} +\end{aligned} +``` + +The $P_{\mathrm{GV}}^{\mathrm{safe}}$ lower bound protects the head divider near a closed gate. A physically consistent operating point satisfies $P_{\mathrm{GV}}>0.01$ so the lower bound is inactive. + +CommonMath defines the helper targets and smooth approximations for [max, clamp, deadband1, and linseg](../../../../CommonMath.md#derived-functions). + +## Initialization + +Initialization is performed by evaluating the steady-state residuals in dependency order. Let subscript $0$ denote initial values and set all internal derivatives to zero. If optional signals are not connected, use: + +```math +\begin{aligned} + \omega_0 &= 0 \\ + P_{\mathrm{aux},0} &= 0 +\end{aligned} +``` + +The connected power-source model initializes mechanical power first; HYGOV reads +the `pmech` signal $P_{\mathrm{mech},0}^{\mathrm{pmech}}$ and converts it to +the component-base value $P_{\mathrm{mech},0}^{\mathrm{hygov}}$. At synchronous +speed the damping term vanishes: + +```math +\begin{aligned} + P_{\mathrm{mech},0}^{\mathrm{hygov}} + &= P_{\mathrm{mech},0}^{\mathrm{pmech}} + \dfrac{S^{\mathrm{pmech}}}{S_{\mathrm{hygov}}^{\mathrm{base}}} \\ + H_0 &= H_{\mathrm{dam}} \\ + q_0 &= q_{\mathrm{NL}} + \dfrac{P_{\mathrm{mech},0}^{\mathrm{hygov}}}{A_tH_0} \\ + P_{\mathrm{GV},0} &= \dfrac{q_0}{\sqrt{H_0}} \\ + P_{\mathrm{GV},0}^{\mathrm{safe}} &= \text{max}\!\left(P_{\mathrm{GV},0}, 0.01\right) \\ + g_0 &= N_{\mathrm{GV}}^{-1}\!\left(P_{\mathrm{GV},0}\right) \\ + c_0 &= g_0 +\end{aligned} +``` + +Then evaluate the governor chain. With $\dot c_0=0$ and $\dot x_{f,0}=0$, the temporary-droop block requires $x_{f,0}=0$: + +```math +\begin{aligned} + \omega_{\mathrm{db},0} &= \text{deadband1}\!\left(\omega_0, -D_{\omega}, D_{\omega}\right) \\ + x_{n,0} &= \omega_{\mathrm{db},0} \\ + y_{\omega,0} &= \omega_{\mathrm{db},0} \\ + x_{f,0} &= 0 \\ + e_{f,0} &= x_{f,0} \\ + f_{c,0} &= 0 \\ + r_{c,0} &= 0 \\ + P_{\mathrm{ref},0} &= e_{f,0} - P_{\mathrm{aux},0} + y_{\omega,0} + R_{\mathrm{perm}}c_0 +\end{aligned} +``` + +For an unsaturated zero-derivative start, require $G^{\min} \le c_0 \le G^{\max}$, $P_{\mathrm{GV},0}>0.01$, and a single-valued inverse $N_{\mathrm{GV}}^{-1}(P_{\mathrm{GV},0})$. Starts that bind a gate limit or cannot be served at the available head are outside these closed-form initialization formulas. + +## Model Outputs + +Output | Units | Description | Note +---------------|--------|-------------------------------------|------ +`pmech` | [p.u.] | Mechanical power output | `mva` signal; oriented as mechanical input to the machine +`leadlag` | [p.u.] | Speed lead-lag denominator state | Not circled in Fig. 1 +`filter` | [p.u.] | Governor error filter output | State 1 +`desiredgate` | [p.u.] | Desired-gate position | State 2 +`gate` | [p.u.] | Gate position | State 3 +`flow` | [p.u.] | Turbine flow | State 4 +`head` | [p.u.] | Turbine head | $(q/P_{\mathrm{GV}}^{\mathrm{safe}})^2$ +`pgv` | [p.u.] | Nonlinear gate-to-power curve output | diff --git a/GridKit/Model/PhasorDynamics/Governor/README.md b/GridKit/Model/PhasorDynamics/Governor/README.md index 8e6305aa8..0a107bcf7 100644 --- a/GridKit/Model/PhasorDynamics/Governor/README.md +++ b/GridKit/Model/PhasorDynamics/Governor/README.md @@ -9,3 +9,4 @@ A governor models the control system that regulates the output power of a machin There are a few standard Governor models - Turbine Governor (See [TGOV1](Tgov1/README.md)) +- Hydro Turbine Governor (See [HYGOV](HYGOV/README.md)) diff --git a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md index 32df36cdc..e51950e9e 100644 --- a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md +++ b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md @@ -145,6 +145,7 @@ are specified: `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` `Tgov1 ` | the TGOV1 governor model | `pmech`, `speed` | `R`, `T1`, `T2`, `T3`, `Pvmax`, `Pvmin`, `Dt` | `none` + `Hygov` | the HYGOV hydro turbine-governor model | `pmech`, `speed`\*, `pref`\*, `paux`\* | `Trate`, `mva`, `Rperm`, `Rtemp`, `Tr`, `Tf`, `Tg`, `Velm`, `Gmax`, `Gmin`, `Tw`, `At`, `Dturb`, `Qnl`, `Tn`, `Tnp`, `db1`, `db2`, `Hdam`, `Gv0`-`Gv5`, `Pgv0`-`Pgv5` | `pmech`, `leadlag`, `filter`, `desiredgate`, `gate`, `flow`, `head`, `pgv` `Ieeet1` | the IEEET1 exciter model | `bus`, `speed`, `efd`, `vs`\* | `Tr`, `Ka`, `Ta`, `Ke`, `Te`, `Kf`, `Tf`, `Vrmin`, `Vrmax`, `E1`, `E2`, `Se1`, `Se2`, `Ispdlim` | `efd`, `ksat` `SexsPti` | the SEXS-PTI simplified exciter model | `bus`, `efd`, `vs`\* | `Ta`, `Tb`, `Te`, `K`, `Efdmax`, `Efdmin` | `efd` `Ieeest` | the IEEEST stabilizer model | `input`, `output` | `A1`, `A2`, `A3`, `A4`, `A5`, `A6`, `T1`, `T2`, `T3`, `T4`, `T5`, `T6`, `Ks`, `Lsmin`, `Lsmax`, `Vcl`, `Vcu`, `Tdelay` | `vss` diff --git a/GridKit/Model/PhasorDynamics/SystemModel.hpp b/GridKit/Model/PhasorDynamics/SystemModel.hpp index d7b4cc9d7..26761302a 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -297,6 +297,38 @@ namespace GridKit addComponent(gov); } + // Add HYGOV governors + for (const auto& govdata : data.hygov) + { + auto* gov = new Hygov(govdata); + + if (govdata.ports.contains(HygovData::Ports::speed)) + { + IdxT speed = govdata.ports.at(HygovData::Ports::speed); + gov->getSignals().template attachSignalNode(getSignal(speed)); + } + + if (govdata.ports.contains(HygovData::Ports::pmech)) + { + IdxT pmech = govdata.ports.at(HygovData::Ports::pmech); + gov->getSignals().template assignSignalNode(getSignal(pmech)); + } + + if (govdata.ports.contains(HygovData::Ports::pref)) + { + IdxT pref = govdata.ports.at(HygovData::Ports::pref); + gov->getSignals().template attachSignalNode(getSignal(pref)); + } + + if (govdata.ports.contains(HygovData::Ports::paux)) + { + IdxT paux = govdata.ports.at(HygovData::Ports::paux); + gov->getSignals().template attachSignalNode(getSignal(paux)); + } + + addComponent(gov); + } + for (const auto& excitedata : data.exciter) { IdxT bus_index = 0; diff --git a/GridKit/Model/PhasorDynamics/SystemModelData.hpp b/GridKit/Model/PhasorDynamics/SystemModelData.hpp index d4deef66e..017e512ca 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelData.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelData.hpp @@ -12,6 +12,7 @@ #include #include #include +#include #include #include #include @@ -41,6 +42,7 @@ namespace GridKit using BusToSignalAdapterDataT = BusToSignalAdapterData; using BusFaultDataT = BusFaultData; using Tgov1DataT = Governor::Tgov1Data; + using HygovDataT = Governor::HygovData; using Ieeet1DataT = Exciter::Ieeet1Data; using SexsPtiDataT = Exciter::SexsPtiData; using IeeestDataT = Stabilizer::IeeestData; @@ -99,6 +101,7 @@ namespace GridKit std::vector load; ///< Loads within the model std::vector loadzip; ///< Loads within the model std::vector gov; ///< Governors within the model + std::vector hygov; ///< HYGOV governors within the model std::vector exciter; ///< Exciters within the model std::vector sexspti; ///< SEXS-PTI exciters within the model std::vector stabilizer; ///< Stabilizers within the model diff --git a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp index 00a6278c2..198a75dc6 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp @@ -142,6 +142,12 @@ namespace GridKit raw_component.get_to(gov); sm.gov.push_back(gov); } + else if (kind == "Hygov") + { + typename SystemModelData::HygovDataT gov; + raw_component.get_to(gov); + sm.hygov.push_back(gov); + } else if (kind == "Ieeet1") { typename SystemModelData::Ieeet1DataT exciter; diff --git a/docs/Figures/PhasorDynamics/HYGOV_diagram.png b/docs/Figures/PhasorDynamics/HYGOV_diagram.png new file mode 100755 index 000000000..fc7d7916d Binary files /dev/null and b/docs/Figures/PhasorDynamics/HYGOV_diagram.png differ diff --git a/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp b/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp index 41f6d2ca1..0f7f659b3 100644 --- a/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp +++ b/tests/UnitTests/Math/SmoothnessIndicatorTests.hpp @@ -55,7 +55,31 @@ namespace GridKit return success.report(__func__); } - TestOutcome deadband() + TestOutcome deadband1() + { + TestStatus success = true; + + const ScalarT lower = scalar(-0.05); + const ScalarT upper = scalar(0.10); + const ScalarT tol = scalar(kSmoothTolerance); + + // Type 1 (no-offset) deadband returns ~0 inside the band ... + const ScalarT inside_input = scalar(0.02); + success *= (std::abs(Math::deadband1(inside_input, lower, upper)) < tol * tol); + + // ... and passes the input through unchanged (no offset) outside the band. + const ScalarT far_above_input = scalar(4.0); + const ScalarT far_below_input = scalar(-4.0); + + success *= std::isfinite(Math::deadband1(far_above_input, lower, upper)); + success *= within(Math::deadband1(far_above_input, lower, upper), far_above_input, tol); + success *= std::isfinite(Math::deadband1(far_below_input, lower, upper)); + success *= within(Math::deadband1(far_below_input, lower, upper), far_below_input, tol); + + return success.report(__func__); + } + + TestOutcome deadband2() { TestStatus success = true; @@ -69,12 +93,12 @@ namespace GridKit const ScalarT expected_below = below_input - lower; const ScalarT expected_above = above_input - upper; - success *= within(Math::deadband(below_input, lower, upper), expected_below, tol); - success *= (std::abs(Math::deadband(inside_input, lower, upper)) < tol * tol); - success *= within(Math::deadband(above_input, lower, upper), expected_above, tol); + success *= within(Math::deadband2(below_input, lower, upper), expected_below, tol); + success *= (std::abs(Math::deadband2(inside_input, lower, upper)) < tol * tol); + success *= within(Math::deadband2(above_input, lower, upper), expected_above, tol); - const ScalarT lower_breakpoint = Math::deadband(lower, lower, upper); - const ScalarT upper_breakpoint = Math::deadband(upper, lower, upper); + const ScalarT lower_breakpoint = Math::deadband2(lower, lower, upper); + const ScalarT upper_breakpoint = Math::deadband2(upper, lower, upper); success *= (lower_breakpoint < scalar(0.0)); success *= (upper_breakpoint > scalar(0.0)); @@ -82,7 +106,7 @@ namespace GridKit success *= (std::abs(upper_breakpoint) < tol); const ScalarT x = scalar(-0.4); - success *= (std::abs(Math::deadband(x, lower, upper) + success *= (std::abs(Math::deadband2(x, lower, upper) - (x - Math::clamp(x, lower, upper))) < scalar(kRoundoffTolerance)); @@ -91,17 +115,17 @@ namespace GridKit const ScalarT expected_far_above = far_above_input - upper; const ScalarT expected_far_below = far_below_input - lower; - success *= std::isfinite(Math::deadband(far_above_input, lower, upper)); - success *= within(Math::deadband(far_above_input, lower, upper), expected_far_above, tol); - success *= std::isfinite(Math::deadband(far_below_input, lower, upper)); - success *= within(Math::deadband(far_below_input, lower, upper), expected_far_below, tol); + success *= std::isfinite(Math::deadband2(far_above_input, lower, upper)); + success *= within(Math::deadband2(far_above_input, lower, upper), expected_far_above, tol); + success *= std::isfinite(Math::deadband2(far_below_input, lower, upper)); + success *= within(Math::deadband2(far_below_input, lower, upper), expected_far_below, tol); const ScalarT point = scalar(0.25); const ScalarT above_point = scalar(0.75); const ScalarT below_point = scalar(-0.25); - success *= (std::abs(Math::deadband(above_point, point, point) - (above_point - point)) + success *= (std::abs(Math::deadband2(above_point, point, point) - (above_point - point)) < scalar(kRoundoffTolerance)); - success *= (std::abs(Math::deadband(below_point, point, point) - (below_point - point)) + success *= (std::abs(Math::deadband2(below_point, point, point) - (below_point - point)) < scalar(kRoundoffTolerance)); return success.report(__func__); diff --git a/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp b/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp index 233c36ddd..aecdc7ae8 100644 --- a/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp +++ b/tests/UnitTests/Math/runSmoothnessIndicatorTests.cpp @@ -7,7 +7,8 @@ int main() GridKit::Testing::SmoothnessIndicatorTests test; result += test.clamp(); - result += test.deadband(); + result += test.deadband1(); + result += test.deadband2(); result += test.limitIndicators(); result += test.slew(); result += test.linseg(); diff --git a/tests/UnitTests/PhasorDynamics/CMakeLists.txt b/tests/UnitTests/PhasorDynamics/CMakeLists.txt index 932cb13b8..de1d7f639 100644 --- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt +++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt @@ -78,6 +78,14 @@ target_link_libraries( GridKit::phasor_dynamics_bus_dependency_tracking GridKit::testing) +add_executable(test_phasor_governor_hygov runGovernorHygovTests.cpp) +target_link_libraries( + test_phasor_governor_hygov + GridKit::definitions + GridKit::phasor_dynamics_components + GridKit::phasor_dynamics_components_dependency_tracking + GridKit::testing) + add_executable(test_phasor_exciter_ieeet1 runExciterIeeet1Tests.cpp) target_link_libraries( test_phasor_exciter_ieeet1 @@ -135,6 +143,7 @@ add_test(NAME PhasorDynamicsBusToSignalAdapterTest COMMAND test_phasor_bustosign add_test(NAME PhasorDynamicsBranchTest COMMAND test_phasor_branch) add_test(NAME PhasorDynamicsGenrouTest COMMAND test_phasor_genrou) add_test(NAME PhasorDynamicsGovernorTgov1Test COMMAND test_phasor_governor_tgov1) +add_test(NAME PhasorDynamicsGovernorHygovTest COMMAND test_phasor_governor_hygov) add_test(NAME PhasorDynamicsExciterIeeet1Test COMMAND test_phasor_exciter_ieeet1) add_test(NAME PhasorDynamicsGensalTest COMMAND test_phasor_gensal) add_test(NAME PhasorDynamicsExciterSexsPtiTest COMMAND test_phasor_exciter_sexspti) @@ -156,6 +165,7 @@ install( test_phasor_loadzip test_phasor_genrou test_phasor_governor_tgov1 + test_phasor_governor_hygov test_phasor_exciter_ieeet1 test_phasor_gensal test_phasor_exciter_sexspti diff --git a/tests/UnitTests/PhasorDynamics/GovernorHygovTests.hpp b/tests/UnitTests/PhasorDynamics/GovernorHygovTests.hpp new file mode 100644 index 000000000..d5495e229 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/GovernorHygovTests.hpp @@ -0,0 +1,363 @@ +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + template + class GovernorHygovTests + { + public: + using ScalarT = scalar_type; + using IdxT = index_type; + using RealT = typename PhasorDynamics::Component::RealT; + + static constexpr ScalarT kTol = static_cast(1.0e-8); + + TestOutcome constructionAndValidation() + { + TestStatus success = true; + + PhasorDynamics::Governor::Hygov hygov(makeHygovData()); + success *= (hygov.size() + == static_cast( + PhasorDynamics::Governor::HygovInternalVariables::MAXIMUM)); + success *= (hygov.getMonitor() != nullptr); + success *= (hygov.verify() == 0); + + auto source_default_hygov = makeHygovSourceDefaultData(); + PhasorDynamics::Governor::Hygov source_default_hygov_model( + source_default_hygov); + success *= (source_default_hygov_model.verify() == 0); + + auto bad_hygov = makeHygovData(); + bad_hygov.parameters[PhasorDynamics::Governor::HygovParameters::db2] = + static_cast(0.01); + PhasorDynamics::Governor::Hygov bad_hygov_model(bad_hygov); + success *= (bad_hygov_model.verify() > 0); + + auto bad_hygov_leadlag = makeHygovData(); + bad_hygov_leadlag.parameters[PhasorDynamics::Governor::HygovParameters::Tn] = + static_cast(0.1); + bad_hygov_leadlag.parameters[PhasorDynamics::Governor::HygovParameters::Tnp] = + static_cast(0.0); + PhasorDynamics::Governor::Hygov bad_hygov_leadlag_model( + bad_hygov_leadlag); + success *= (bad_hygov_leadlag_model.verify() > 0); + + auto bad_hygov_curve = makeHygovData(); + bad_hygov_curve.parameters[PhasorDynamics::Governor::HygovParameters::Gv2] = + static_cast(0.2); + PhasorDynamics::Governor::Hygov bad_hygov_curve_model(bad_hygov_curve); + success *= (bad_hygov_curve_model.verify() > 0); + + auto bad_hygov_trate = makeHygovData(); + bad_hygov_trate.parameters[PhasorDynamics::Governor::HygovParameters::Trate] = + static_cast(0.0); + PhasorDynamics::Governor::Hygov bad_hygov_trate_model(bad_hygov_trate); + success *= (bad_hygov_trate_model.verify() > 0); + + auto bad_hygov_base = makeHygovData(); + bad_hygov_base.parameters[PhasorDynamics::Governor::HygovParameters::mva] = + static_cast(0.0); + PhasorDynamics::Governor::Hygov bad_hygov_base_model(bad_hygov_base); + success *= (bad_hygov_base_model.verify() > 0); + + return success.report(__func__); + } + + TestOutcome signals() + { + using Var = PhasorDynamics::Governor::HygovInternalVariables; + + TestStatus success = true; + + PhasorDynamics::SignalNode pmech_node; + PhasorDynamics::SignalNode omega_node; + ScalarT pmech_value{0.40}; + ScalarT omega_value{0.0}; + IdxT pmech_index = 4; + IdxT omega_index = 5; + pmech_node.set(&pmech_value, &pmech_index); + omega_node.set(&omega_value, &omega_index); + + PhasorDynamics::Governor::Hygov hygov(makeHygovData()); + auto& hygov_signals = hygov.getSignals(); + hygov_signals + .template assignSignalNode( + &pmech_node); + hygov_signals + .template attachSignalNode( + &omega_node); + + success *= (hygov.allocate() == 0); + hygov.y()[index(Var::PMECH)] = pmech_value; + success *= (hygov.verify() == 0); + success *= (hygov.initialize() == 0); + success *= (hygov.tagDifferentiable() == 0); + success *= (hygov.evaluateResidual() == 0); + + success *= isEqual(hygov.y()[index(Var::Q)], static_cast(0.50), kTol); + success *= isEqual(hygov.y()[index(Var::G)], static_cast(0.50), kTol); + success *= isEqual(hygov.y()[index(Var::C)], static_cast(0.50), kTol); + success *= (hygov.tag()[index(Var::G)] == false); + + for (size_t i = 0; i < hygov.getResidual().size(); ++i) + { + success *= isEqual(hygov.getResidual()[i], static_cast(0.0), kTol); + success *= isEqual(hygov.yp()[i], static_cast(0.0), kTol); + } + + return success.report(__func__); + } + + TestOutcome sourceDefault() + { + using Var = PhasorDynamics::Governor::HygovInternalVariables; + + TestStatus success = true; + + PhasorDynamics::SignalNode pmech_node; + ScalarT pmech_value{0.40}; + IdxT pmech_index = 4; + pmech_node.set(&pmech_value, &pmech_index); + + PhasorDynamics::Governor::Hygov hygov(makeHygovSourceDefaultData()); + hygov.getSignals() + .template assignSignalNode( + &pmech_node); + + success *= (hygov.allocate() == 0); + hygov.y()[index(Var::PMECH)] = pmech_value; + success *= (hygov.verify() == 0); + success *= (hygov.initialize() == 0); + success *= (hygov.tagDifferentiable() == 0); + success *= (hygov.evaluateResidual() == 0); + + success *= isEqual(hygov.y()[index(Var::Q)], static_cast(0.50), kTol); + success *= isEqual(hygov.y()[index(Var::PGV)], static_cast(0.50), kTol); + success *= isEqual(hygov.y()[index(Var::G)], static_cast(0.50), kTol); + success *= (hygov.tag()[index(Var::XN)] == false); + + for (size_t i = 0; i < hygov.getResidual().size(); ++i) + { + success *= isEqual(hygov.getResidual()[i], static_cast(0.0), kTol); + success *= isEqual(hygov.yp()[i], static_cast(0.0), kTol); + } + + return success.report(__func__); + } + + TestOutcome baseConversion() + { + using Var = PhasorDynamics::Governor::HygovInternalVariables; + using Params = PhasorDynamics::Governor::HygovParameters; + + TestStatus success = true; + + PhasorDynamics::SignalNode pmech_node; + ScalarT pmech_value{0.40}; + IdxT pmech_index = 6; + pmech_node.set(&pmech_value, &pmech_index); + + auto data = makeHygovData(); + data.parameters[Params::Trate] = static_cast(50.0); + + PhasorDynamics::Governor::Hygov hygov(data); + hygov.getSignals().template assignSignalNode(&pmech_node); + + success *= (hygov.allocate() == 0); + hygov.y()[index(Var::PMECH)] = pmech_value; + success *= (hygov.verify() == 0); + success *= (hygov.initialize() == 0); + success *= (hygov.evaluateResidual() == 0); + + success *= isEqual(hygov.y()[index(Var::Q)], static_cast(0.90), kTol); + success *= isEqual(hygov.y()[index(Var::G)], static_cast(0.90), kTol); + success *= isEqual(hygov.y()[index(Var::PMECH)], static_cast(0.40), kTol); + success *= isEqual(pmech_node.read(), static_cast(0.40), kTol); + + for (size_t i = 0; i < hygov.getResidual().size(); ++i) + { + success *= isEqual(hygov.getResidual()[i], static_cast(0.0), kTol); + success *= isEqual(hygov.yp()[i], static_cast(0.0), kTol); + } + + return success.report(__func__); + } + + TestOutcome jsonParseAndSystemAssembly() + { + TestStatus success = true; + + std::istringstream input(R"json( +{ + "header": { + "format_version": 0, + "format_revision": 1, + "case_name": "hydro governor", + "case_description": "HYGOV parser test", + "case_comments": "", + "freq_base": 60.0, + "va_base": 100000000.0 + }, + "buses": [ + { + "number": 1, + "class": "bus", + "name": "Bus 1", + "init": { "Vr": 1.0, "Vi": 0.0 }, + "v_base": 1.0 + } + ], + "signals": [ + { "signal_id": 10, "name": "Pmech" } + ], + "devices": [ + { + "class": "Genrou", + "ports": { "bus": 1, "pmech": 10 }, + "id": "GEN1", + "params": { + "p0": 0.3, "q0": 0.0, "H": 3.0, "D": 0.0, "Ra": 0.0, + "Tdop": 7.0, "Tdopp": 0.04, "Tqop": 0.75, "Tqopp": 0.05, + "Xd": 2.1, "Xdp": 0.2, "Xdpp": 0.18, "Xq": 0.5, "Xqp": 0.5, + "Xqpp": 0.18, "Xl": 0.15, "S10": 0.0, "S12": 0.0, + "mva": 100.0 + } + }, + { + "class": "Hygov", + "ports": { "pmech": 10 }, + "id": "HYG1", + "params": { + "Trate": 50.0, "mva": 100.0, "Rperm": 0.05, "Rtemp": 0.4, "Tr": 5.0, + "Tf": 0.2, "Tg": 0.0, "Velm": 0.5, "Gmax": 1.0, "Gmin": 0.0, + "Tw": 1.0, "At": 1.0, "Dturb": 0.0, "Qnl": 0.1, + "Tn": 0.0, "Tnp": 1.0, "db1": 0.0, "db2": 0.0, "Hdam": 1.0, + "Gv0": 0.0, "Gv1": 0.2, "Gv2": 0.4, "Gv3": 0.6, "Gv4": 0.8, "Gv5": 1.0, + "Pgv0": 0.0, "Pgv1": 0.2, "Pgv2": 0.4, "Pgv3": 0.6, "Pgv4": 0.8, "Pgv5": 1.0 + } + } + ] +} +)json"); + + auto data = PhasorDynamics::parseSystemModelData(input); + success *= (data.hygov.size() == 1); + const auto trate_param = + data.hygov[0].parameters.at(PhasorDynamics::Governor::HygovParameters::Trate); + const auto mva_param = + data.hygov[0].parameters.at(PhasorDynamics::Governor::HygovParameters::mva); + success *= (std::get(trate_param) == static_cast(50.0)); + success *= (std::get(mva_param) == static_cast(100.0)); + PhasorDynamics::SystemModel system(data); + success *= (system.allocate() == 0); + success *= (system.initialize() == 0); + const auto hygov_size = + static_cast(PhasorDynamics::Governor::HygovInternalVariables::MAXIMUM); + const auto hygov_offset = static_cast(system.size() - hygov_size); + success *= isEqual( + system.y()[hygov_offset + index(PhasorDynamics::Governor::HygovInternalVariables::Q)], + static_cast(0.70), + kTol); + success *= isEqual( + system.y()[hygov_offset + index(PhasorDynamics::Governor::HygovInternalVariables::G)], + static_cast(0.70), + kTol); + success *= (system.evaluateResidual() == 0); + success *= (system.size() == 35); + + return success.report(__func__); + } + + private: + static size_t index(PhasorDynamics::Governor::HygovInternalVariables variable) + { + return static_cast(variable); + } + + auto makeHygovData() -> PhasorDynamics::Governor::HygovData + { + using Params = PhasorDynamics::Governor::HygovParameters; + using Mon = PhasorDynamics::Governor::HygovMonitorableVariables; + + PhasorDynamics::Governor::HygovData data; + data.device_class = "Hygov"; + data.disambiguation_string = "hygov_test"; + data.monitored_variables.insert(Mon::pmech); + data.monitored_variables.insert(Mon::gate); + + data.parameters[Params::Trate] = static_cast(100.0); + data.parameters[Params::mva] = static_cast(100.0); + data.parameters[Params::Rperm] = static_cast(0.05); + data.parameters[Params::Rtemp] = static_cast(0.4); + data.parameters[Params::Tr] = static_cast(5.0); + data.parameters[Params::Tf] = static_cast(0.2); + data.parameters[Params::Tg] = static_cast(0.0); + data.parameters[Params::Velm] = static_cast(0.5); + data.parameters[Params::Gmax] = static_cast(1.0); + data.parameters[Params::Gmin] = static_cast(0.0); + data.parameters[Params::Tw] = static_cast(1.0); + data.parameters[Params::At] = static_cast(1.0); + data.parameters[Params::Dturb] = static_cast(0.0); + data.parameters[Params::Qnl] = static_cast(0.1); + data.parameters[Params::Tn] = static_cast(0.0); + data.parameters[Params::Tnp] = static_cast(1.0); + data.parameters[Params::db1] = static_cast(0.0); + data.parameters[Params::db2] = static_cast(0.0); + data.parameters[Params::Hdam] = static_cast(1.0); + data.parameters[Params::Gv0] = static_cast(0.0); + data.parameters[Params::Gv1] = static_cast(0.2); + data.parameters[Params::Gv2] = static_cast(0.4); + data.parameters[Params::Gv3] = static_cast(0.6); + data.parameters[Params::Gv4] = static_cast(0.8); + data.parameters[Params::Gv5] = static_cast(1.0); + data.parameters[Params::Pgv0] = static_cast(0.0); + data.parameters[Params::Pgv1] = static_cast(0.2); + data.parameters[Params::Pgv2] = static_cast(0.4); + data.parameters[Params::Pgv3] = static_cast(0.6); + data.parameters[Params::Pgv4] = static_cast(0.8); + data.parameters[Params::Pgv5] = static_cast(1.0); + + return data; + } + + auto makeHygovSourceDefaultData() -> PhasorDynamics::Governor::HygovData + { + using Params = PhasorDynamics::Governor::HygovParameters; + + auto data = makeHygovData(); + data.parameters[Params::Tn] = static_cast(0.0); + data.parameters[Params::Tnp] = static_cast(0.0); + data.parameters[Params::Gv0] = static_cast(0.0); + data.parameters[Params::Gv1] = static_cast(0.0); + data.parameters[Params::Gv2] = static_cast(0.0); + data.parameters[Params::Gv3] = static_cast(0.0); + data.parameters[Params::Gv4] = static_cast(0.0); + data.parameters[Params::Gv5] = static_cast(0.0); + data.parameters[Params::Pgv0] = static_cast(0.0); + data.parameters[Params::Pgv1] = static_cast(0.0); + data.parameters[Params::Pgv2] = static_cast(0.0); + data.parameters[Params::Pgv3] = static_cast(0.0); + data.parameters[Params::Pgv4] = static_cast(0.0); + data.parameters[Params::Pgv5] = static_cast(0.0); + + return data; + } + }; + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/runGovernorHygovTests.cpp b/tests/UnitTests/PhasorDynamics/runGovernorHygovTests.cpp new file mode 100644 index 000000000..45154c841 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/runGovernorHygovTests.cpp @@ -0,0 +1,16 @@ +#include "GovernorHygovTests.hpp" + +int main() +{ + GridKit::Testing::TestingResults result; + + GridKit::Testing::GovernorHygovTests test; + + result += test.constructionAndValidation(); + result += test.signals(); + result += test.sourceDefault(); + result += test.baseConversion(); + result += test.jsonParseAndSystemAssembly(); + + return result.summary(); +}