diff --git a/CHANGELOG.md b/CHANGELOG.md index 066786e09..ffe21010b 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 `GASTPTI` governor model implementation for PhasorDynamics. ## v0.1 diff --git a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp index e48a8a9c5..0b574fb75 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..2ec5aaa55 100644 --- a/GridKit/Model/PhasorDynamics/Governor/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/Governor/CMakeLists.txt @@ -4,3 +4,4 @@ # ]] add_subdirectory(Tgov1) +add_subdirectory(GASTPTI) diff --git a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/CMakeLists.txt new file mode 100644 index 000000000..41dacae76 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/CMakeLists.txt @@ -0,0 +1,54 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +set(_install_headers GastPti.hpp GastPtiData.hpp) + +if(GRIDKIT_ENABLE_ENZYME) + gridkit_add_library( + phasor_dynamics_governor_gastpti + SOURCES GastPtiEnzyme.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_gastpti + SOURCES GastPti.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_gastpti_dependency_tracking + SOURCES GastPtiDependencyTracking.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_gastpti) +target_link_libraries( + phasor_dynamics_components_dependency_tracking + INTERFACE GridKit::phasor_dynamics_governor_gastpti_dependency_tracking) diff --git a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPti.cpp b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPti.cpp new file mode 100644 index 000000000..e0eabfe22 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPti.cpp @@ -0,0 +1,27 @@ +/** + * @file GastPti.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Non-Enzyme instantiation for the GASTPTI governor model. + */ + +#include "GastPtiImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + template + int GastPti::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for GastPti..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class GastPti; + template class GastPti; + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPti.hpp b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPti.hpp new file mode 100644 index 000000000..8bb9aeebc --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPti.hpp @@ -0,0 +1,131 @@ +/** + * @file GastPti.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Declaration of the GASTPTI governor model. + */ + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + template + class SignalNode; + + namespace Governor + { + /// Internal variables of a `GastPti`. + enum class GastPtiInternalVariables : size_t + { + XVALVE, ///< Fuel-valve state + XFLOW, ///< Fuel-flow state + XTEMP, ///< Exhaust-temperature feedback state + VLOAD, ///< Speed/load fuel demand + VTEMP, ///< Temperature-limit fuel demand + VLV, ///< Low-value gate output + PMECH, ///< Mechanical-power output + MAXIMUM, + }; + + /// External variables of a `GastPti`. + enum class GastPtiExternalVariables : size_t + { + OMEGA, ///< Machine speed deviation + PREF, ///< Active-power/load reference + MAXIMUM, + }; + + template + class GastPti : 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::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 ModelDataT = GastPtiData; + using SignalT = SignalNode; + using MonitorT = Model::VariableMonitor; + + GastPti(); + GastPti(const ModelDataT& data); + ~GastPti() override; + + 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 toComponentBase(ScalarT value) const; + ScalarT toSystemBase(ScalarT value) const; + + RealT R_{0}; + RealT T1_{0}; + RealT T2_{0}; + RealT T3_{0}; + RealT At_{0}; + RealT Kt_{0}; + RealT Vmax_{0}; + RealT Vmin_{0}; + RealT Dturb_{0}; + RealT Trate_{0}; + RealT va_component_base_{0}; + + int parameter_error_count_{0}; + + ScalarT pref_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/GASTPTI/GastPtiData.hpp b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiData.hpp new file mode 100644 index 000000000..1ac6ca6f7 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiData.hpp @@ -0,0 +1,68 @@ +/** + * @file GastPtiData.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Modeling data for the GASTPTI governor model. + */ + +#pragma once + +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + /// Parameter keys for the GASTPTI governor model. + enum class GastPtiParameters + { + R, ///< Permanent droop + T1, ///< Fuel-valve time constant + T2, ///< Fuel-flow time constant + T3, ///< Exhaust-temperature time constant + At, ///< Ambient-temperature load limit + Kt, ///< Exhaust-temperature feedback gain + Vmax, ///< Maximum fuel-valve/turbine-power limit + Vmin, ///< Minimum fuel-valve/turbine-power limit + Dturb, ///< Turbine damping coefficient + Trate ///< Turbine-rating power base + }; + + /// Ports for the GASTPTI governor model. + enum class GastPtiPorts + { + 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 + }; + + /// Variables available through the monitor interface. + enum class GastPtiMonitorableVariables + { + pmech, ///< Mechanical power output + fuelvalve, ///< Fuel-valve state + fuelflow, ///< Fuel-flow state + exhausttemp, ///< Exhaust-temperature feedback state + vload, ///< Speed/load fuel demand + vtemp, ///< Temperature-limit fuel demand + vlv ///< Low-value gate output + }; + + template + struct GastPtiData : public ComponentData + { + GastPtiData() = default; + + using Parameters = GastPtiParameters; + using Ports = GastPtiPorts; + using MonitorableVariables = GastPtiMonitorableVariables; + }; + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiDependencyTracking.cpp new file mode 100644 index 000000000..8c588e794 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiDependencyTracking.cpp @@ -0,0 +1,27 @@ +/** + * @file GastPtiDependencyTracking.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Dependency-tracking instantiations for the GASTPTI governor model. + */ + +#include "GastPtiImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + template + int GastPti::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for GastPti..." << std::endl; + Log::misc() << "Jacobian evaluation is not implemented!" << std::endl; + return 0; + } + + template class GastPti; + template class GastPti; + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiEnzyme.cpp b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiEnzyme.cpp new file mode 100644 index 000000000..10a7313d1 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiEnzyme.cpp @@ -0,0 +1,75 @@ +/** + * @file GastPtiEnzyme.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Enzyme sparse Jacobian for the GASTPTI governor model. + */ + +#include + +#include "GastPtiImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + template + int GastPti::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for GastPti..." << 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::GastPti; + 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_); + + GridKit::Enzyme::Sparse::DfDws::eval(this, + f_.size(), + ws_.size(), + (this->getResidualIndices()).data(), + ws_indices_.data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + J_); + return 0; + } + + template class GastPti; + template class GastPti; + } // namespace Governor + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiImpl.hpp b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiImpl.hpp new file mode 100644 index 000000000..594ea6aa8 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/GastPtiImpl.hpp @@ -0,0 +1,349 @@ +/** + * @file GastPtiImpl.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Definition of the GASTPTI governor model. + */ + +#pragma once + +#include +#include + +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Governor + { + using Log = ::GridKit::Utilities::Logger; + + template + GastPti::GastPti() + { + size_ = static_cast(GastPtiInternalVariables::MAXIMUM); + } + + template + GastPti::GastPti(const ModelDataT& data) + : monitor_(std::make_unique(data)) + { + initModelParams(data); + initializeMonitor(); + size_ = static_cast(GastPtiInternalVariables::MAXIMUM); + } + + template + GastPti::~GastPti() + { + } + + template + scalar_type GastPti::toComponentBase(scalar_type value) const + { + return value * va_system_base_ / va_component_base_; + } + + template + scalar_type GastPti::toSystemBase(scalar_type value) const + { + return value / toComponentBase(static_cast(ONE)); + } + + template + void GastPti::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() << "GastPti: 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() << "GastPti: parameter '" << name << "' must be numeric\n"; + ++parameter_error_count_; + } + }; + + load_required_real(Params::R, R_, "R"); + load_required_real(Params::T1, T1_, "T1"); + load_required_real(Params::T2, T2_, "T2"); + load_required_real(Params::T3, T3_, "T3"); + load_required_real(Params::At, At_, "At"); + load_required_real(Params::Kt, Kt_, "Kt"); + load_required_real(Params::Vmax, Vmax_, "Vmax"); + load_required_real(Params::Vmin, Vmin_, "Vmin"); + load_required_real(Params::Dturb, Dturb_, "Dturb"); + load_required_real(Params::Trate, Trate_, "Trate"); + + va_component_base_ = Trate_ * static_cast(1.0e6); + } + + template + const Model::VariableMonitorBase* GastPti::getMonitor() const + { + return monitor_.get(); + } + + template + void GastPti::initializeMonitor() + { + using Variable = typename ModelDataT::MonitorableVariables; + auto index = [](GastPtiInternalVariables variable) + { + return static_cast(variable); + }; + + monitor_->set(Variable::pmech, [this, index] + { return y_[index(GastPtiInternalVariables::PMECH)]; }); + monitor_->set(Variable::fuelvalve, [this, index] + { return y_[index(GastPtiInternalVariables::XVALVE)]; }); + monitor_->set(Variable::fuelflow, [this, index] + { return y_[index(GastPtiInternalVariables::XFLOW)]; }); + monitor_->set(Variable::exhausttemp, [this, index] + { return y_[index(GastPtiInternalVariables::XTEMP)]; }); + monitor_->set(Variable::vload, [this, index] + { return y_[index(GastPtiInternalVariables::VLOAD)]; }); + monitor_->set(Variable::vtemp, [this, index] + { return y_[index(GastPtiInternalVariables::VTEMP)]; }); + monitor_->set(Variable::vlv, [this, index] + { return y_[index(GastPtiInternalVariables::VLV)]; }); + } + + template + int GastPti::setGridKitComponentID(IdxT component_id) + { + gridkit_component_id_ = component_id; + return 0; + } + + template + int GastPti::allocate() + { + size_ = static_cast(GastPtiInternalVariables::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(GastPtiExternalVariables::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(GastPtiInternalVariables::PMECH)], + &(this->getVariableIndex(static_cast(GastPtiInternalVariables::PMECH)))); + } + + return 0; + } + + template + int GastPti::verify() const + { + int ret = parameter_error_count_; + + auto check = [&](bool condition, const char* message) + { + if (!condition) + { + Log::error() << "GastPti: " << message << '\n'; + ret += 1; + } + }; + + check(R_ > ZERO, "R must be positive"); + check(T1_ >= ZERO, "T1 must be non-negative"); + check(T2_ >= ZERO, "T2 must be non-negative"); + check(T3_ >= ZERO, "T3 must be non-negative"); + check(At_ >= ZERO, "At must be non-negative"); + check(Kt_ >= ZERO, "Kt must be non-negative"); + check(Vmin_ <= Vmax_, "Vmin must be less than or equal to Vmax"); + check(Dturb_ >= ZERO, "Dturb must be non-negative"); + check(Trate_ > ZERO, "Trate must be positive"); + + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "GastPti: omega signal attached with no linked source\n"; + ret += 1; + } + + if (signals_.template isAttached() + && !signals_.template isLinked()) + { + Log::error() << "GastPti: pref signal attached with no linked source\n"; + ret += 1; + } + + return ret; + } + + template + int GastPti::initialize() + { + if (parameter_error_count_ > 0 || verify() > 0) + { + Log::error() << "GastPti: cannot initialize with invalid configuration\n"; + return 1; + } + + const auto XVALVE = static_cast(GastPtiInternalVariables::XVALVE); + const auto XFLOW = static_cast(GastPtiInternalVariables::XFLOW); + const auto XTEMP = static_cast(GastPtiInternalVariables::XTEMP); + const auto VLOAD = static_cast(GastPtiInternalVariables::VLOAD); + const auto VTEMP = static_cast(GastPtiInternalVariables::VTEMP); + const auto VLV = static_cast(GastPtiInternalVariables::VLV); + const auto PMECH = static_cast(GastPtiInternalVariables::PMECH); + + ScalarT omega0{ZERO}; + if (signals_.template isAttached()) + { + omega0 = signals_.template readExternalVariable(); + } + + // Machine initialization seeds the assigned PMECH signal before + // governor initialization reads it. + const ScalarT pmech0 = toComponentBase(y_[PMECH]); + const ScalarT xflow0 = pmech0 + Dturb_ * omega0; + + y_[XFLOW] = xflow0; + y_[XVALVE] = xflow0; + y_[XTEMP] = xflow0; + y_[VTEMP] = At_ + Kt_ * (At_ - y_[XTEMP]); + y_[VLV] = y_[XVALVE]; + y_[VLOAD] = y_[XVALVE]; + y_[PMECH] = toSystemBase(pmech0); + + pref_set_ = y_[VLOAD] + omega0 / R_; + if (signals_.template isAttached()) + { + pref_set_ = signals_.template readExternalVariable(); + y_[VLOAD] = pref_set_ - omega0 / R_; + } + + std::fill(yp_.begin(), yp_.end(), ZERO); + return 0; + } + + template + int GastPti::tagDifferentiable() + { + std::fill(tag_.begin(), tag_.end(), false); + tag_[static_cast(GastPtiInternalVariables::XVALVE)] = (T1_ > ZERO); + tag_[static_cast(GastPtiInternalVariables::XFLOW)] = (T2_ > ZERO); + tag_[static_cast(GastPtiInternalVariables::XTEMP)] = (T3_ > ZERO); + return 0; + } + + template + __attribute__((always_inline)) inline int GastPti::evaluateInternalResidual( + ScalarT* y, + ScalarT* yp, + [[maybe_unused]] ScalarT* wb, + ScalarT* ws, + ScalarT* f) + { + const auto XVALVE = static_cast(GastPtiInternalVariables::XVALVE); + const auto XFLOW = static_cast(GastPtiInternalVariables::XFLOW); + const auto XTEMP = static_cast(GastPtiInternalVariables::XTEMP); + const auto VLOAD = static_cast(GastPtiInternalVariables::VLOAD); + const auto VTEMP = static_cast(GastPtiInternalVariables::VTEMP); + const auto VLV = static_cast(GastPtiInternalVariables::VLV); + const auto PMECH = static_cast(GastPtiInternalVariables::PMECH); + + const auto OMEGA = static_cast(GastPtiExternalVariables::OMEGA); + const auto PREF = static_cast(GastPtiExternalVariables::PREF); + + const ScalarT xvalve = y[XVALVE]; + const ScalarT xflow = y[XFLOW]; + const ScalarT xtemp = y[XTEMP]; + const ScalarT vload = y[VLOAD]; + const ScalarT vtemp = y[VTEMP]; + const ScalarT vlv = y[VLV]; + const ScalarT pmech = y[PMECH]; + + const ScalarT omega = ws[OMEGA]; + const ScalarT pref = ws[PREF]; + + const ScalarT fvalve = vlv - xvalve; + + f[XVALVE] = -T1_ * yp[XVALVE] + Math::antiwindup(xvalve, fvalve, Vmin_, Vmax_); + f[XFLOW] = -T2_ * yp[XFLOW] - xflow + xvalve; + f[XTEMP] = -T3_ * yp[XTEMP] - xtemp + xflow; + f[VLOAD] = -omega + R_ * (pref - vload); + f[VTEMP] = -vtemp + At_ + Kt_ * (At_ - xtemp); + f[VLV] = -vlv + Math::min(vload, vtemp); + const RealT system_base = va_system_base_ / static_cast(1.0e6); + f[PMECH] = -system_base * pmech + Trate_ * (xflow - Dturb_ * omega); + + return 0; + } + + template + int GastPti::evaluateResidual() + { + const auto OMEGA = static_cast(GastPtiExternalVariables::OMEGA); + const auto PREF = static_cast(GastPtiExternalVariables::PREF); + + ws_[OMEGA] = ZERO; + ws_[PREF] = pref_set_; + ws_indices_[OMEGA] = INVALID_INDEX; + ws_indices_[PREF] = 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(); + } + + 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/GASTPTI/README.md b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/README.md new file mode 100644 index 000000000..e898d70d5 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Governor/GASTPTI/README.md @@ -0,0 +1,221 @@ +# **Gas Turbine-Governor Model (GASTPTI)** + +GASTPTI is a gas turbine-governor model for thermal generating units. In +GridKit it is represented as a governor model that reads machine speed +deviation and supplies mechanical power to the machine through a fuel-valve, +fuel-flow, and exhaust-temperature limiting chain. + +Notes: +- GASTPTI uses `Trate` as its turbine-power component base. +- The shared `pmech` signal is treated as system-base power. GridKit converts + between system base and `Trate` at the GASTPTI signal boundary. +- The diagram shows the GASTD speed deadband block (`dbL`/`dbH`). That + block is only used by GASTD; GASTPTI uses the speed-deviation input + $\omega$ directly. + +## Block Diagram + +Standard GASTPTI block diagram. + +
+ + + Figure 1: Governor GASTPTI model. Figure courtesy of [PowerWorld](https://www.powerworld.com/WebHelp/) +
+ +## Model Parameters + +Symbol | Units | JSON | Description | Typical Value | Note +--------------------------------|----------|----------|----------------------------------------------|---------------|------ +$R$ | [p.u.] | `R` | Permanent droop | 0.05 | Required nonzero +$T_1$ | [sec] | `T1` | Fuel-valve time constant | 0.4 | If zero, $x_{\mathrm{valve}}$ is algebraic +$T_2$ | [sec] | `T2` | Fuel-flow time constant | 0.1 | If zero, $x_{\mathrm{flow}}$ is algebraic +$T_3$ | [sec] | `T3` | Exhaust-temperature time constant | 3.0 | If zero, $x_{\mathrm{temp}}$ is algebraic +$A_{\mathrm{T}}$ | [p.u.] | `At` | Ambient-temperature load limit | 1.0 | +$K_{\mathrm{T}}$ | [p.u.] | `Kt` | Exhaust-temperature feedback gain | 2.0 | +$V^{\max}$ | [p.u.] | `Vmax` | Maximum fuel-valve/turbine-power limit | 1.0 | +$V^{\min}$ | [p.u.] | `Vmin` | Minimum fuel-valve/turbine-power limit | 0.0 | +$D_{\mathrm{turb}}$ | [p.u.] | `Dturb` | Turbine damping coefficient | 0.0 | Multiplies speed deviation +$T_{\mathrm{rate}}$ | [MW] | `Trate` | Turbine-rating power base | 100.0 | Required positive value + +### Parameter Validation + +Invalid GASTPTI parameter sets are rejected by the following checks. + +The required checks are: + +```math +\begin{aligned} + &R > 0 \\ + &T_1, T_2, T_3 \ge 0 \\ + &A_{\mathrm{T}} \ge 0 \\ + &K_{\mathrm{T}} \ge 0 \\ + &V^{\min} \le V^{\max} \\ + &D_{\mathrm{turb}} \ge 0 \\ + &T_{\mathrm{rate}} > 0 +\end{aligned} +``` + +### Model Derived Parameters + +The system power base $S_{\mathrm{sys}}$ is derived from the system model. + +## Model Variables + +### Internal Variables + +#### Differential + +Symbol | Units | Description | Note +------------------------|--------|------------------------------------|------ +$x_{\mathrm{valve}}$ | [p.u.] | Fuel-valve state | State 1 in Fig. 1; source label: `Fuel Valve`; algebraic when $T_1 = 0$ +$x_{\mathrm{flow}}$ | [p.u.] | Fuel-flow state | State 2 in Fig. 1; source label: `Fuel Flow`; algebraic when $T_2 = 0$ +$x_{\mathrm{temp}}$ | [p.u.] | Exhaust-temperature feedback state | State 3 in Fig. 1; source label: `Exhaust Temperature`; algebraic when $T_3 = 0$ + +#### Algebraic + +Symbol | Units | Description | Note +--------------------------------|--------|--------------------------------------|------ +$V_{\mathrm{load}}$ | [p.u.] | Speed/load fuel demand | $P_{\mathrm{ref}}$ less droop feedback; LV gate input +$V_{\mathrm{temp}}$ | [p.u.] | Temperature-limit fuel demand | Exhaust-temperature branch; LV gate input +$V_{\mathrm{LV}}$ | [p.u.] | LV gate output | Lesser of $V_{\mathrm{load}}$ and $V_{\mathrm{temp}}$; drives the fuel-valve lag +$P_{\text{m}}$ | [p.u.] | Mechanical power to generator | On system base + +### 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 | On `Trate` component base; source label: `Pref` + +## Model Equations + +### Differential Equations + +The lag residuals are written in descriptor form. When $T_1$, $T_2$, or +$T_3$ is zero, the corresponding residual becomes algebraic. + +```math +\begin{aligned} + 0 &= + -T_1 \dot x_{\mathrm{valve}} + + \text{antiwindup}( + x_{\mathrm{valve}}, + V_{\mathrm{LV}} - x_{\mathrm{valve}}, + V^{\min}, + V^{\max} + ) \\ + 0 &= -T_2 \dot x_{\mathrm{flow}} - x_{\mathrm{flow}} + x_{\mathrm{valve}} \\ + 0 &= -T_3 \dot x_{\mathrm{temp}} - x_{\mathrm{temp}} + x_{\mathrm{flow}} +\end{aligned} +``` + +CommonMath defines the [Anti-Windup](../../../../CommonMath.md#anti-windup-indicator) +target and smooth approximation. + +### Algebraic Equations + +The algebraic targets form the two LV gate inputs, select between them with the +low-value gate, and assemble the mechanical power output. + +```math +\begin{aligned} + 0 &= -\omega + R(P_{\mathrm{ref}} - V_{\mathrm{load}}) \\ + 0 &= -V_{\mathrm{temp}} + + A_{\mathrm{T}} + + K_{\mathrm{T}}(A_{\mathrm{T}} - x_{\mathrm{temp}}) \\ + 0 &= -V_{\mathrm{LV}} + \min(V_{\mathrm{load}}, V_{\mathrm{temp}}) \\ + 0 &= -S_{\mathrm{sys}} P_{\text{m}} + + T_{\mathrm{rate}}(x_{\mathrm{flow}} - D_{\mathrm{turb}}\omega) +\end{aligned} +``` + +CommonMath defines the helper targets and smooth approximations for +[min](../../../../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. A standard power-flow start uses synchronous machine +speed: + +```math +\begin{aligned} + \omega_0 &= 0 +\end{aligned} +``` + +The machine initializes mechanical power first; GASTPTI reads it as +$P_{\text{m},0}$ on the system base and converts it to the `Trate` +component base: + +```math +\begin{aligned} + S_{\mathrm{sys}} P_{\text{m},0} + &= T_{\mathrm{rate}}P^{\mathrm{gastpti}}_{\text{m},0} +\end{aligned} +``` + +The output residual fixes the fuel-flow state, and the two plain lag residuals +pass that value through to the fuel valve and exhaust-temperature states: + +```math +\begin{aligned} + x_{\mathrm{flow},0} &= + P^{\mathrm{gastpti}}_{\mathrm{mech},0} + + D_{\mathrm{turb}}\omega_0 \\ + x_{\mathrm{valve},0} &= x_{\mathrm{flow},0} \\ + x_{\mathrm{temp},0} &= x_{\mathrm{flow},0} +\end{aligned} +``` + +Then evaluate the temperature-limit branch and the LV-gate target: + +```math +\begin{aligned} + V_{\mathrm{temp},0} + &= A_{\mathrm{T}} + + K_{\mathrm{T}}(A_{\mathrm{T}} - x_{\mathrm{temp},0}) \\ + V_{\mathrm{LV},0} &= x_{\mathrm{valve},0} +\end{aligned} +``` + +For an unsaturated start where the speed/load demand is selected by the LV +gate, choose or verify: + +```math +\begin{aligned} + V_{\mathrm{load},0} &= x_{\mathrm{valve},0} \\ + P_{\mathrm{ref},0} + &= V_{\mathrm{load},0} + \dfrac{\omega_0}{R} +\end{aligned} +``` + +This start requires $V^{\min} \le x_{\mathrm{valve},0} \le V^{\max}$ and +$V_{\mathrm{temp},0} \ge V_{\mathrm{load},0}$ so the LV gate selects the +speed/load demand and the fuel-valve anti-windup residual holds. If the +supplied $P_{\mathrm{ref}}$ is connected, use its initial value and verify the +same residuals. + +Starts that bind the fuel-valve limit or the temperature gate are outside these +closed-form initialization formulas. The zero-derivative lag relations still +require $x_{\mathrm{temp},0}=x_{\mathrm{flow},0}=x_{\mathrm{valve},0}$. + +## Model Outputs + +Output | Units | Description | Note +---------------|--------|------------------------------------|------ +`pmech` | [p.u.] | Mechanical power output | On system base; oriented as mechanical input to the machine +`fuelvalve` | [p.u.] | Fuel-valve state | State 1 +`fuelflow` | [p.u.] | Fuel-flow state | State 2 +`exhausttemp` | [p.u.] | Exhaust-temperature feedback state | State 3 +`vload` | [p.u.] | Speed/load fuel demand | LV gate input +`vtemp` | [p.u.] | Temperature-limit fuel demand | LV gate input +`vlv` | [p.u.] | LV gate output | Selected fuel demand diff --git a/GridKit/Model/PhasorDynamics/Governor/README.md b/GridKit/Model/PhasorDynamics/Governor/README.md index 8e6305aa8..b0a5c5a40 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)) +- Gas Turbine Governor (See [GASTPTI](GASTPTI/README.md)) diff --git a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md index 32df36cdc..40925c7e9 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` + `GastPti` | the GASTPTI gas turbine-governor model | `pmech`, `speed`\*, `pref`\* | `R`, `T1`, `T2`, `T3`, `At`, `Kt`, `Vmax`, `Vmin`, `Dturb`, `Trate` | `pmech`, `fuelvalve`, `fuelflow`, `exhausttemp`, `vload`, `vtemp`, `vlv` `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..449c3c2da 100644 --- a/GridKit/Model/PhasorDynamics/SystemModel.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModel.hpp @@ -297,6 +297,33 @@ namespace GridKit addComponent(gov); } + // Add GASTPTI governors + for (const auto& govdata : data.gastpti) + { + auto* gov = new GastPti(govdata); + + if (govdata.ports.contains(GastPtiData::Ports::speed)) + { + IdxT speed = govdata.ports.at(GastPtiData::Ports::speed); + gov->getSignals().template attachSignalNode(getSignal(speed)); + } + + if (govdata.ports.contains(GastPtiData::Ports::pmech)) + { + IdxT pmech = govdata.ports.at(GastPtiData::Ports::pmech); + gov->getSignals().template assignSignalNode(getSignal(pmech)); + } + + if (govdata.ports.contains(GastPtiData::Ports::pref)) + { + IdxT pref = govdata.ports.at(GastPtiData::Ports::pref); + gov->getSignals().template attachSignalNode( + getSignal(pref)); + } + + 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..e415f5102 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 GastPtiDataT = Governor::GastPtiData; 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 gastpti; ///< GASTPTI 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..7d50cc498 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 == "GastPti") + { + typename SystemModelData::GastPtiDataT gov; + raw_component.get_to(gov); + sm.gastpti.push_back(gov); + } else if (kind == "Ieeet1") { typename SystemModelData::Ieeet1DataT exciter; diff --git a/docs/Figures/PhasorDynamics/GASTPTI_diagram.png b/docs/Figures/PhasorDynamics/GASTPTI_diagram.png new file mode 100644 index 000000000..3bc194bf7 Binary files /dev/null and b/docs/Figures/PhasorDynamics/GASTPTI_diagram.png differ diff --git a/tests/UnitTests/PhasorDynamics/CMakeLists.txt b/tests/UnitTests/PhasorDynamics/CMakeLists.txt index 932cb13b8..465daeba5 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_gastpti runGovernorGastPtiTests.cpp) +target_link_libraries( + test_phasor_governor_gastpti + 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 PhasorDynamicsGovernorGastPtiTest COMMAND test_phasor_governor_gastpti) 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_gastpti test_phasor_exciter_ieeet1 test_phasor_gensal test_phasor_exciter_sexspti diff --git a/tests/UnitTests/PhasorDynamics/GovernorGastPtiTests.hpp b/tests/UnitTests/PhasorDynamics/GovernorGastPtiTests.hpp new file mode 100644 index 000000000..c13805984 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/GovernorGastPtiTests.hpp @@ -0,0 +1,572 @@ +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + template + class GovernorGastPtiTests + { + public: + using ScalarT = scalar_type; + using IdxT = index_type; + using RealT = typename PhasorDynamics::Component::RealT; + using Gov = PhasorDynamics::Governor::GastPti; + using Data = PhasorDynamics::Governor::GastPtiData; + using Var = PhasorDynamics::Governor::GastPtiInternalVariables; + using Ext = PhasorDynamics::Governor::GastPtiExternalVariables; + using Params = PhasorDynamics::Governor::GastPtiParameters; + using Mon = PhasorDynamics::Governor::GastPtiMonitorableVariables; + + GovernorGastPtiTests() = default; + ~GovernorGastPtiTests() = default; + + TestOutcome constructor() + { + TestStatus success = true; + + Gov model(makeTestData()); + + const auto* monitor = model.getMonitor(); + success *= (model.size() == static_cast(Var::MAXIMUM)); + success *= (monitor != nullptr); + if (monitor != nullptr) + { + success *= (!monitor->empty()); + } + success *= (model.verify() == 0); + + return success.report(__func__); + } + + TestOutcome zeroInitialResidual() + { + TestStatus success = true; + + Gov model(makeTestData()); + + PhasorDynamics::SignalNode pmech_node; + PhasorDynamics::SignalNode omega_node; + const ScalarT pmech0 = scalar(kInitialPmech); + ScalarT pmech_value{0.0}; + ScalarT omega_value = scalar(kInitialOmega); + IdxT pmech_index = INVALID_INDEX; + IdxT omega_index = 9; + pmech_node.set(&pmech_value, &pmech_index); + omega_node.set(&omega_value, &omega_index); + + model.getSignals().template assignSignalNode(&pmech_node); + model.getSignals().template attachSignalNode(&omega_node); + + success *= (model.allocate() == 0); + pmech_node.init(pmech0); + success *= pmech_node.linked(); + success *= (pmech_node.getVariableIndex() == static_cast(Var::PMECH)); + + success *= (model.verify() == 0); + success *= (model.initialize() == 0); + success *= (model.tagDifferentiable() == 0); + success *= (model.evaluateResidual() == 0); + + const ScalarT xflow0 = pmech0 + scalar(kDturb) * omega_value; + const ScalarT vtemp0 = scalar(kAt) + scalar(kKt) * (scalar(kAt) - xflow0); + + success *= isEqual(model.y()[index(Var::XVALVE)], xflow0, scalar(kTolerance)); + success *= isEqual(model.y()[index(Var::XFLOW)], xflow0, scalar(kTolerance)); + success *= isEqual(model.y()[index(Var::XTEMP)], xflow0, scalar(kTolerance)); + success *= isEqual(model.y()[index(Var::VLOAD)], xflow0, scalar(kTolerance)); + success *= isEqual(model.y()[index(Var::VTEMP)], vtemp0, scalar(kTolerance)); + success *= isEqual(model.y()[index(Var::VLV)], xflow0, scalar(kTolerance)); + success *= (model.tag()[index(Var::XVALVE)] == true); + success *= (model.tag()[index(Var::XFLOW)] == true); + success *= (model.tag()[index(Var::XTEMP)] == true); + + checkZeroResidual(model, success); + + return success.report(__func__); + } + + TestOutcome baseConversion() + { + TestStatus success = true; + + auto data = makeTestData(); + data.parameters[Params::Trate] = static_cast(kConversionTrate); + Gov model(data); + model.setSystemBase(static_cast(kSystemFrequency), + static_cast(kConversionSystemBase * 1.0e6)); + + PhasorDynamics::SignalNode pmech_node; + PhasorDynamics::SignalNode omega_node; + const ScalarT pmech0 = scalar(kConversionInitialPmech); + ScalarT pmech_value{0.0}; + ScalarT omega_value = scalar(kInitialOmega); + IdxT pmech_index = INVALID_INDEX; + IdxT omega_index = 9; + pmech_node.set(&pmech_value, &pmech_index); + omega_node.set(&omega_value, &omega_index); + + model.getSignals().template assignSignalNode(&pmech_node); + model.getSignals().template attachSignalNode(&omega_node); + + success *= (model.allocate() == 0); + pmech_node.init(pmech0); + success *= (model.verify() == 0); + success *= (model.initialize() == 0); + success *= (model.evaluateResidual() == 0); + + const ScalarT pmech_component = + pmech0 * scalar(kConversionSystemBase / kConversionTrate); + const ScalarT xflow0 = pmech_component + scalar(kDturb) * omega_value; + const ScalarT vtemp0 = scalar(kAt) + scalar(kKt) * (scalar(kAt) - xflow0); + + success *= isEqual(model.y()[index(Var::XVALVE)], xflow0, scalar(kTolerance)); + success *= isEqual(model.y()[index(Var::XFLOW)], xflow0, scalar(kTolerance)); + success *= isEqual(model.y()[index(Var::XTEMP)], xflow0, scalar(kTolerance)); + success *= isEqual(model.y()[index(Var::VLOAD)], xflow0, scalar(kTolerance)); + success *= isEqual(model.y()[index(Var::VTEMP)], vtemp0, scalar(kTolerance)); + success *= isEqual(model.y()[index(Var::VLV)], xflow0, scalar(kTolerance)); + success *= isEqual(model.y()[index(Var::PMECH)], pmech0, scalar(kTolerance)); + + checkZeroResidual(model, success); + + return success.report(__func__); + } + + TestOutcome prefSignal() + { + TestStatus success = true; + + Gov model(makeTestData()); + + PhasorDynamics::SignalNode pmech_node; + PhasorDynamics::SignalNode omega_node; + PhasorDynamics::SignalNode pref_node; + const ScalarT pmech0 = scalar(kInitialPmech); + ScalarT pmech_value{0.0}; + ScalarT omega_value = scalar(kInitialOmega); + ScalarT pref_value = prefForInitialPoint(pmech0, omega_value); + IdxT pmech_index = INVALID_INDEX; + IdxT omega_index = 9; + IdxT pref_index = 10; + pmech_node.set(&pmech_value, &pmech_index); + omega_node.set(&omega_value, &omega_index); + pref_node.set(&pref_value, &pref_index); + + model.getSignals().template assignSignalNode(&pmech_node); + model.getSignals().template attachSignalNode(&omega_node); + model.getSignals().template attachSignalNode(&pref_node); + + success *= (model.allocate() == 0); + pmech_node.init(pmech0); + success *= (model.verify() == 0); + success *= (model.initialize() == 0); + success *= (model.evaluateResidual() == 0); + checkZeroResidual(model, success); + + pref_value += scalar(kPrefStep); + success *= (model.evaluateResidual() == 0); + success *= isEqual(model.getResidual()[index(Var::VLOAD)], + scalar(kR * kPrefStep), + scalar(kTolerance)); + + return success.report(__func__); + } + + TestOutcome residual() + { + TestStatus success = true; + + Gov model(makeTestData()); + + PhasorDynamics::SignalNode omega_node; + PhasorDynamics::SignalNode pref_node; + ScalarT omega_value = scalar(kResidualOmega); + ScalarT pref_value = scalar(kResidualPref); + IdxT omega_index = 7; + IdxT pref_index = 8; + omega_node.set(&omega_value, &omega_index); + pref_node.set(&pref_value, &pref_index); + + model.getSignals().template attachSignalNode(&omega_node); + model.getSignals().template attachSignalNode(&pref_node); + + success *= (model.allocate() == 0); + + const ScalarT xvalve = scalar(kResidualXvalve); + const ScalarT xflow = scalar(kResidualXflow); + const ScalarT xtemp = scalar(kResidualXtemp); + const ScalarT vload = scalar(kResidualVload); + const ScalarT vtemp = scalar(kResidualVtemp); + const ScalarT vlv = scalar(kResidualVlv); + const ScalarT pmech = scalar(kResidualPmech); + const ScalarT xvalve_dot = scalar(kResidualXvalveDot); + const ScalarT xflow_dot = scalar(kResidualXflowDot); + const ScalarT xtemp_dot = scalar(kResidualXtempDot); + + model.y()[index(Var::XVALVE)] = xvalve; + model.y()[index(Var::XFLOW)] = xflow; + model.y()[index(Var::XTEMP)] = xtemp; + model.y()[index(Var::VLOAD)] = vload; + model.y()[index(Var::VTEMP)] = vtemp; + model.y()[index(Var::VLV)] = vlv; + model.y()[index(Var::PMECH)] = pmech; + + model.yp()[index(Var::XVALVE)] = xvalve_dot; + model.yp()[index(Var::XFLOW)] = xflow_dot; + model.yp()[index(Var::XTEMP)] = xtemp_dot; + + success *= (model.verify() == 0); + success *= (model.evaluateResidual() == 0); + + const ScalarT valve_target = vlv - xvalve; + const ScalarT selected_vlv = vload; + const std::vector expected = { + -scalar(kT1) * xvalve_dot + valve_target, + -scalar(kT2) * xflow_dot - xflow + xvalve, + -scalar(kT3) * xtemp_dot - xtemp + xflow, + -omega_value + scalar(kR) * (pref_value - vload), + -vtemp + scalar(kAt) + scalar(kKt) * (scalar(kAt) - xtemp), + -vlv + selected_vlv, + -scalar(kSystemBase) * pmech + scalar(kTrate) * (xflow - scalar(kDturb) * omega_value), + }; + + checkResidual(model, expected, success); + + return success.report(__func__); + } + + TestOutcome antiWindupLimiter() + { + TestStatus success = true; + + Gov model(makeTestData()); + success *= (model.allocate() == 0); + + auto check_valve = [&](ScalarT xvalve, ScalarT vlv, ScalarT expected) + { + std::fill(model.y().begin(), model.y().end(), ScalarT{0}); + std::fill(model.yp().begin(), model.yp().end(), ScalarT{0}); + model.y()[index(Var::XVALVE)] = xvalve; + model.y()[index(Var::VLV)] = vlv; + success *= (model.evaluateResidual() == 0); + success *= isEqual(model.getResidual()[index(Var::XVALVE)], + expected, + scalar(kSmoothTolerance)); + }; + + check_valve(scalar(2.2), scalar(3.2), scalar(0.0)); + check_valve(scalar(2.2), scalar(1.2), scalar(-1.0)); + check_valve(scalar(-1.0), scalar(-2.0), scalar(0.0)); + check_valve(scalar(-1.0), scalar(0.0), scalar(1.0)); + + return success.report(__func__); + } + + TestOutcome timeConstantTags() + { + TestStatus success = true; + + auto data = makeTestData(); + data.parameters[Params::T1] = static_cast(0.0); + data.parameters[Params::T2] = static_cast(kT2); + data.parameters[Params::T3] = static_cast(0.0); + + Gov model(data); + success *= (model.allocate() == 0); + success *= (model.tagDifferentiable() == 0); + success *= (model.tag()[index(Var::XVALVE)] == false); + success *= (model.tag()[index(Var::XFLOW)] == true); + success *= (model.tag()[index(Var::XTEMP)] == false); + + return success.report(__func__); + } + + TestOutcome parameterValidation() + { + TestStatus success = true; + + auto missing = makeTestData(); + missing.parameters.erase(Params::R); + Gov missing_model(missing); + success *= (missing_model.verify() > 0); + + auto missing_trate = makeTestData(); + missing_trate.parameters.erase(Params::Trate); + Gov missing_trate_model(missing_trate); + success *= (missing_trate_model.verify() > 0); + + auto negative_time = makeTestData(); + negative_time.parameters[Params::T2] = static_cast(-kT2); + Gov negative_time_model(negative_time); + success *= (negative_time_model.verify() > 0); + + auto invalid_limits = makeTestData(); + invalid_limits.parameters[Params::Vmin] = static_cast(kVmax + 0.1); + invalid_limits.parameters[Params::Vmax] = static_cast(kVmax); + Gov invalid_limits_model(invalid_limits); + success *= (invalid_limits_model.verify() > 0); + + auto invalid_trate = makeTestData(); + invalid_trate.parameters[Params::Trate] = true; + Gov invalid_trate_model(invalid_trate); + success *= (invalid_trate_model.verify() > 0); + + auto zero_trate = makeTestData(); + zero_trate.parameters[Params::Trate] = static_cast(0.0); + Gov zero_trate_model(zero_trate); + success *= (zero_trate_model.verify() > 0); + + return success.report(__func__); + } + + TestOutcome signalValidation() + { + TestStatus success = true; + + PhasorDynamics::SignalNode omega_node; + Gov omega_model(makeTestData()); + omega_model.getSignals().template attachSignalNode(&omega_node); + success *= (omega_model.verify() > 0); + + PhasorDynamics::SignalNode pref_node; + Gov pref_model(makeTestData()); + pref_model.getSignals().template attachSignalNode(&pref_node); + success *= (pref_model.verify() > 0); + + return success.report(__func__); + } + +#ifdef GRIDKIT_ENABLE_ENZYME + TestOutcome jacobian() + { + TestStatus success = true; + + auto dependency_tracking_jacobian = dependencyTrackingJacobian(); + auto enzyme_jacobian = enzymeJacobian(); + + success *= (dependency_tracking_jacobian.size() == enzyme_jacobian.size()); + for (size_t i = 0; i < dependency_tracking_jacobian.size(); ++i) + { + success *= isEqual(dependency_tracking_jacobian[i], enzyme_jacobian[i], kTolerance); + } + + return success.report(__func__); + } +#endif + + private: + static constexpr RealT kTolerance = 100.0 * std::numeric_limits::epsilon(); + static constexpr RealT kSmoothTolerance = 1.0e-2; + + static constexpr RealT kR = 0.05; + static constexpr RealT kT1 = 0.4; + static constexpr RealT kT2 = 0.5; + static constexpr RealT kT3 = 0.25; + static constexpr RealT kAt = 2.0; + static constexpr RealT kKt = 0.3; + static constexpr RealT kVmax = 1.2; + static constexpr RealT kVmin = 0.0; + static constexpr RealT kDturb = 0.1; + static constexpr RealT kTrate = 100.0; + static constexpr RealT kSystemBase = 100.0; + + static constexpr RealT kInitialPmech = 0.75; + static constexpr RealT kInitialOmega = 0.02; + static constexpr RealT kPrefStep = 0.1; + + static constexpr RealT kSystemFrequency = 60.0; + static constexpr RealT kConversionTrate = 50.0; + static constexpr RealT kConversionSystemBase = 100.0; + static constexpr RealT kConversionInitialPmech = 0.40; + + static constexpr RealT kResidualOmega = 0.02; + static constexpr RealT kResidualPref = 1.25; + static constexpr RealT kResidualXvalve = 0.7; + static constexpr RealT kResidualXflow = 0.6; + static constexpr RealT kResidualXtemp = 0.5; + static constexpr RealT kResidualVload = 0.9; + static constexpr RealT kResidualVtemp = 2.5; + static constexpr RealT kResidualVlv = 0.8; + static constexpr RealT kResidualPmech = 0.55; + static constexpr RealT kResidualXvalveDot = 0.05; + static constexpr RealT kResidualXflowDot = -0.1; + static constexpr RealT kResidualXtempDot = 0.2; + + static ScalarT scalar(RealT value) + { + return static_cast(value); + } + + template + static value_type value(RealT value) + { + return value_type{value}; + } + + static size_t index(Var variable) + { + return static_cast(variable); + } + + template + static value_type prefForInitialPoint(const value_type& pmech, const value_type& omega) + { + return pmech + value(kDturb) * omega + omega / value(kR); + } + + void checkResidual(const Gov& model, + const std::vector& expected, + TestStatus& success) const + { + const auto& residual = model.getResidual(); + for (size_t i = 0; i < expected.size(); ++i) + { + if (!isEqual(residual[i], expected[i], scalar(kTolerance))) + { + std::cout << "Unexpected GASTPTI residual at index " << i << ": " + << residual[i] << " != " << expected[i] << "\n"; + success = false; + } + } + } + + void checkZeroResidual(const Gov& model, TestStatus& success) const + { + std::vector expected(static_cast(Var::MAXIMUM), ScalarT{0}); + checkResidual(model, expected, success); + } + + Data makeTestData() + { + Data data; + data.device_class = "GastPti"; + data.disambiguation_string = "gastpti_test"; + data.monitored_variables.insert(Mon::pmech); + data.monitored_variables.insert(Mon::fuelvalve); + + data.parameters[Params::R] = static_cast(kR); + data.parameters[Params::T1] = static_cast(kT1); + data.parameters[Params::T2] = static_cast(kT2); + data.parameters[Params::T3] = static_cast(kT3); + data.parameters[Params::At] = static_cast(kAt); + data.parameters[Params::Kt] = static_cast(kKt); + data.parameters[Params::Vmax] = static_cast(kVmax); + data.parameters[Params::Vmin] = static_cast(kVmin); + data.parameters[Params::Dturb] = static_cast(kDturb); + data.parameters[Params::Trate] = static_cast(kTrate); + + return data; + } + +#ifdef GRIDKIT_ENABLE_ENZYME + using DependencyMap = DependencyTracking::Variable::DependencyMap; + + std::vector dependencyTrackingJacobian() + { + using ADScalarT = DependencyTracking::Variable; + using ADGov = PhasorDynamics::Governor::GastPti; + + ADGov model(makeTestData()); + + PhasorDynamics::SignalNode omega_node; + PhasorDynamics::SignalNode pref_node; + ADScalarT omega_value{kResidualOmega}; + ADScalarT pref_value{kResidualPref}; + IdxT omega_index = static_cast(Var::MAXIMUM); + IdxT pref_index = omega_index + 1; + omega_node.set(&omega_value, &omega_index); + pref_node.set(&pref_value, &pref_index); + + model.getSignals().template attachSignalNode(&omega_node); + model.getSignals().template attachSignalNode(&pref_node); + model.allocate(); + model.updateTime(0.0, 1.0); + + model.y()[index(Var::XVALVE)] = ADScalarT{kResidualXvalve}; + model.y()[index(Var::XFLOW)] = ADScalarT{kResidualXflow}; + model.y()[index(Var::XTEMP)] = ADScalarT{kResidualXtemp}; + model.y()[index(Var::VLOAD)] = ADScalarT{kResidualVload}; + model.y()[index(Var::VTEMP)] = ADScalarT{kResidualVtemp}; + model.y()[index(Var::VLV)] = ADScalarT{kResidualVlv}; + model.y()[index(Var::PMECH)] = ADScalarT{kResidualPmech}; + + model.yp()[index(Var::XVALVE)] = ADScalarT{kResidualXvalveDot}; + model.yp()[index(Var::XFLOW)] = ADScalarT{kResidualXflowDot}; + model.yp()[index(Var::XTEMP)] = ADScalarT{kResidualXtempDot}; + + for (size_t i = 0; i < model.y().size(); ++i) + { + model.y()[i].setVariableNumber(i); + model.yp()[i].setVariableNumber(i); + } + omega_value.setVariableNumber(static_cast(omega_index)); + pref_value.setVariableNumber(static_cast(pref_index)); + + model.evaluateResidual(); + + std::vector dependencies(model.getResidual().size()); + for (size_t i = 0; i < dependencies.size(); ++i) + { + dependencies[i] = model.getResidual()[i].getDependencies(); + } + + return dependencies; + } + + std::vector enzymeJacobian() + { + Gov model(makeTestData()); + + PhasorDynamics::SignalNode omega_node; + PhasorDynamics::SignalNode pref_node; + ScalarT omega_value = scalar(kResidualOmega); + ScalarT pref_value = scalar(kResidualPref); + IdxT omega_index = static_cast(Var::MAXIMUM); + IdxT pref_index = omega_index + 1; + omega_node.set(&omega_value, &omega_index); + pref_node.set(&pref_value, &pref_index); + + model.getSignals().template attachSignalNode(&omega_node); + model.getSignals().template attachSignalNode(&pref_node); + model.allocate(); + model.updateTime(0.0, 1.0); + + model.y()[index(Var::XVALVE)] = scalar(kResidualXvalve); + model.y()[index(Var::XFLOW)] = scalar(kResidualXflow); + model.y()[index(Var::XTEMP)] = scalar(kResidualXtemp); + model.y()[index(Var::VLOAD)] = scalar(kResidualVload); + model.y()[index(Var::VTEMP)] = scalar(kResidualVtemp); + model.y()[index(Var::VLV)] = scalar(kResidualVlv); + model.y()[index(Var::PMECH)] = scalar(kResidualPmech); + + model.yp()[index(Var::XVALVE)] = scalar(kResidualXvalveDot); + model.yp()[index(Var::XFLOW)] = scalar(kResidualXflowDot); + model.yp()[index(Var::XTEMP)] = scalar(kResidualXtempDot); + + model.evaluateResidual(); + model.evaluateJacobian(); + + auto& model_jacobian = model.getJacobian(); + model_jacobian.deduplicate(); + + return MapFromCOO(model_jacobian); + } +#endif + }; + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/runGovernorGastPtiTests.cpp b/tests/UnitTests/PhasorDynamics/runGovernorGastPtiTests.cpp new file mode 100644 index 000000000..3932dd4a8 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/runGovernorGastPtiTests.cpp @@ -0,0 +1,23 @@ +#include "GovernorGastPtiTests.hpp" + +int main() +{ + GridKit::Testing::TestingResults result; + + GridKit::Testing::GovernorGastPtiTests test; + + result += test.constructor(); + result += test.zeroInitialResidual(); + result += test.baseConversion(); + result += test.prefSignal(); + result += test.residual(); + result += test.antiWindupLimiter(); + result += test.timeConstantTags(); + result += test.parameterValidation(); + result += test.signalValidation(); +#ifdef GRIDKIT_ENABLE_ENZYME + result += test.jacobian(); +#endif + + return result.summary(); +}