diff --git a/CHANGELOG.md b/CHANGELOG.md index ba0f1dde6..6761bb82f 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 `HYGOV` governor model implementation for PhasorDynamics. ## v0.1 diff --git a/GridKit/CommonMath.hpp b/GridKit/CommonMath.hpp index e4a7ec87f..8f9f8c5cf 100644 --- a/GridKit/CommonMath.hpp +++ b/GridKit/CommonMath.hpp @@ -162,7 +162,31 @@ namespace GridKit } /** - * @brief Smooth two-sided deadband function + * @brief Smooth Type 1 no-offset two-sided deadband function + * + * Smooth approximation to a deadband that returns zero inside the band and + * passes the input through unchanged outside the band. + * + * @tparam ScalarT - scalar data type + * @tparam RealT - Real data type (see GridKit::ScalarTraits::RealT) + * + * @param[in] x - Input signal + * @param[in] lower - Lower breakpoint + * @param[in] upper - Upper breakpoint + * @return Smooth no-offset deadbanded value + */ + template + __attribute__((always_inline)) inline ScalarT deadband1( + const ScalarT x, + const RealT lower, + const RealT upper) + { + assert(lower <= upper); + return x * (sigmoid(lower - x) + sigmoid(x - upper)); + } + + /** + * @brief Smooth Type 2 offset two-sided deadband function * * Smooth approximation to x - min(max(x, lower), upper), composed from the * smooth ramp function. @@ -173,10 +197,10 @@ namespace GridKit * @param[in] x - Input signal * @param[in] lower - Lower breakpoint * @param[in] upper - Upper breakpoint - * @return Smooth deadbanded value + * @return Smooth offset deadbanded value */ template - __attribute__((always_inline)) inline ScalarT deadband( + __attribute__((always_inline)) inline ScalarT deadband2( const ScalarT x, const RealT lower, const RealT upper) diff --git a/GridKit/Model/PhasorDynamics/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..64b81f1b2 100644 --- a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp +++ b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp @@ -7,6 +7,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..bd75dc940 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/Hygov.hpp @@ -0,0 +1,157 @@ +/** + * @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 + 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 + 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::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 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 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 setDerivedParams(); + + ScalarT gatePower(ScalarT gate) const; + RealT invertGatePower(RealT pgv) const; + ScalarT toComponentBase(ScalarT value) const; + ScalarT toSystemBase(ScalarT value) const; + + static constexpr RealT TIME_CONSTANT_MINIMUM = static_cast(1.0e-3); + + RealT Trate_{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_{}; + + RealT va_component_base_{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..5c44a03a4 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovData.hpp @@ -0,0 +1,88 @@ +/** + * @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 + 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 + filter, ///< Governor error filter output + desiredgate, ///< Desired-gate position + gate, ///< Gate position + flow, ///< Turbine flow + head ///< Turbine head + }; + + 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..cc6755968 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovEnzyme.cpp @@ -0,0 +1,90 @@ +/** + * @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 + { + template + int Hygov::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Hygov..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + + if (J_rows_buffer_ == nullptr) + { + auto size = static_cast(size_); + auto signal_size = static_cast(ws_.size()); + auto buffer_size = 2 * size * 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]; + } + + nnz_ = 0; + + GridKit::Enzyme::Sparse::DfDy, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::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, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::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::DfDws, + GridKit::Enzyme::Sparse::MemberFunctions::InternalResidualWithSignal>::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 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..cfb48f7cf --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/HygovImpl.hpp @@ -0,0 +1,603 @@ +/** + * @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 + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + using Log = ::GridKit::Utilities::Logger; + + template + Hygov::Hygov() + { + size_ = static_cast(HygovInternalVariables::MAXIMUM); + setDerivedParams(); + } + + template + Hygov::Hygov(const ModelDataT& data) + : monitor_(std::make_unique(data)) + { + initModelParams(data); + initializeMonitor(); + + size_ = static_cast(HygovInternalVariables::MAXIMUM); + setDerivedParams(); + } + + template + Hygov::~Hygov() + { + } + + template + void Hygov::setDerivedParams() + { + va_component_base_ = Trate_ * static_cast(1.0e6); + } + + template + scalar_type Hygov::toComponentBase(scalar_type value) const + { + return value * va_system_base_ / va_component_base_; + } + + template + scalar_type Hygov::toSystemBase(scalar_type value) const + { + return value / toComponentBase(static_cast(ONE)); + } + + template + void Hygov::initModelParams(const ModelDataT& data) + { + using Params = typename ModelDataT::Parameters; + + Trate_ = ZERO; + Rperm_ = static_cast(0.04); + Rtemp_ = static_cast(0.3); + Tr_ = static_cast(5.0); + Tf_ = static_cast(0.05); + Tg_ = static_cast(0.5); + Velm_ = static_cast(0.2); + Gmax_ = ONE; + Gmin_ = ZERO; + Tw_ = ONE; + At_ = static_cast(1.2); + Dturb_ = static_cast(0.5); + Qnl_ = static_cast(0.05); + Tn_ = ZERO; + Tnp_ = ZERO; + db1_ = ZERO; + db2_ = ZERO; + Hdam_ = ONE; + Gv_.fill(ZERO); + Pgv_.fill(ZERO); + + if (data.parameters.contains(Params::Trate)) + { + Trate_ = std::get(data.parameters.at(Params::Trate)); + } + if (data.parameters.contains(Params::Rperm)) + { + Rperm_ = std::get(data.parameters.at(Params::Rperm)); + } + if (data.parameters.contains(Params::Rtemp)) + { + Rtemp_ = std::get(data.parameters.at(Params::Rtemp)); + } + if (data.parameters.contains(Params::Tr)) + { + Tr_ = std::get(data.parameters.at(Params::Tr)); + } + if (data.parameters.contains(Params::Tf)) + { + Tf_ = std::get(data.parameters.at(Params::Tf)); + } + if (data.parameters.contains(Params::Tg)) + { + Tg_ = std::get(data.parameters.at(Params::Tg)); + } + if (data.parameters.contains(Params::Velm)) + { + Velm_ = std::get(data.parameters.at(Params::Velm)); + } + if (data.parameters.contains(Params::Gmax)) + { + Gmax_ = std::get(data.parameters.at(Params::Gmax)); + } + if (data.parameters.contains(Params::Gmin)) + { + Gmin_ = std::get(data.parameters.at(Params::Gmin)); + } + if (data.parameters.contains(Params::Tw)) + { + Tw_ = std::get(data.parameters.at(Params::Tw)); + } + if (data.parameters.contains(Params::At)) + { + At_ = std::get(data.parameters.at(Params::At)); + } + if (data.parameters.contains(Params::Dturb)) + { + Dturb_ = std::get(data.parameters.at(Params::Dturb)); + } + if (data.parameters.contains(Params::Qnl)) + { + Qnl_ = std::get(data.parameters.at(Params::Qnl)); + } + if (data.parameters.contains(Params::Tn)) + { + Tn_ = std::get(data.parameters.at(Params::Tn)); + } + if (data.parameters.contains(Params::Tnp)) + { + Tnp_ = std::get(data.parameters.at(Params::Tnp)); + } + if (data.parameters.contains(Params::db1)) + { + db1_ = std::get(data.parameters.at(Params::db1)); + } + if (data.parameters.contains(Params::db2)) + { + db2_ = std::get(data.parameters.at(Params::db2)); + } + if (data.parameters.contains(Params::Hdam)) + { + Hdam_ = std::get(data.parameters.at(Params::Hdam)); + } + if (data.parameters.contains(Params::Gv0)) + { + Gv_[0] = std::get(data.parameters.at(Params::Gv0)); + } + if (data.parameters.contains(Params::Gv1)) + { + Gv_[1] = std::get(data.parameters.at(Params::Gv1)); + } + if (data.parameters.contains(Params::Gv2)) + { + Gv_[2] = std::get(data.parameters.at(Params::Gv2)); + } + if (data.parameters.contains(Params::Gv3)) + { + Gv_[3] = std::get(data.parameters.at(Params::Gv3)); + } + if (data.parameters.contains(Params::Gv4)) + { + Gv_[4] = std::get(data.parameters.at(Params::Gv4)); + } + if (data.parameters.contains(Params::Gv5)) + { + Gv_[5] = std::get(data.parameters.at(Params::Gv5)); + } + if (data.parameters.contains(Params::Pgv0)) + { + Pgv_[0] = std::get(data.parameters.at(Params::Pgv0)); + } + if (data.parameters.contains(Params::Pgv1)) + { + Pgv_[1] = std::get(data.parameters.at(Params::Pgv1)); + } + if (data.parameters.contains(Params::Pgv2)) + { + Pgv_[2] = std::get(data.parameters.at(Params::Pgv2)); + } + if (data.parameters.contains(Params::Pgv3)) + { + Pgv_[3] = std::get(data.parameters.at(Params::Pgv3)); + } + if (data.parameters.contains(Params::Pgv4)) + { + Pgv_[4] = std::get(data.parameters.at(Params::Pgv4)); + } + if (data.parameters.contains(Params::Pgv5)) + { + Pgv_[5] = std::get(data.parameters.at(Params::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_; + } + + Tr_ = std::max(Tr_, TIME_CONSTANT_MINIMUM); + Tf_ = std::max(Tf_, TIME_CONSTANT_MINIMUM); + Tg_ = std::max(Tg_, TIME_CONSTANT_MINIMUM); + Tw_ = std::max(Tw_, TIME_CONSTANT_MINIMUM); + Tnp_ = std::max(Tnp_, TIME_CONSTANT_MINIMUM); + + 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::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)]; }); + } + + 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_.resize(size); + y_.resize(size); + yp_.resize(size); + tag_.resize(size); + abs_tol_.resize(size); + variable_indices_.resize(size); + residual_indices_.resize(size); + + auto signal_size = static_cast(HygovExternalVariables::MAXIMUM); + ws_.resize(signal_size); + ws_indices_.resize(signal_size); + std::fill(ws_.begin(), ws_.end(), ScalarT{0}); + std::fill(ws_indices_.begin(), ws_indices_.end(), 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 = 0; + + auto check = [&](bool condition, const char* message) + { + if (!condition) + { + Log::error() << "Hygov: " << message << '\n'; + ret += 1; + } + }; + + check(Trate_ > ZERO, "Trate must be positive"); + check(Rtemp_ != ZERO, "Rtemp must be nonzero"); + check(Velm_ >= ZERO, "Velm must be non-negative"); + check(Gmin_ <= Gmax_, "Gmin must be less than or equal to Gmax"); + check(db1_ >= ZERO, "db1 must be non-negative"); + 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"); + } + + 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 (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 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 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]); + + 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_[XF] = ZERO; + y_[EF] = ZERO; + y_[FC] = ZERO; + y_[RC] = ZERO; + y_[PMECH] = toSystemBase(pmech0); + + const ScalarT yomega = y_[XN] + leadlag_gain_ * (y_[OMEGADB] - y_[XN]); + pref_set_ = y_[EF] - paux_set_ + 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)] = true; + tag_[static_cast(HygovInternalVariables::XF)] = true; + tag_[static_cast(HygovInternalVariables::C)] = true; + tag_[static_cast(HygovInternalVariables::G)] = true; + tag_[static_cast(HygovInternalVariables::Q)] = true; + return 0; + } + + template + int Hygov::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + + template + __attribute__((always_inline)) inline int + Hygov::evaluateInternalResidual( + const ScalarT* y, + const ScalarT* yp, + [[maybe_unused]] const ScalarT* wb, + const 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 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 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 ef = y[EF]; + const ScalarT fc = y[FC]; + const ScalarT rc = y[RC]; + const ScalarT pgv = y[PGV]; + 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]; + + const ScalarT yomega = xn + leadlag_gain_ * (omegadb - xn); + + f[XN] = -yp[XN] + (omegadb - xn) / Tnp_; + f[XF] = -yp[XF] + (ef - xf) / Tf_; + f[C] = -yp[C] + Math::antiwindup(c, rc, Gmin_, Gmax_); + f[G] = -yp[G] + (c - g) / Tg_; + f[Q] = -yp[Q] + (Hdam_ - head) / Tw_; + f[OMEGADB] = -omegadb + Math::deadband1(omega, -db1_, db1_); + f[EF] = -ef + pref + paux - yomega - Rperm_ * c; + f[FC] = -fc + (xf / Tr_ + (ef - xf) / Tf_) / Rtemp_; + f[RC] = -rc + Math::clamp(fc, -Velm_, Velm_); + f[PGV] = -pgv + gatePower(g); + f[H] = -q * q + head * pgv * pgv; + f[PMECH] = -toComponentBase(pmech) + At_ * head * (q - Qnl_) - Dturb_ * omega * g; + + 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..07df81da7 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/HYGOV/README.md @@ -0,0 +1,258 @@ +# **Hydro Turbine-Governor Model (HYGOV)** + +HYGOV is a hydro turbine-governor model with temporary droop, a gate servo, and +a nonlinear single-penstock turbine. It reads machine speed deviation and +outputs mechanical power. + +Notes: +- PowerWorld uses the connected machine base when `Trate = 0`; GridKit requires + `Trate` to be set explicitly. +- HYGOVD `dbL`/`dbH`, `db2` backlash, and Kaplan blade-servo fields are not + modeled. + +## Block Diagram + +Standard HYGOV block diagram. + +![](../../../../../docs/Figures/PhasorDynamics/HYGOV_diagram.png) + +Figure 1: Governor HYGOV model. Figure courtesy of [PowerWorld](https://www.powerworld.com/WebHelp/) + +## Model Parameters + +Symbol | Units | JSON | Description | Typical Value | Note +--------------------------------|----------|----------|----------------------------------------------|---------------|------ +$T_{\mathrm{rate}}$ | [MW] | `Trate` | Component-base turbine rating | 100.0 | Required positive value +$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 +$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 | +$T_{np}$ | [sec] | `Tnp` | Speed lead-lag denominator time constant | 0.0 | +$D_{\omega}$ | [p.u.] | `db1` | Type 1 speed deadband threshold | 0.0 | Uses CommonMath `deadband1` +$H_{\mathrm{dam}}$ | [p.u.] | `Hdam` | Head available at dam | 1.0 | +$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$ + +All-zero `Gv` and `Pgv` source points select the identity curve. + +### Parameter Validation + +Invalid HYGOV 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_r,T_f,T_g,T_w,T_{np}\} \\ + T_{\mathrm{rate}}, H_{\mathrm{dam}} + &> 0 \\ + R_{\mathrm{temp}} + &\ne 0 \\ + V_{\mathrm{elm}}, D_{\omega} + &\ge 0 \\ + G^{\min} &\le G^{\max} \\ + G_V^{(k)} &< G_V^{(k+1)} + \quad k\in\{0,\ldots,4\} \\ + P_{\mathrm{GV}}^{(k)} &\le P_{\mathrm{GV}}^{(k+1)} + \quad k\in\{0,\ldots,4\} +\end{aligned} +``` + +Initialization also requires $N_{\mathrm{GV}}^{-1}$ to be single-valued at the +initial operating point. + +### Model Derived Parameters + +The speed lead-lag coefficient and nonlinear gate-to-power curve are: + +```math +\begin{aligned} + k_n &= + \dfrac{T_n}{T_{np}} \\ + 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} +``` + +CommonMath defines the [linear segment](../../../../CommonMath.md#derived-functions) +helper used by $N_{\mathrm{GV}}$. + +## 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 +$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` +$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)$ +$H$ | [p.u.] | Turbine head | Implicit water-column head +$P_{\text{m}}$ | [p.u.] | Mechanical power to generator | Read by a machine model + +### 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 + +### Differential Equations + +```math +\begin{aligned} + 0 &= + -\dot{x}_n + + \dfrac{1}{T_{np}} + \left(\omega_{\mathrm{db}} - x_n\right) \\ + 0 &= + -\dot{x}_f + + \dfrac{1}{T_f} + \left(e_f - x_f\right) \\ + 0 &= + -\dot{c} + + \text{antiwindup} + \left(c, r_c;\, G^{\min}, G^{\max}\right) \\ + 0 &= + -\dot{g} + + \dfrac{1}{T_g} + \left(c - g\right) \\ + 0 &= + -\dot{q} + + \dfrac{1}{T_w} + \left(H_{\mathrm{dam}} - H\right) +\end{aligned} +``` + +CommonMath defines the [Anti-Windup](../../../../CommonMath.md#antiwindup) +helper. + +### Algebraic Equations + +```math +\begin{aligned} + 0 &= + -\omega_{\mathrm{db}} + + \text{deadband1} + \left(\omega;\, -D_{\omega}, D_{\omega}\right) \\ + 0 &= + -e_f + + P_{\mathrm{ref}} + + P_{\mathrm{aux}} + - x_n + - k_n\left(\omega_{\mathrm{db}} - x_n\right) + - R_{\mathrm{perm}}c \\ + 0 &= + -f_c + + \dfrac{1}{R_{\mathrm{temp}}} + \left[ + \dfrac{x_f}{T_r} + + \dfrac{e_f - x_f}{T_f} + \right] \\ + 0 &= + -r_c + + \text{clamp} + \left(f_c;\, -V_{\mathrm{elm}}, V_{\mathrm{elm}}\right) \\ + 0 &= + -P_{\mathrm{GV}} + + N_{\mathrm{GV}}(g) \\ + 0 &= + -q^2 + + H P_{\mathrm{GV}}^2 \\ + 0 &= + -\dfrac{S_{\mathrm{sys}}}{T_{\mathrm{rate}}}P_{\text{m}} + + \left[ + A_t H\left(q - q_{\mathrm{NL}}\right) + - D_{\mathrm{turb}}\omega g + \right] +\end{aligned} +``` + +CommonMath defines helper targets and smooth approximations for +[deadband1 and clamp](../../../../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: + +```math +\begin{aligned} + \omega_0 &= 0 \\ + P_{\mathrm{aux},0} &= 0 \\ + H_0 &= H_{\mathrm{dam}} \\ + q_0 + &= q_{\mathrm{NL}} + + \dfrac{S_{\mathrm{sys}}P_{\text{m},0}} + {T_{\mathrm{rate}}A_tH_0} \\ + P_{\mathrm{GV},0} &= \dfrac{q_0}{\sqrt{H_0}} \\ + g_0 &= N_{\mathrm{GV}}^{-1}\!\left(P_{\mathrm{GV},0}\right) \\ + c_0 &= g_0 \\ + \omega_{\mathrm{db},0} + &= \text{deadband1}\!\left(\omega_0, -D_{\omega}, D_{\omega}\right) \\ + x_{n,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} + + x_{n,0} + + k_n\left(\omega_{\mathrm{db},0} - x_{n,0}\right) + + R_{\mathrm{perm}}c_0 +\end{aligned} +``` + +This closed-form start requires $G^{\min} \le c_0 \le G^{\max}$ and a +single-valued inverse $N_{\mathrm{GV}}^{-1}(P_{\mathrm{GV},0})$. + +## Monitorable Output + +Output | Units | Description | Note +---------------|--------|-------------------------------------|------ +`pmech` | [p.u.] | Mechanical-power output | $P_{\text{m}}$ +`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 | Implicit water-column head diff --git a/GridKit/Model/PhasorDynamics/Governor/README.md b/GridKit/Model/PhasorDynamics/Governor/README.md index fc46e5eb3..6d3983adb 100644 --- a/GridKit/Model/PhasorDynamics/Governor/README.md +++ b/GridKit/Model/PhasorDynamics/Governor/README.md @@ -9,5 +9,6 @@ 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)) - IEEE Type G1 Turbine Governor (See [IEEEG1](IEEEG1/README.md)) - General Governor (See [GGOV1](GGOV1/README.md)) diff --git a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md index 5bff41ad4..8abd23e7a 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` | `Trate`, `R`, `T1`, `T2`, `T3`, `Pvmax`, `Pvmin`, `Dt` | `none` + `Hygov` | the HYGOV hydro turbine-governor model | `pmech`, `speed`\*, `pref`\*, `paux`\* | `Trate`, `Rperm`, `Rtemp`, `Tr`, `Tf`, `Tg`, `Velm`, `Gmax`, `Gmin`, `Tw`, `At`, `Dturb`, `Qnl`, `Tn`, `Tnp`, `db1`, `db2`, `Hdam`, `Gv0`-`Gv5`, `Pgv0`-`Pgv5` | `pmech`, `filter`, `desiredgate`, `gate`, `flow`, `head` `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/SystemModelData.hpp b/GridKit/Model/PhasorDynamics/SystemModelData.hpp index 8a6e73056..7ce549029 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelData.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelData.hpp @@ -12,6 +12,7 @@ #include #include #include +#include #include #include #include @@ -42,6 +43,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; @@ -106,6 +108,21 @@ namespace GridKit std::vector stabilizer; ///< Stabilizers within the model std::vector constant_source; ///< Constant signal sources within the model std::vector signal; ///< Signal nodes + std::vector bus; ///< Buses within the model + 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 genrou; ///< GENROU instances within the model + std::vector gensal; ///< GENSAL instances within the model + std::vector genclassical; ///< Classical generator instances within the model + std::vector loadz; ///< LoadZ instances within the model + std::vector loadzip; ///< LoadZIP instances 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 + std::vector signal; ///< Signal nodes /// Monitor sink specs std::vector monitor_sink; diff --git a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp index a3682feca..3227aa95c 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/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index 287b328dc..9367531fb 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -259,6 +259,42 @@ namespace GridKit addComponent(gov); } + // Add HYGOV governors + for (const auto& govdata : data.hygov) + { + auto* gov = new Hygov(govdata); + + if (govdata.ports.contains(HygovPorts::speed)) + { + IdxT speed = govdata.ports.at(HygovPorts::speed); + constexpr auto OMEGA = HygovExternalVariables::OMEGA; + gov->getSignals().template attachSignalNode(getSignal(speed)); + } + + if (govdata.ports.contains(HygovPorts::pmech)) + { + IdxT pmech = govdata.ports.at(HygovPorts::pmech); + constexpr auto PMECH = HygovInternalVariables::PMECH; + gov->getSignals().template assignSignalNode(getSignal(pmech)); + } + + if (govdata.ports.contains(HygovPorts::pref)) + { + IdxT pref = govdata.ports.at(HygovPorts::pref); + constexpr auto PREF = HygovExternalVariables::PREF; + gov->getSignals().template attachSignalNode(getSignal(pref)); + } + + if (govdata.ports.contains(HygovPorts::paux)) + { + IdxT paux = govdata.ports.at(HygovPorts::paux); + constexpr auto PAUX = HygovExternalVariables::PAUX; + gov->getSignals().template attachSignalNode(getSignal(paux)); + } + + addComponent(gov); + } + for (const auto& excitedata : data.exciter) { IdxT bus_index = 0; 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/docs/GridKit/Model/PhasorDynamics/Governor/HYGOV/README.md b/docs/GridKit/Model/PhasorDynamics/Governor/HYGOV/README.md new file mode 100644 index 000000000..0e4dd411a --- /dev/null +++ b/docs/GridKit/Model/PhasorDynamics/Governor/HYGOV/README.md @@ -0,0 +1,6 @@ +# HYGOV + +```{include} ../../../../../../GridKit/Model/PhasorDynamics/Governor/HYGOV/README.md +:start-line: 1 +:relative-images: +``` diff --git a/docs/GridKit/Model/PhasorDynamics/Governor/README.md b/docs/GridKit/Model/PhasorDynamics/Governor/README.md index fe3d5b5b7..4e38abd24 100644 --- a/docs/GridKit/Model/PhasorDynamics/Governor/README.md +++ b/docs/GridKit/Model/PhasorDynamics/Governor/README.md @@ -8,6 +8,7 @@ TGOV1 IEEEG1 GGOV1 +HYGOV ``` ```{include} ../../../../../GridKit/Model/PhasorDynamics/Governor/README.md 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..42d64bbf2 100644 --- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt +++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt @@ -88,6 +88,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_systemmodel + GridKit::phasor_dynamics_systemmodel_dependency_tracking + GridKit::testing) + add_executable(test_phasor_exciter_ieeet1 runExciterIeeet1Tests.cpp) target_link_libraries( test_phasor_exciter_ieeet1 @@ -146,6 +154,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) @@ -168,6 +177,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..261298b47 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/GovernorHygovTests.hpp @@ -0,0 +1,415 @@ +#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() + { + using Params = PhasorDynamics::Governor::HygovParameters; + + 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_curve = makeHygovData(); + bad_hygov_curve.parameters[Params::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[Params::Trate] = static_cast(0.0); + PhasorDynamics::Governor::Hygov bad_hygov_trate_model(bad_hygov_trate); + success *= (bad_hygov_trate_model.verify() > 0); + + auto missing_hygov_trate = makeHygovData(); + missing_hygov_trate.parameters.erase(Params::Trate); + PhasorDynamics::Governor::Hygov missing_hygov_trate_model( + missing_hygov_trate); + success *= (missing_hygov_trate_model.verify() > 0); + + auto bad_hygov_hdam = makeHygovData(); + bad_hygov_hdam.parameters[Params::Hdam] = static_cast(0.0); + PhasorDynamics::Governor::Hygov bad_hygov_hdam_model( + bad_hygov_hdam); + success *= (bad_hygov_hdam_model.verify() > 0); + + auto bad_hygov_rtemp = makeHygovData(); + bad_hygov_rtemp.parameters[Params::Rtemp] = static_cast(0.0); + PhasorDynamics::Governor::Hygov bad_hygov_rtemp_model( + bad_hygov_rtemp); + success *= (bad_hygov_rtemp_model.verify() > 0); + + auto bad_hygov_limits = makeHygovData(); + bad_hygov_limits.parameters[Params::Velm] = static_cast(-0.1); + bad_hygov_limits.parameters[Params::db1] = static_cast(-0.1); + PhasorDynamics::Governor::Hygov bad_hygov_limits_model( + bad_hygov_limits); + success *= (bad_hygov_limits_model.verify() > 0); + + auto bad_hygov_gate_limits = makeHygovData(); + bad_hygov_gate_limits.parameters[Params::Gmin] = static_cast(1.1); + PhasorDynamics::Governor::Hygov bad_hygov_gate_limits_model( + bad_hygov_gate_limits); + success *= (bad_hygov_gate_limits_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)] == true); + + 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)] == true); + + 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 zeroTimeConstants() + { + using Params = PhasorDynamics::Governor::HygovParameters; + 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); + + auto data = makeHygovData(); + data.parameters[Params::Tr] = static_cast(0.0); + data.parameters[Params::Tf] = static_cast(0.0); + data.parameters[Params::Tg] = static_cast(0.0); + data.parameters[Params::Tw] = static_cast(0.0); + data.parameters[Params::Tnp] = static_cast(0.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.tagDifferentiable() == 0); + success *= (hygov.evaluateResidual() == 0); + + success *= (hygov.tag()[index(Var::XN)] == true); + success *= (hygov.tag()[index(Var::XF)] == true); + success *= (hygov.tag()[index(Var::C)] == true); + success *= (hygov.tag()[index(Var::G)] == true); + success *= (hygov.tag()[index(Var::Q)] == true); + + 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, "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); + success *= (std::get(trate_param) == static_cast(50.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() == 33); + + 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::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..8eb6e82f2 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/runGovernorHygovTests.cpp @@ -0,0 +1,17 @@ +#include "GovernorHygovTests.hpp" + +int main() +{ + GridKit::Testing::TestingResults result; + + GridKit::Testing::GovernorHygovTests test; + + result += test.constructionAndValidation(); + result += test.signals(); + result += test.sourceDefault(); + result += test.zeroTimeConstants(); + result += test.baseConversion(); + result += test.jsonParseAndSystemAssembly(); + + return result.summary(); +}