diff --git a/CHANGELOG.md b/CHANGELOG.md index ba0f1dde6..1028ca570 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -52,6 +52,7 @@ - Added component model developer checklist to a README file. - Added `IEEEST` Stabilizer Model - Added `SEXS-PTI` Exciter Model +- Added `ESDC1A` Exciter Model - Added `GENSAL` Machine Model - Added 200 Bus Synthetic Illinois Case - Added node objects to `PowerElectronics` module & updated all examples to make use of them. 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..b862d41e3 100644 --- a/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp +++ b/GridKit/Model/PhasorDynamics/ComponentLibrary.hpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include diff --git a/GridKit/Model/PhasorDynamics/Exciter/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Exciter/CMakeLists.txt index 1120e20c2..7e5e4b9e3 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/CMakeLists.txt +++ b/GridKit/Model/PhasorDynamics/Exciter/CMakeLists.txt @@ -3,5 +3,6 @@ # - Luke Lowery # ]] +add_subdirectory(ESDC1A) add_subdirectory(IEEET1) add_subdirectory(SEXS-PTI) diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/CMakeLists.txt b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/CMakeLists.txt new file mode 100644 index 000000000..53807c5df --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/CMakeLists.txt @@ -0,0 +1,46 @@ +# [[ +# Author(s): +# - Luke Lowery +# ]] + +set(_install_headers Esdc1a.hpp Esdc1aData.hpp) + +if(GRIDKIT_ENABLE_ENZYME) + gridkit_add_library( + phasor_dynamics_exciter_esdc1a + SOURCES Esdc1aEnzyme.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_exciter_esdc1a + SOURCES Esdc1a.cpp + HEADERS ${_install_headers} + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES GridKit::phasor_dynamics_core GridKit::phasor_dynamics_signal) +endif() + +gridkit_add_library( + phasor_dynamics_exciter_esdc1a_dependency_tracking + SOURCES Esdc1aDependencyTracking.cpp + INCLUDE_DIRECTORIES PRIVATE ${GRIDKIT_THIRD_PARTY_DIR}/magic-enum/include + LINK_LIBRARIES GridKit::phasor_dynamics_core GridKit::phasor_dynamics_signal_dependency_tracking) + +target_link_libraries( + phasor_dynamics_components + INTERFACE GridKit::phasor_dynamics_exciter_esdc1a) +target_link_libraries( + phasor_dynamics_components_dependency_tracking + INTERFACE GridKit::phasor_dynamics_exciter_esdc1a_dependency_tracking) diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.cpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.cpp new file mode 100644 index 000000000..6c5c0f89c --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.cpp @@ -0,0 +1,27 @@ +/** + * @file Esdc1a.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Non-Enzyme instantiation for the ESDC1A exciter model. + */ + +#include "Esdc1aImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Exciter + { + template + int Esdc1a::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Esdc1a..." << std::endl; + Log::misc() << "Jacobian evaluation not implemented!" << std::endl; + return 0; + } + + template class Esdc1a; + template class Esdc1a; + } // namespace Exciter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp new file mode 100644 index 000000000..beeed6538 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1a.hpp @@ -0,0 +1,156 @@ +/** + * @file Esdc1a.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Declaration of the ESDC1A exciter model. + */ + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + template + class BusBase; + + template + class SignalNode; + + namespace Exciter + { + /// Internal variables of an `Esdc1a`. + enum class Esdc1aInternalVariables : size_t + { + EFDP, ///< Field-voltage state before optional speed multiplier + VC, ///< Sensed compensated voltage + VR, ///< Voltage-regulator output + VF, ///< Stabilizing feedback output + XLL, ///< Lead-lag state + EV, ///< Voltage-regulator input error + VLL, ///< Lead-lag block output + VHV, ///< High-value gate output + SE, ///< Saturation coefficient + VFE, ///< Exciter feedback signal + EFD, ///< Field-voltage output + MAXIMUM, + }; + + /// External variables of an `Esdc1a`. + enum class Esdc1aExternalVariables : size_t + { + OMEGA, ///< Machine speed deviation + VS, ///< Stabilizer input signal + VUEL, ///< Under-excitation limiter input + MAXIMUM, + }; + + template + class Esdc1a : public Component + { + using Component::abs_tol_; + using Component::alpha_; + 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::variable_indices_; + using Component::wb_; + using Component::y_; + using Component::yp_; + + public: + using RealT = typename Component::RealT; + using bus_type = BusBase; + using signal_type = SignalNode; + using model_data_type = Esdc1aData; + using MonitorT = Model::VariableMonitor; + + Esdc1a(bus_type* bus); + Esdc1a(bus_type* bus, const model_data_type& data); + ~Esdc1a(); + + 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 rel_tol) 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 model_data_type& data); + void setDerivedParams(); + void initializeMonitor(); + + static constexpr RealT TIME_CONSTANT_MINIMUM = static_cast(1.0e-3); + + bus_type* bus_{nullptr}; + + RealT Tr_{0.0}; + RealT Ka_{40.0}; + RealT Ta_{0.1}; + RealT Tb_{0.0}; + RealT Tc_{0.0}; + RealT Vrmax_{1.0}; + RealT Vrmin_{-1.0}; + RealT Ke_{0.1}; + RealT Te_{0.5}; + RealT Kf_{0.05}; + RealT Tf1_{0.7}; + RealT spdmlt_{0.0}; + RealT E1_{2.8}; + RealT Se1_{0.08}; + RealT E2_{3.7}; + RealT Se2_{0.33}; + IdxT UEL_{0}; + RealT exclim_{1.0}; + + IdxT parameter_error_count_{0}; + + RealT sUEL_{0}; + RealT sUELoff_{1}; + RealT slim_{0}; + RealT slim_off_{1}; + RealT SA_{0}; + RealT SB_{0}; + + ScalarT vref_{0}; + + ComponentSignals signals_; + std::unique_ptr monitor_; + + std::vector ws_; + std::vector ws_indices_; + }; + } // namespace Exciter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aData.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aData.hpp new file mode 100644 index 000000000..0c6a23797 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aData.hpp @@ -0,0 +1,76 @@ +/** + * @file Esdc1aData.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Modeling data for the ESDC1A exciter model. + */ + +#pragma once + +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Exciter + { + /// Parameter keys for the ESDC1A exciter model. + enum class Esdc1aParameters + { + Tr, ///< Transducer time constant + Ka, ///< Voltage-regulator gain + Ta, ///< Voltage-regulator time constant + Tb, ///< Lead-lag denominator time constant + Tc, ///< Lead-lag numerator time constant + Vrmax, ///< Maximum voltage-regulator output + Vrmin, ///< Minimum voltage-regulator output + Ke, ///< Exciter field-resistance line-slope margin + Te, ///< Exciter field time constant + Kf, ///< Stabilizing feedback gain + Tf1, ///< Feedback lead time constant + Spdmlt, ///< Speed multiplier flag + E1, ///< First saturation voltage point + Se1, ///< Saturation value at E1 + E2, ///< Second saturation voltage point + Se2, ///< Saturation value at E2 + UEL, ///< UEL input-location selector + exclim ///< Exciter feedback lower-limit flag + }; + + /// Ports for the ESDC1A exciter model. + enum class Esdc1aPorts + { + bus, ///< Unique ID of the terminal bus + speed, ///< Unique ID of the generator speed-deviation signal + efd, ///< Unique ID of the output EFD signal + vs, ///< Unique ID of the optional stabilizer input signal + vuel ///< Unique ID of the optional UEL input signal + }; + + /// Variables available through the monitor interface. + enum class Esdc1aMonitorableVariables + { + efd, ///< Field-voltage output + vc, ///< Sensed compensated voltage + vr, ///< Voltage-regulator output + vf, ///< Stabilizing feedback state + se, ///< Saturation coefficient + vfe ///< Exciter feedback signal + }; + + template + struct Esdc1aData : public ComponentData + { + Esdc1aData() = default; + + using Parameters = Esdc1aParameters; + using Ports = Esdc1aPorts; + using MonitorableVariables = Esdc1aMonitorableVariables; + }; + } // namespace Exciter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aDependencyTracking.cpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aDependencyTracking.cpp new file mode 100644 index 000000000..c7c25aec6 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aDependencyTracking.cpp @@ -0,0 +1,27 @@ +/** + * @file Esdc1aDependencyTracking.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Dependency-tracking instantiations for the ESDC1A exciter model. + */ + +#include "Esdc1aImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Exciter + { + template + int Esdc1a::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Esdc1a..." << std::endl; + Log::misc() << "Jacobian evaluation not implemented!" << std::endl; + return 0; + } + + template class Esdc1a; + template class Esdc1a; + } // namespace Exciter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aEnzyme.cpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aEnzyme.cpp new file mode 100644 index 000000000..f74326132 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aEnzyme.cpp @@ -0,0 +1,105 @@ +/** + * @file Esdc1aEnzyme.cpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Enzyme sparse Jacobian for the ESDC1A exciter model. + */ + +#include + +#include "Esdc1aImpl.hpp" + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Exciter + { + template + int Esdc1a::evaluateJacobian() + { + Log::misc() << "Evaluate Jacobian for Esdc1a..." << std::endl; + Log::misc() << "Jacobian evaluation is experimental!" << std::endl; + + if (J_rows_buffer_ == nullptr) + { + auto size = static_cast(size_); + auto bus_size = static_cast(bus_->size()); + auto signal_size = static_cast(ws_.size()); + auto buffer_size = 2 * size * size + size * bus_size + size * signal_size; + J_rows_buffer_ = new IdxT[buffer_size]; + J_cols_buffer_ = new IdxT[buffer_size]; + J_vals_buffer_ = new RealT[buffer_size]; + } + + using ModelT = GridKit::PhasorDynamics::Exciter::Esdc1a; + using Fn = GridKit::Enzyme::Sparse::MemberFunctions; + + nnz_ = 0; + + GridKit::Enzyme::Sparse::DfDy::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + GridKit::Enzyme::Sparse::DfDyp::eval(this, + f_.size(), + y_.size(), + (this->getResidualIndices()).data(), + (this->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + alpha_, + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + GridKit::Enzyme::Sparse::DfDwb::eval(this, + f_.size(), + static_cast(bus_->size()), + (this->getResidualIndices()).data(), + (bus_->getVariableIndices()).data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + GridKit::Enzyme::Sparse::DfDws::eval(this, + f_.size(), + ws_.size(), + (this->getResidualIndices()).data(), + ws_indices_.data(), + y_.data(), + yp_.data(), + wb_.data(), + ws_.data(), + J_rows_buffer_, + J_cols_buffer_, + J_vals_buffer_, + nnz_); + + this->constructCoo(); + + return 0; + } + + template class Esdc1a; + template class Esdc1a; + } // namespace Exciter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp new file mode 100644 index 000000000..10962f251 --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/Esdc1aImpl.hpp @@ -0,0 +1,568 @@ +/** + * @file Esdc1aImpl.hpp + * @author Luke Lowery (lukel@tamu.edu) + * @brief Definition of the ESDC1A exciter model. + */ + +#pragma once + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace PhasorDynamics + { + namespace Exciter + { + using Log = ::GridKit::Utilities::Logger; + + template + Esdc1a::Esdc1a(bus_type* bus) + : bus_(bus) + { + setDerivedParams(); + size_ = static_cast(Esdc1aInternalVariables::MAXIMUM); + } + + template + Esdc1a::Esdc1a(bus_type* bus, const model_data_type& data) + : bus_(bus), + monitor_(std::make_unique(data)) + { + initModelParams(data); + setDerivedParams(); + initializeMonitor(); + size_ = static_cast(Esdc1aInternalVariables::MAXIMUM); + } + + template + Esdc1a::~Esdc1a() + { + } + + template + void Esdc1a::initModelParams(const model_data_type& data) + { + using Params = typename model_data_type::Parameters; + + parameter_error_count_ = 0; + + auto load_real = [&](auto key, RealT& target, const char* name) + { + if (!data.parameters.contains(key)) + { + return; + } + + const auto& value = data.parameters.at(key); + if (const auto* real_value = std::get_if(&value)) + { + target = *real_value; + } + else if (const auto* index_value = std::get_if(&value)) + { + target = static_cast(*index_value); + } + else + { + Log::error() << "Esdc1a: parameter '" << name << "' must be numeric\n"; + ++parameter_error_count_; + } + }; + + auto load_switch = [&](auto key, RealT& target, const char* name) + { + if (!data.parameters.contains(key)) + { + return; + } + + const auto& value = data.parameters.at(key); + if (const auto* bool_value = std::get_if(&value)) + { + target = *bool_value ? ONE : ZERO; + } + else if (const auto* index_value = std::get_if(&value); + index_value && (*index_value == 0 || *index_value == 1)) + { + target = static_cast(*index_value); + } + else if (const auto* real_value = std::get_if(&value); + real_value && (*real_value == ZERO || *real_value == ONE) ) + { + target = *real_value; + } + else + { + Log::error() << "Esdc1a: parameter '" << name << "' must be bool or 0/1\n"; + ++parameter_error_count_; + } + }; + + auto load_selector = [&](auto key, IdxT& target, const char* name) + { + if (!data.parameters.contains(key)) + { + return; + } + + const auto& value = data.parameters.at(key); + if (const auto* index_value = std::get_if(&value)) + { + target = *index_value; + } + else if (const auto* real_value = std::get_if(&value)) + { + const RealT rounded = std::round(*real_value); + if (*real_value >= ZERO && *real_value == rounded) + { + target = static_cast(rounded); + } + else + { + Log::error() << "Esdc1a: parameter '" << name << "' must be an integer selector\n"; + ++parameter_error_count_; + } + } + else + { + Log::error() << "Esdc1a: parameter '" << name << "' must be an integer selector\n"; + ++parameter_error_count_; + } + }; + + load_real(Params::Tr, Tr_, "Tr"); + load_real(Params::Ka, Ka_, "Ka"); + load_real(Params::Ta, Ta_, "Ta"); + load_real(Params::Tb, Tb_, "Tb"); + load_real(Params::Tc, Tc_, "Tc"); + load_real(Params::Vrmax, Vrmax_, "Vrmax"); + load_real(Params::Vrmin, Vrmin_, "Vrmin"); + load_real(Params::Ke, Ke_, "Ke"); + load_real(Params::Te, Te_, "Te"); + load_real(Params::Kf, Kf_, "Kf"); + load_real(Params::Tf1, Tf1_, "Tf1"); + load_switch(Params::Spdmlt, spdmlt_, "Spdmlt"); + load_real(Params::E1, E1_, "E1"); + load_real(Params::Se1, Se1_, "Se1"); + load_real(Params::E2, E2_, "E2"); + load_real(Params::Se2, Se2_, "Se2"); + load_selector(Params::UEL, UEL_, "UEL"); + load_switch(Params::exclim, exclim_, "exclim"); + } + + template + void Esdc1a::setDerivedParams() + { + Tr_ = std::max(Tr_, TIME_CONSTANT_MINIMUM); + Tb_ = std::max(Tb_, TIME_CONSTANT_MINIMUM); + Tf1_ = std::max(Tf1_, TIME_CONSTANT_MINIMUM); + + sUEL_ = UEL_ >= static_cast(2) ? ONE : ZERO; + sUELoff_ = ONE - sUEL_; + slim_ = exclim_; + slim_off_ = ONE - slim_; + + if (Se1_ == ZERO && Se2_ == ZERO) + { + SA_ = ZERO; + SB_ = ZERO; + return; + } + if (E1_ <= ZERO || E2_ <= ZERO || E1_ == E2_ + || Se1_ <= ZERO || Se2_ <= ZERO || Se1_ == Se2_) + { + SA_ = ZERO; + SB_ = ZERO; + return; + } + + const RealT C = std::sqrt(Se2_ / Se1_); + SA_ = (C * E1_ - E2_) / (C - ONE); + SB_ = Se1_ / ((E1_ - SA_) * (E1_ - SA_)); + } + + template + const Model::VariableMonitorBase* Esdc1a::getMonitor() const + { + return monitor_.get(); + } + + template + void Esdc1a::initializeMonitor() + { + using Variable = typename model_data_type::MonitorableVariables; + auto index = [](Esdc1aInternalVariables variable) + { + return static_cast(variable); + }; + + monitor_->set(Variable::efd, [this, index] + { return y_[index(Esdc1aInternalVariables::EFD)]; }); + monitor_->set(Variable::vc, [this, index] + { return y_[index(Esdc1aInternalVariables::VC)]; }); + monitor_->set(Variable::vr, [this, index] + { return y_[index(Esdc1aInternalVariables::VR)]; }); + monitor_->set(Variable::vf, [this, index] + { return y_[index(Esdc1aInternalVariables::VF)]; }); + monitor_->set(Variable::se, [this, index] + { return y_[index(Esdc1aInternalVariables::SE)]; }); + monitor_->set(Variable::vfe, [this, index] + { return y_[index(Esdc1aInternalVariables::VFE)]; }); + } + + template + int Esdc1a::setGridKitComponentID(IdxT component_id) + { + gridkit_component_id_ = component_id; + return 0; + } + + template + int Esdc1a::allocate() + { + size_ = static_cast(Esdc1aInternalVariables::MAXIMUM); + auto size = static_cast(size_); + + f_.assign(size, ScalarT{0}); + y_.assign(size, ScalarT{0}); + yp_.assign(size, ScalarT{0}); + tag_.assign(size, false); + abs_tol_.assign(size, ScalarT{0}); + variable_indices_.resize(size); + residual_indices_.resize(size); + + wb_.assign(2, ScalarT{0}); + + auto signal_size = static_cast(Esdc1aExternalVariables::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(Esdc1aInternalVariables::EFD)], + &(this->getVariableIndex(static_cast(Esdc1aInternalVariables::EFD)))); + } + + return 0; + } + + template + int Esdc1a::verify() const + { + int ret = static_cast(parameter_error_count_); + + auto check = [&](bool condition, const char* message) + { + if (!condition) + { + Log::error() << "Esdc1a: " << message << '\n'; + ret += 1; + } + }; + + if (bus_ == nullptr) + { + Log::error() << "Esdc1a: bus pointer is null\n"; + ret += 1; + } + + check(Ka_ > ZERO, "Ka must be positive"); + check(Ta_ > ZERO, "Ta must be positive"); + check(Tc_ >= ZERO, "Tc must be non-negative"); + check(Te_ > ZERO, "Te must be positive"); + check(Vrmin_ <= Vrmax_, "Vrmin must be less than or equal to Vrmax"); + check(spdmlt_ == ZERO || spdmlt_ == ONE, "Spdmlt must be 0 or 1"); + check(exclim_ == ZERO || exclim_ == ONE, "exclim must be 0 or 1"); + check(UEL_ >= static_cast(0) && UEL_ <= static_cast(3), + "UEL must be 0, 1, 2, or 3"); + + if (!(Se1_ == ZERO && Se2_ == ZERO) ) + { + check(E1_ > ZERO, "E1 must be positive when saturation is enabled"); + check(E2_ > ZERO, "E2 must be positive when saturation is enabled"); + check(Se1_ > ZERO, "Se1 must be positive when saturation is enabled"); + check(Se2_ > ZERO, "Se2 must be positive when saturation is enabled"); + check(E1_ != E2_, "E1 and E2 must differ when saturation is enabled"); + check(Se1_ != Se2_, "Se1 and Se2 must differ when saturation is enabled"); + } + + if (!signals_.template isAssigned()) + { + Log::error() << "Esdc1a: required EFD signal is not assigned\n"; + ret += 1; + } + + if (spdmlt_ == ONE + && !signals_.template isAttached()) + { + Log::error() << "Esdc1a: speed signal is required when Spdmlt is enabled\n"; + ret += 1; + } + + auto check_attached_signal = [&](bool attached, bool linked, const char* name) + { + if (attached && !linked) + { + Log::error() << "Esdc1a: " << name << " signal attached with no linked source\n"; + ret += 1; + } + }; + + check_attached_signal( + signals_.template isAttached(), + signals_.template isAttached() + && signals_.template isLinked(), + "speed"); + check_attached_signal( + signals_.template isAttached(), + signals_.template isAttached() + && signals_.template isLinked(), + "VS"); + check_attached_signal( + signals_.template isAttached(), + signals_.template isAttached() + && signals_.template isLinked(), + "VUEL"); + + return ret; + } + + template + int Esdc1a::initialize() + { + if (verify() > 0) + { + Log::error() << "Esdc1a: cannot initialize with invalid configuration\n"; + return 1; + } + + const auto EFDP = static_cast(Esdc1aInternalVariables::EFDP); + const auto VC = static_cast(Esdc1aInternalVariables::VC); + const auto VR = static_cast(Esdc1aInternalVariables::VR); + const auto VF = static_cast(Esdc1aInternalVariables::VF); + const auto XLL = static_cast(Esdc1aInternalVariables::XLL); + const auto EV = static_cast(Esdc1aInternalVariables::EV); + const auto VLL = static_cast(Esdc1aInternalVariables::VLL); + const auto VHV = static_cast(Esdc1aInternalVariables::VHV); + const auto SE = static_cast(Esdc1aInternalVariables::SE); + const auto VFE = static_cast(Esdc1aInternalVariables::VFE); + const auto EFD = static_cast(Esdc1aInternalVariables::EFD); + + ScalarT omega0{ZERO}; + if (signals_.template isAttached()) + { + omega0 = signals_.template readExternalVariable(); + } + + ScalarT vs0{ZERO}; + if (signals_.template isAttached()) + { + vs0 = signals_.template readExternalVariable(); + } + + ScalarT vuel0{ZERO}; + if (signals_.template isAttached()) + { + vuel0 = signals_.template readExternalVariable(); + } + + const ScalarT d0 = ONE + spdmlt_ * omega0; + if (d0 == ZERO) + { + Log::error() << "Esdc1a: speed multiplier denominator is zero at initialization\n"; + return 1; + } + + const ScalarT Ec0 = std::sqrt(bus_->Vr() * bus_->Vr() + bus_->Vi() * bus_->Vi()); + + const ScalarT efd0 = y_[EFD]; + const ScalarT efdp0 = efd0 / d0; + const ScalarT se0 = SB_ * Math::qramp(efdp0 - SA_); + const ScalarT vfe0 = slim_off_ * (Ke_ + se0) * efdp0 + + slim_ * Math::ramp((Ke_ + se0) * efdp0); + const ScalarT vr0 = vfe0; + const ScalarT vhv0 = vr0 / Ka_; + auto inverse_ramp = [](RealT y) + { + const RealT scaled_y = Math::MU * y; + if (scaled_y > static_cast(50.0)) + { + return y; + } + return std::log(std::expm1(scaled_y)) / Math::MU; + }; + + ScalarT gate_input0 = vhv0; + if (sUEL_ == ZERO) + { + const RealT ramp_target = static_cast(vhv0 - vuel0); + if (ramp_target <= ZERO) + { + Log::error() << "Esdc1a: smooth high-value gate is active at initialization\n"; + return 1; + } + gate_input0 = vuel0 + inverse_ramp(ramp_target); + } + + const ScalarT vc0 = Ec0; + const ScalarT vf0 = ScalarT{ZERO}; + const ScalarT ev0 = gate_input0; + const ScalarT xll0 = gate_input0; + const ScalarT vll0 = gate_input0; + + if (vr0 < Vrmin_ || vr0 > Vrmax_) + { + Log::error() << "Esdc1a: initialized VR is outside limits\n"; + return 1; + } + + vref_ = ev0 + vc0 + vf0 - vs0 - sUEL_ * vuel0; + + y_[EFDP] = efdp0; + y_[VC] = vc0; + y_[VR] = vr0; + y_[VF] = vf0; + y_[XLL] = xll0; + y_[EV] = ev0; + y_[VLL] = vll0; + y_[VHV] = vhv0; + y_[SE] = se0; + y_[VFE] = vfe0; + y_[EFD] = efd0; + + std::fill(yp_.begin(), yp_.end(), ZERO); + return 0; + } + + template + int Esdc1a::tagDifferentiable() + { + std::fill(tag_.begin(), tag_.end(), false); + tag_[static_cast(Esdc1aInternalVariables::EFDP)] = true; + tag_[static_cast(Esdc1aInternalVariables::VC)] = true; + tag_[static_cast(Esdc1aInternalVariables::VR)] = true; + tag_[static_cast(Esdc1aInternalVariables::VF)] = true; + tag_[static_cast(Esdc1aInternalVariables::XLL)] = true; + return 0; + } + + template + int Esdc1a::setAbsoluteTolerance(RealT rel_tol) + { + std::fill(abs_tol_.begin(), abs_tol_.end(), rel_tol); + return 0; + } + + template + __attribute__((always_inline)) inline int Esdc1a::evaluateInternalResidual( + const ScalarT* y, + const ScalarT* yp, + const ScalarT* wb, + const ScalarT* ws, + ScalarT* f) + { + const auto EFDP = static_cast(Esdc1aInternalVariables::EFDP); + const auto VC = static_cast(Esdc1aInternalVariables::VC); + const auto VR = static_cast(Esdc1aInternalVariables::VR); + const auto VF = static_cast(Esdc1aInternalVariables::VF); + const auto XLL = static_cast(Esdc1aInternalVariables::XLL); + const auto EV = static_cast(Esdc1aInternalVariables::EV); + const auto VLL = static_cast(Esdc1aInternalVariables::VLL); + const auto VHV = static_cast(Esdc1aInternalVariables::VHV); + const auto SE = static_cast(Esdc1aInternalVariables::SE); + const auto VFE = static_cast(Esdc1aInternalVariables::VFE); + const auto EFD = static_cast(Esdc1aInternalVariables::EFD); + + const auto OMEGA = static_cast(Esdc1aExternalVariables::OMEGA); + const auto VS = static_cast(Esdc1aExternalVariables::VS); + const auto VUEL = static_cast(Esdc1aExternalVariables::VUEL); + + const ScalarT efdp = y[EFDP]; + const ScalarT vc = y[VC]; + const ScalarT vr = y[VR]; + const ScalarT vf = y[VF]; + const ScalarT xll = y[XLL]; + const ScalarT ev = y[EV]; + const ScalarT vll = y[VLL]; + const ScalarT vhv = y[VHV]; + const ScalarT se = y[SE]; + const ScalarT vfe = y[VFE]; + const ScalarT efd = y[EFD]; + + const ScalarT omega = ws[OMEGA]; + const ScalarT vs = ws[VS]; + const ScalarT vuel = ws[VUEL]; + + const ScalarT Ec = std::sqrt(wb[0] * wb[0] + wb[1] * wb[1]); + const ScalarT ev_target = vref_ + vs + sUEL_ * vuel - vc - vf; + + f[EFDP] = -yp[EFDP] + (vr - vfe) / Te_; + f[VC] = -yp[VC] + (Ec - vc) / Tr_; + f[VR] = -yp[VR] + Math::antiwindup(vr, -vr + Ka_ * vhv, Vrmin_, Vrmax_) / Ta_; + f[VF] = -yp[VF] + (-vf + Kf_ * (vr - vfe) / Te_) / Tf1_; + f[XLL] = -yp[XLL] + (ev - xll) / Tb_; + f[EV] = -ev + ev_target; + f[VLL] = -vll + xll + (Tc_ / Tb_) * (ev - xll); + f[VHV] = -vhv + sUEL_ * vll + sUELoff_ * Math::max(vll, vuel); + f[SE] = -se + SB_ * Math::qramp(efdp - SA_); + f[VFE] = -vfe + slim_off_ * (Ke_ + se) * efdp + slim_ * Math::ramp((Ke_ + se) * efdp); + f[EFD] = -efd + (ONE + spdmlt_ * omega) * efdp; + + return 0; + } + + template + int Esdc1a::evaluateResidual() + { + const auto OMEGA = static_cast(Esdc1aExternalVariables::OMEGA); + const auto VS = static_cast(Esdc1aExternalVariables::VS); + const auto VUEL = static_cast(Esdc1aExternalVariables::VUEL); + + std::fill(ws_.begin(), ws_.end(), ScalarT{ZERO}); + 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_[VS] = signals_.template readExternalVariable(); + ws_indices_[VS] = + signals_.template readExternalVariableIndex(); + } + if (signals_.template isAttached()) + { + ws_[VUEL] = signals_.template readExternalVariable(); + ws_indices_[VUEL] = + signals_.template readExternalVariableIndex(); + } + + wb_[0] = bus_->Vr(); + wb_[1] = bus_->Vi(); + + evaluateInternalResidual(y_.data(), yp_.data(), wb_.data(), ws_.data(), f_.data()); + return 0; + } + } // namespace Exciter + } // namespace PhasorDynamics +} // namespace GridKit diff --git a/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md new file mode 100644 index 000000000..5bc7ea52e --- /dev/null +++ b/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md @@ -0,0 +1,329 @@ +# **IEEE DC1A Excitation System Model (ESDC1A)** + +ESDC1A is an IEEE DC1A excitation-system model. In GridKit it reads the +connected bus voltage, optional stabilizer and under-excitation limiter signals, +and publishes field voltage. + +## Notes + +- Internal voltage signals are on model base unless otherwise stated. +- The connected bus supplies $E_C=\sqrt{V_{\mathrm{r}}^2+V_{\mathrm{i}}^2}$. +- The source diagram labels the optional multiplier input as `Speed`; GridKit + uses machine speed deviation, so the enabled multiplier is $1+\omega$. +- The PowerWorld parameter table names the UEL selector `UEL`; `UEL >= 2` + routes the UEL input through the input-error summing junction and `UEL < 2` + routes it through the high-value gate. +- `efd` is a required output signal. `speed` is required only when + $s_{\mathrm{spd}}=1$; `vs` and `vuel` are optional and default to zero. + +## Block Diagram + +Standard ESDC1A block diagram. + +![](../../../../../docs/Figures/PhasorDynamics/ESDC1A_diagram.png) + +Figure 1: ESDC1A block diagram. Figure courtesy of [PowerWorld](https://www.powerworld.com/WebHelp/) + +## Model Parameters + +Symbol | Units | JSON | Description | Typical Value | Note +------------------------------------|-----------|-----------|--------------------------------------------------|---------------|------ +$T_R$ | [sec] | `Tr` | Transducer time constant | 0.0 | State 2 +$K_A$ | [p.u.] | `Ka` | Voltage-regulator gain | 40.0 | +$T_A$ | [sec] | `Ta` | Voltage-regulator time constant | 0.1 | State 3 +$T_B$ | [sec] | `Tb` | Lead-lag denominator time constant | 0.0 | State 5 +$T_C$ | [sec] | `Tc` | Lead-lag numerator time constant | 0.0 | +$V_R^{\max}$ | [p.u.] | `Vrmax` | Maximum voltage-regulator output | 1.0 | +$V_R^{\min}$ | [p.u.] | `Vrmin` | Minimum voltage-regulator output | -1.0 | +$K_E$ | [p.u.] | `Ke` | Exciter field-resistance line-slope margin | 0.1 | +$T_E$ | [sec] | `Te` | Exciter field time constant | 0.5 | State 1 +$K_F$ | [p.u.] | `Kf` | Stabilizing feedback gain | 0.05 | +$T_{F1}$ | [sec] | `Tf1` | Feedback lead time constant | 0.7 | State 4 +$s_{\mathrm{spd}}$ | [binary] | `Spdmlt` | Speed multiplier flag | 0.0 | 1 enables the speed multiplier +$E_1$ | [p.u.] | `E1` | First saturation voltage point | 2.8 | +$S_E(E_1)$ | [p.u.] | `Se1` | Saturation value at $E_1$ | 0.08 | +$E_2$ | [p.u.] | `E2` | Second saturation voltage point | 3.7 | +$S_E(E_2)$ | [p.u.] | `Se2` | Saturation value at $E_2$ | 0.33 | +$I_{\mathrm{UEL}}$ | [integer] | `UEL` | Under-excitation limiter input-location selector | 0 | 0/1 = high-value gate, 2/3 = input-error summing junction +$s_{\mathrm{lim}}$ | [binary] | `exclim` | Exciter feedback lower-limit flag | 1.0 | 1 enables the zero lower limit on $V_{\mathrm{FE}}$ + +### Parameter Validation + +Invalid ESDC1A parameter sets are rejected by the following checks. Let $\epsilon_T=10^{-3}$. + +```math +\begin{aligned} + T &\leftarrow \max\!\left(T, \epsilon_T\right) + \quad T\in\{T_R,T_B,T_{F1}\} \\ + K_A + &> 0 \\ + T_A, T_E + &> 0 \\ + T_C + &\ge 0 \\ + V_R^{\min} + &\le V_R^{\max} \\ + s_{\mathrm{spd}}, s_{\mathrm{lim}} + &\in \{0,1\} \\ + I_{\mathrm{UEL}} + &\in \{0,1,2,3\} \\ + \left(S_E(E_1), S_E(E_2)\right) + &=(0,0) + \quad\text{or}\quad + \begin{gathered} + E_1, E_2, S_E(E_1), S_E(E_2) > 0 \\ + E_1 \ne E_2 \\ + S_E(E_1) \ne S_E(E_2) + \end{gathered} +\end{aligned} +``` + +### Model Derived Parameters + +```math +\begin{aligned} + s_{\mathrm{UEL}} + &= + \begin{cases} + 1 & I_{\mathrm{UEL}} \ge 2 \\ + 0 & I_{\mathrm{UEL}} < 2 + \end{cases} \\ + s_{\mathrm{UEL}}^\mathrm{off} + &= 1 - s_{\mathrm{UEL}} \\ + s_{\mathrm{lim}}^\mathrm{off} + &= 1 - s_{\mathrm{lim}} +\end{aligned} +``` + +When saturation is disabled, $S_A=0$ and $S_B=0$. Otherwise, + +```math +\begin{aligned} + C &= \sqrt{\dfrac{S_E(E_2)}{S_E(E_1)}} \\ + S_A &= \dfrac{C E_1 - E_2}{C - 1} \\ + S_B &= \dfrac{S_E(E_1)}{(E_1 - S_A)^2} +\end{aligned} +``` + +## Model Variables + +### Internal Variables + +#### Differential + +Symbol | Units | Description | Note +------------------------------------|--------|------------------------------------------------------|------ +$E_{\mathrm{fd}}'$ | [p.u.] | Field-voltage state before optional speed multiplier | State 1 in Fig. 1; source label: `EFD` +$V_C$ | [p.u.] | Sensed compensated voltage | State 2 in Fig. 1; source label: `Sensed Vt` +$V_R$ | [p.u.] | Voltage-regulator output | State 3 in Fig. 1; source label: `VR` +$V_F$ | [p.u.] | Stabilizing feedback output | State 4 in Fig. 1; source label: `VF` +$x_{\mathrm{LL}}$ | [p.u.] | Lead-lag block state | State 5 in Fig. 1; source label: `Lead-Lag` + +#### Algebraic + +Symbol | Units | Description | Note +------------------------------------|--------|--------------------------------------------------|------ +$e_V$ | [p.u.] | Voltage-regulator input error | +$V_{\mathrm{LL}}$ | [p.u.] | Lead-lag block output | Input to high-value gate +$V_{\mathrm{HV}}$ | [p.u.] | High-value gate output | Input to voltage regulator +$S_E$ | [p.u.] | Saturation coefficient evaluated at $E_{\mathrm{fd}}'$ | +$V_{\mathrm{FE}}$ | [p.u.] | Exciter feedback signal after optional lower limit | +$E_{\mathrm{fd}}$ | [p.u.] | Field-voltage output | Published through `efd` + +### External Variables + +#### Differential + +None. + +#### Algebraic + +Symbol | Units | Type | Description | Note +------------------------------------|--------|---------|-------------------------------------|------ +$V_{\mathrm{r}}$ | [p.u.] | Known | Terminal-bus voltage, real component | Bus input +$V_{\mathrm{i}}$ | [p.u.] | Known | Terminal-bus voltage, imaginary component | Bus input +$\omega$ | [p.u.] | Known | Machine speed deviation | Optional signal port `speed`; required when $s_{\mathrm{spd}}=1$ +$V_S$ | [p.u.] | Known | Stabilizer input signal | Optional signal port `vs`; defaults to zero +$V_{\mathrm{UEL}}$ | [p.u.] | Known | Under-excitation limiter input | Optional signal port `vuel`; defaults to zero +$V_{\mathrm{ref}}$ | [p.u.] | Unknown | Voltage-control reference | Initialized constant setpoint; source label: `VREF` + +## Model Equations + +### Differential Equations + +```math +\begin{aligned} + 0 &= + -\dot{E}_{\mathrm{fd}}' + + \dfrac{1}{T_E} + \left(V_R - V_{\mathrm{FE}}\right) \\ + 0 &= + -\dot{V}_C + + \dfrac{1}{T_R} + \left( + \sqrt{V_{\mathrm{r}}^2+V_{\mathrm{i}}^2} + - V_C + \right) \\ + 0 &= + -\dot{V}_R + + \dfrac{1}{T_A} + \text{antiwindup} + \left( + V_R,\, + -V_R + K_A V_{\mathrm{HV}};\, + V_R^{\min}, V_R^{\max} + \right) \\ + 0 &= + -\dot{V}_F + + \dfrac{1}{T_{F1}} + \left[ + -V_F + + \dfrac{K_F}{T_E} + \left(V_R - V_{\mathrm{FE}}\right) + \right] \\ + 0 &= + -\dot{x}_{\mathrm{LL}} + + \dfrac{1}{T_B} + \left(e_V - x_{\mathrm{LL}}\right) +\end{aligned} +``` + +CommonMath defines the [Anti-Windup](../../../../CommonMath.md#anti-windup-indicator) +target and smooth approximation. + +### Algebraic Equations + +```math +\begin{aligned} + 0 &= + -e_V + + V_{\mathrm{ref}} + + V_S + + s_{\mathrm{UEL}}V_{\mathrm{UEL}} + - V_C + - V_F \\ + 0 &= + -V_{\mathrm{LL}} + + x_{\mathrm{LL}} + + \dfrac{T_C}{T_B} + \left(e_V - x_{\mathrm{LL}}\right) \\ + 0 &= + -V_{\mathrm{HV}} + + s_{\mathrm{UEL}}V_{\mathrm{LL}} + + s_{\mathrm{UEL}}^\mathrm{off} + \text{max}\left(V_{\mathrm{LL}}, V_{\mathrm{UEL}}\right) \\ + 0 &= + -S_E + + S_B q\left(E_{\mathrm{fd}}' - S_A\right) \\ + 0 &= + -V_{\mathrm{FE}} + + s_{\mathrm{lim}}^\mathrm{off} + \left(K_E + S_E\right)E_{\mathrm{fd}}' + + s_{\mathrm{lim}}\rho + \left(\left(K_E + S_E\right)E_{\mathrm{fd}}'\right) \\ + 0 &= + -E_{\mathrm{fd}} + + \left(1 + s_{\mathrm{spd}}\omega\right)E_{\mathrm{fd}}' +\end{aligned} +``` + +CommonMath defines helper targets and smooth approximations for +[max](../../../../CommonMath.md#derived-functions), the [ramp](../../../../CommonMath.md#primitives) +$\rho$, and the [quadratic ramp](../../../../CommonMath.md#primitives) $q$. + +## Initialization + +### External Priors + +```math +\begin{aligned} + V_{\mathrm{r},0}, V_{\mathrm{i},0} + &\leftarrow \text{terminal-bus voltage} \\ + E_{\mathrm{fd},0} + &\leftarrow \text{field-voltage signal start} \\ + \omega_0 + &\leftarrow \text{speed-deviation input or }0 \\ + V_{S,0} + &\leftarrow \text{stabilizer input or }0 \\ + V_{\mathrm{UEL},0} + &\leftarrow \text{under-excitation limiter input or }0 +\end{aligned} +``` + +### Internal 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} + E_{C,0} + &= \sqrt{V_{\mathrm{r},0}^2+V_{\mathrm{i},0}^2} \\ + d_0 + &= 1 + s_{\mathrm{spd}}\omega_0 \\ + E_{\mathrm{fd},0}' + &= \dfrac{E_{\mathrm{fd},0}}{d_0} \\ + S_{E,0} + &= S_B q\left(E_{\mathrm{fd},0}' - S_A\right) \\ + V_{\mathrm{FE},0} + &= + s_{\mathrm{lim}}^\mathrm{off} + \left(K_E + S_{E,0}\right)E_{\mathrm{fd},0}' + + s_{\mathrm{lim}}\rho + \left(\left(K_E + S_{E,0}\right)E_{\mathrm{fd},0}'\right) \\ + V_{R,0} + &= V_{\mathrm{FE},0} \\ + V_{\mathrm{HV},0} + &= \dfrac{V_{R,0}}{K_A} \\ + g_0 + &= + \begin{cases} + V_{\mathrm{HV},0}, + & s_{\mathrm{UEL}}=1 \\ + V_{\mathrm{UEL},0} + + \rho^{-1} + \left(V_{\mathrm{HV},0}-V_{\mathrm{UEL},0}\right), + & s_{\mathrm{UEL}}=0 + \end{cases} \\ + V_{C,0} + &= E_{C,0} \\ + V_{F,0} + &= 0 \\ + e_{V,0} + &= g_0 \\ + x_{\mathrm{LL},0} + &= g_0 \\ + V_{\mathrm{LL},0} + &= g_0 +\end{aligned} +``` + +Initialization rejects $d_0=0$, $V_{R,0}$ outside +$[V_R^{\min},V_R^{\max}]$, and high-value-gate active starts with +$s_{\mathrm{UEL}}=0$ and $V_{\mathrm{HV},0}\le V_{\mathrm{UEL},0}$. + +### External Solved + +```math +\begin{aligned} + V_{\mathrm{ref},0} + &\leftarrow + e_{V,0} + + V_{C,0} + + V_{F,0} + - V_{S,0} + - s_{\mathrm{UEL}}V_{\mathrm{UEL},0} +\end{aligned} +``` + +## Monitorable Outputs + +Output | Units | Description | Note +----------------|--------|-------------------------------------|------ +`efd` | [p.u.] | Field-voltage output | $E_{\mathrm{fd}}$ +`vc` | [p.u.] | Sensed compensated voltage | $V_C$ +`vr` | [p.u.] | Voltage-regulator output | $V_R$ +`vf` | [p.u.] | Stabilizing feedback output | $V_F$ +`se` | [p.u.] | Saturation coefficient | $S_E$ +`vfe` | [p.u.] | Exciter feedback signal | $V_{\mathrm{FE}}$ diff --git a/GridKit/Model/PhasorDynamics/Exciter/README.md b/GridKit/Model/PhasorDynamics/Exciter/README.md index 5391918f0..63ea3bd87 100644 --- a/GridKit/Model/PhasorDynamics/Exciter/README.md +++ b/GridKit/Model/PhasorDynamics/Exciter/README.md @@ -1,7 +1,7 @@ # **Exciter Models** > [!NOTE] -> IEEET1 and SEXS-PTI exciters are currently implemented. +> EXDC1 is not currently implemented. ## Introduction @@ -15,6 +15,7 @@ There are a few standard Exciter models - ESAC6A Excitation Model (See [ESAC6A](ESAC6A/README.md)) - IEEE Type 1 Excitation Model (See [IEEET1](IEEET1/README.md)) - IEEE DC1 Excitation Model (See [EXDC1](EXDC1/README.md)) +- IEEE DC1A Excitation Model (See [ESDC1A](ESDC1A/README.md)) - ESDC2A Excitation Model (See [ESDC2A](ESDC2A/README.md)) - EXAC1 Excitation Model (See [EXAC1](EXAC1/README.md)) - IEEE ST4B Excitation Model (See [ESST4B](ESST4B/README.md)) diff --git a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md index 5bff41ad4..ce4f4038a 100644 --- a/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md +++ b/GridKit/Model/PhasorDynamics/INPUT_FORMAT.md @@ -146,6 +146,7 @@ are specified: `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` `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` + `Esdc1a` | the ESDC1A exciter model | `bus`, `efd`, `speed`\*, `vs`\*, `vuel`\* | `Tr`, `Ka`, `Ta`, `Tb`, `Tc`, `Vrmax`, `Vrmin`, `Ke`, `Te`, `Kf`, `Tf1`, `Spdmlt`, `E1`, `Se1`, `E2`, `Se2`, `UEL`, `exclim` | `efd`, `vc`, `vr`, `vf`, `se`, `vfe` `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` `BusFault` | simple impedance-based fault at a bus | `bus`, `status`\* | `state0`, `R`, `X` | `state`, `ir`, `ii` diff --git a/GridKit/Model/PhasorDynamics/SystemModelData.hpp b/GridKit/Model/PhasorDynamics/SystemModelData.hpp index 8a6e73056..037eebefc 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelData.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelData.hpp @@ -10,6 +10,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 Esdc1aDataT = Exciter::Esdc1aData; using Ieeet1DataT = Exciter::Ieeet1Data; using SexsPtiDataT = Exciter::SexsPtiData; using IeeestDataT = Stabilizer::IeeestData; @@ -101,6 +103,7 @@ namespace GridKit 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 esdc1a; ///< ESDC1A exciters 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 a3682feca..66cd5aa6f 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelDataJSONParser.hpp @@ -148,6 +148,12 @@ namespace GridKit raw_component.get_to(exciter); sm.exciter.push_back(exciter); } + else if (kind == "Esdc1a") + { + typename SystemModelData::Esdc1aDataT exciter; + raw_component.get_to(exciter); + sm.esdc1a.push_back(exciter); + } else if (kind == "SexsPti") { typename SystemModelData::SexsPtiDataT exciter; diff --git a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp index 287b328dc..efe16e812 100644 --- a/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp +++ b/GridKit/Model/PhasorDynamics/SystemModelImpl.hpp @@ -293,6 +293,47 @@ namespace GridKit addComponent(exciter); } + for (const auto& excitedata : data.esdc1a) + { + IdxT bus_index = 0; + if (excitedata.ports.contains(Esdc1aPorts::bus)) + { + bus_index = excitedata.ports.at(Esdc1aPorts::bus); + } + + auto* exciter = new Esdc1a(getBus(bus_index), excitedata); + + if (excitedata.ports.contains(Esdc1aPorts::speed)) + { + IdxT speed = excitedata.ports.at(Esdc1aPorts::speed); + constexpr auto OMEGA = Esdc1aExternalVariables::OMEGA; + exciter->getSignals().template attachSignalNode(getSignal(speed)); + } + + if (excitedata.ports.contains(Esdc1aPorts::efd)) + { + IdxT efd = excitedata.ports.at(Esdc1aPorts::efd); + constexpr auto EFD = Esdc1aInternalVariables::EFD; + exciter->getSignals().template assignSignalNode(getSignal(efd)); + } + + if (excitedata.ports.contains(Esdc1aPorts::vs)) + { + IdxT vs = excitedata.ports.at(Esdc1aPorts::vs); + constexpr auto VS = Esdc1aExternalVariables::VS; + exciter->getSignals().template attachSignalNode(getSignal(vs)); + } + + if (excitedata.ports.contains(Esdc1aPorts::vuel)) + { + IdxT vuel = excitedata.ports.at(Esdc1aPorts::vuel); + constexpr auto VUEL = Esdc1aExternalVariables::VUEL; + exciter->getSignals().template attachSignalNode(getSignal(vuel)); + } + + addComponent(exciter); + } + for (const auto& excitedata : data.sexspti) { IdxT bus_index = 0; diff --git a/docs/Figures/PhasorDynamics/ESDC1A_diagram.png b/docs/Figures/PhasorDynamics/ESDC1A_diagram.png new file mode 100644 index 000000000..c77a100a6 Binary files /dev/null and b/docs/Figures/PhasorDynamics/ESDC1A_diagram.png differ diff --git a/docs/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md b/docs/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md new file mode 100644 index 000000000..aa9a4c776 --- /dev/null +++ b/docs/GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md @@ -0,0 +1,6 @@ +# ESDC1A + +```{include} ../../../../../../GridKit/Model/PhasorDynamics/Exciter/ESDC1A/README.md +:start-line: 1 +:relative-images: +``` diff --git a/docs/GridKit/Model/PhasorDynamics/Exciter/README.md b/docs/GridKit/Model/PhasorDynamics/Exciter/README.md index 2980fd53c..8e11df505 100644 --- a/docs/GridKit/Model/PhasorDynamics/Exciter/README.md +++ b/docs/GridKit/Model/PhasorDynamics/Exciter/README.md @@ -6,6 +6,7 @@ :hidden: ESAC6A +ESDC1A IEEET1 EXDC1 ESDC2A diff --git a/tests/UnitTests/PhasorDynamics/CMakeLists.txt b/tests/UnitTests/PhasorDynamics/CMakeLists.txt index 3d2973b62..a3633dbfa 100644 --- a/tests/UnitTests/PhasorDynamics/CMakeLists.txt +++ b/tests/UnitTests/PhasorDynamics/CMakeLists.txt @@ -96,6 +96,14 @@ target_link_libraries( GridKit::phasor_dynamics_components_dependency_tracking GridKit::testing) +add_executable(test_phasor_exciter_esdc1a runExciterEsdc1aTests.cpp) +target_link_libraries( + test_phasor_exciter_esdc1a + GridKit::definitions + GridKit::phasor_dynamics_components + GridKit::phasor_dynamics_components_dependency_tracking + GridKit::testing) + add_executable(test_phasor_exciter_sexspti runExciterSexsPtiTests.cpp) target_link_libraries( test_phasor_exciter_sexspti @@ -147,6 +155,7 @@ 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 PhasorDynamicsExciterIeeet1Test COMMAND test_phasor_exciter_ieeet1) +add_test(NAME PhasorDynamicsExciterEsdc1aTest COMMAND test_phasor_exciter_esdc1a) add_test(NAME PhasorDynamicsGensalTest COMMAND test_phasor_gensal) add_test(NAME PhasorDynamicsExciterSexsPtiTest COMMAND test_phasor_exciter_sexspti) add_test(NAME PhasorDynamicsStabilizerIeeestTest COMMAND test_phasor_stabilizer_ieeest) @@ -169,6 +178,7 @@ install( test_phasor_genrou test_phasor_governor_tgov1 test_phasor_exciter_ieeet1 + test_phasor_exciter_esdc1a test_phasor_gensal test_phasor_exciter_sexspti test_phasor_stabilizer_ieeest diff --git a/tests/UnitTests/PhasorDynamics/ExciterEsdc1aTests.hpp b/tests/UnitTests/PhasorDynamics/ExciterEsdc1aTests.hpp new file mode 100644 index 000000000..29d2d8505 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/ExciterEsdc1aTests.hpp @@ -0,0 +1,685 @@ +#pragma once + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace GridKit +{ + namespace Testing + { + template + class ExciterEsdc1aTests + { + public: + using RealT = typename PhasorDynamics::Component::RealT; + + ExciterEsdc1aTests() = default; + ~ExciterEsdc1aTests() = default; + + static constexpr ScalarT kTol = static_cast(1.0e-12); + + TestOutcome constructor() + { + using namespace PhasorDynamics::Exciter; + + TestStatus success = true; + + PhasorDynamics::Bus bus(3.0, 4.0); + + Esdc1a default_exciter(&bus); + success *= (default_exciter.size() == static_cast(Esdc1aInternalVariables::MAXIMUM)); + success *= (default_exciter.getMonitor() == nullptr); + + auto data = makeDefaultData(); + data.monitored_variables.insert(Esdc1aMonitorableVariables::efd); + Esdc1a data_exciter(&bus, data); + success *= (data_exciter.size() == static_cast(Esdc1aInternalVariables::MAXIMUM)); + success *= (data_exciter.getMonitor() != nullptr); + + PhasorDynamics::SignalNode efd_node; + ScalarT efd_value{0.0}; + IdxT efd_index = INVALID_INDEX; + efd_node.set(&efd_value, &efd_index); + + data_exciter.getSignals().template assignSignalNode(&efd_node); + data_exciter.allocate(); + data_exciter.tagDifferentiable(); + + success *= (data_exciter.verify() == 0); + success *= (data_exciter.tag()[idx(Esdc1aInternalVariables::EFDP)] == true); + success *= (data_exciter.tag()[idx(Esdc1aInternalVariables::VC)] == true); + success *= (data_exciter.tag()[idx(Esdc1aInternalVariables::VR)] == true); + success *= (data_exciter.tag()[idx(Esdc1aInternalVariables::VF)] == true); + success *= (data_exciter.tag()[idx(Esdc1aInternalVariables::XLL)] == true); + + return success.report(__func__); + } + + TestOutcome zeroInitialResidual() + { + using namespace PhasorDynamics::Exciter; + + TestStatus success = true; + + Fixture fixture(makeDefaultData()); + success *= fixture.allocateAndInitialize(1.2); + + const auto& f = fixture.exciter.getResidual(); + for (size_t i = 0; i < f.size(); ++i) + { + if (!isEqual(f[i], static_cast(0.0), kTol)) + { + std::cout << "Non-zero ESDC1A residual at index " << i << ": " << f[i] << "\n"; + success = false; + } + } + + success *= fixture.efd_node.linked(); + success *= (fixture.efd_node.getVariableIndex() + == static_cast(idx(Esdc1aInternalVariables::EFD))); + success *= isEqual(fixture.efd_node.read(), static_cast(1.2), kTol); + success *= smoothHighValueGateInitialResidual(); + + return success.report(__func__); + } + + TestOutcome blockDiagramSemantics() + { + TestStatus success = true; + + success *= voltageErrorSummingJunction(); + success *= speedMultiplierSelector(); + success *= leadLagBlockSemantics(); + success *= timeConstantClampSemantics(); + success *= uelRoutingSelector(); + success *= exciterFeedbackLimiter(); + + return success.report(__func__); + } + + TestOutcome parameterValidation() + { + using namespace PhasorDynamics::Exciter; + + TestStatus success = true; + + PhasorDynamics::Bus bus(1.0, 0.0); + PhasorDynamics::SignalNode efd_node; + ScalarT efd_value{0.0}; + IdxT efd_index = INVALID_INDEX; + efd_node.set(&efd_value, &efd_index); + + auto valid = makeDefaultData(); + Esdc1a valid_model(&bus, valid); + valid_model.getSignals().template assignSignalNode(&efd_node); + valid_model.allocate(); + success *= (valid_model.verify() == 0); + + auto invalid_ta = makeDefaultData(); + invalid_ta.parameters[Esdc1aParameters::Ta] = 0.0; + Esdc1a invalid_ta_model(&bus, invalid_ta); + invalid_ta_model.getSignals().template assignSignalNode(&efd_node); + invalid_ta_model.allocate(); + success *= (invalid_ta_model.verify() > 0); + + auto zero_tb_nonzero_tc = makeDefaultData(); + zero_tb_nonzero_tc.parameters[Esdc1aParameters::Tc] = 0.1; + Esdc1a zero_tb_nonzero_tc_model(&bus, zero_tb_nonzero_tc); + zero_tb_nonzero_tc_model.getSignals().template assignSignalNode(&efd_node); + zero_tb_nonzero_tc_model.allocate(); + success *= (zero_tb_nonzero_tc_model.verify() == 0); + + auto invalid_tc = makeDefaultData(); + invalid_tc.parameters[Esdc1aParameters::Tc] = -0.1; + Esdc1a invalid_tc_model(&bus, invalid_tc); + invalid_tc_model.getSignals().template assignSignalNode(&efd_node); + invalid_tc_model.allocate(); + success *= (invalid_tc_model.verify() > 0); + + auto invalid_saturation = makeDefaultData(); + invalid_saturation.parameters[Esdc1aParameters::Se1] = 0.0; + invalid_saturation.parameters[Esdc1aParameters::Se2] = 0.33; + Esdc1a invalid_saturation_model(&bus, invalid_saturation); + invalid_saturation_model.getSignals().template assignSignalNode(&efd_node); + invalid_saturation_model.allocate(); + success *= (invalid_saturation_model.verify() > 0); + + auto missing_efd = makeDefaultData(); + Esdc1a missing_efd_model(&bus, missing_efd); + missing_efd_model.allocate(); + success *= (missing_efd_model.verify() > 0); + + auto missing_speed = makeDefaultData(); + missing_speed.parameters[Esdc1aParameters::Spdmlt] = 1.0; + Esdc1a missing_speed_model(&bus, missing_speed); + missing_speed_model.getSignals().template assignSignalNode(&efd_node); + missing_speed_model.allocate(); + success *= (missing_speed_model.verify() > 0); + + return success.report(__func__); + } + +#ifdef GRIDKIT_ENABLE_ENZYME + TestOutcome jacobianStructureAndValues() + { + TestStatus success = true; + + const auto tol = static_cast(1.0e-9); + + auto data = makeDefaultData(); + data.parameters[Params::Spdmlt] = 1.0; + data.parameters[Params::UEL] = static_cast(2); + + auto dependency_tracking_jacobian = dependencyTrackingJacobian(data); + auto enzyme_jacobian = enzymeJacobian(data); + + 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], tol); + } + + return success.report(__func__); + } +#endif + + private: + using Internal = PhasorDynamics::Exciter::Esdc1aInternalVariables; + using External = PhasorDynamics::Exciter::Esdc1aExternalVariables; + using Params = PhasorDynamics::Exciter::Esdc1aParameters; + using DataT = PhasorDynamics::Exciter::Esdc1aData; + + static size_t idx(Internal variable) + { + return static_cast(variable); + } + + auto makeDefaultData() -> DataT + { + DataT data; + data.device_class = "exciter"; + data.disambiguation_string = "esdc1a_test"; + + data.parameters[Params::Tr] = 0.0; + data.parameters[Params::Ka] = 40.0; + data.parameters[Params::Ta] = 0.1; + data.parameters[Params::Tb] = 0.0; + data.parameters[Params::Tc] = 0.0; + data.parameters[Params::Vrmax] = 1.0; + data.parameters[Params::Vrmin] = -1.0; + data.parameters[Params::Ke] = 0.1; + data.parameters[Params::Te] = 0.5; + data.parameters[Params::Kf] = 0.05; + data.parameters[Params::Tf1] = 0.7; + data.parameters[Params::Spdmlt] = 0.0; + data.parameters[Params::E1] = 2.8; + data.parameters[Params::Se1] = 0.08; + data.parameters[Params::E2] = 3.7; + data.parameters[Params::Se2] = 0.33; + data.parameters[Params::UEL] = static_cast(0); + data.parameters[Params::exclim] = 1.0; + + return data; + } + + struct Fixture + { + using BusT = PhasorDynamics::Bus; + using SignalT = PhasorDynamics::SignalNode; + using ExciterT = PhasorDynamics::Exciter::Esdc1a; + + DataT data; + BusT bus; + SignalT efd_node; + SignalT omega_node; + SignalT vs_node; + SignalT vuel_node; + + ScalarT efd_value{0.0}; + ScalarT omega_value{0.0}; + ScalarT vs_value{0.0}; + ScalarT vuel_value{-2.0}; + + IdxT efd_index{INVALID_INDEX}; + IdxT omega_index{20}; + IdxT vs_index{21}; + IdxT vuel_index{22}; + + ExciterT exciter; + + explicit Fixture(const DataT& data_in) + : data(data_in), + bus(3.0, 4.0), + exciter(&bus, data) + { + efd_node.set(&efd_value, &efd_index); + omega_node.set(&omega_value, &omega_index); + vs_node.set(&vs_value, &vs_index); + vuel_node.set(&vuel_value, &vuel_index); + + exciter.getSignals().template assignSignalNode(&efd_node); + exciter.getSignals().template attachSignalNode(&omega_node); + exciter.getSignals().template attachSignalNode(&vs_node); + exciter.getSignals().template attachSignalNode(&vuel_node); + } + + bool allocateAndInitialize(ScalarT efd0) + { + bus.allocate(); + bus.initialize(); + exciter.allocate(); + efd_node.init(efd0); + return exciter.verify() == 0 + && exciter.initialize() == 0 + && exciter.evaluateResidual() == 0; + } + }; + + bool voltageErrorSummingJunction() + { + Fixture fixture(makeDefaultData()); + if (!fixture.allocateAndInitialize(1.2)) + { + return false; + } + + fixture.vs_value += 0.1; + fixture.exciter.evaluateResidual(); + bool success = fixture.exciter.getResidual()[idx(Internal::EV)] > static_cast(0.0); + + fixture.vs_value -= 0.1; + fixture.exciter.y()[idx(Internal::VC)] += 0.1; + fixture.exciter.evaluateResidual(); + success = success && fixture.exciter.getResidual()[idx(Internal::EV)] < static_cast(0.0); + + fixture.exciter.y()[idx(Internal::VC)] -= 0.1; + fixture.exciter.y()[idx(Internal::VF)] += 0.1; + fixture.exciter.evaluateResidual(); + success = success && fixture.exciter.getResidual()[idx(Internal::EV)] < static_cast(0.0); + + return success; + } + + bool speedMultiplierSelector() + { + auto disabled_data = makeDefaultData(); + Fixture disabled(disabled_data); + if (!disabled.allocateAndInitialize(1.2)) + { + return false; + } + + disabled.omega_value = 0.05; + disabled.exciter.evaluateResidual(); + bool success = isEqual(disabled.exciter.getResidual()[idx(Internal::EFD)], + static_cast(0.0), + kTol); + + auto enabled_data = makeDefaultData(); + enabled_data.parameters[Params::Spdmlt] = 1.0; + Fixture enabled(enabled_data); + if (!enabled.allocateAndInitialize(1.2)) + { + return false; + } + + enabled.omega_value = 0.05; + enabled.exciter.evaluateResidual(); + success = success && enabled.exciter.getResidual()[idx(Internal::EFD)] > static_cast(0.0); + + return success; + } + + bool leadLagBlockSemantics() + { + Fixture clamped(makeDefaultData()); + if (!clamped.allocateAndInitialize(1.2)) + { + return false; + } + + clamped.exciter.y()[idx(Internal::VLL)] += 0.1; + clamped.exciter.evaluateResidual(); + bool success = clamped.exciter.getResidual()[idx(Internal::VLL)] < static_cast(0.0); + + auto active_data = makeDefaultData(); + active_data.parameters[Params::Tb] = 0.5; + active_data.parameters[Params::Tc] = 0.2; + Fixture active(active_data); + if (!active.allocateAndInitialize(1.2)) + { + return false; + } + + active.exciter.y()[idx(Internal::EV)] += 0.1; + active.exciter.evaluateResidual(); + success = success && active.exciter.getResidual()[idx(Internal::VLL)] > static_cast(0.0); + + active.exciter.y()[idx(Internal::EV)] -= 0.1; + active.exciter.y()[idx(Internal::VLL)] += 0.1; + active.exciter.evaluateResidual(); + success = success && active.exciter.getResidual()[idx(Internal::VLL)] < static_cast(0.0); + + return success; + } + + bool smoothHighValueGateInitialResidual() + { + auto data = makeDefaultData(); + + PhasorDynamics::Bus bus(3.0, 4.0); + PhasorDynamics::SignalNode efd_node; + PhasorDynamics::SignalNode omega_node; + PhasorDynamics::SignalNode vs_node; + + ScalarT efd_value{0.0}; + ScalarT omega_value{0.0}; + ScalarT vs_value{0.0}; + + IdxT efd_index{INVALID_INDEX}; + IdxT omega_index{20}; + IdxT vs_index{21}; + + efd_node.set(&efd_value, &efd_index); + omega_node.set(&omega_value, &omega_index); + vs_node.set(&vs_value, &vs_index); + + PhasorDynamics::Exciter::Esdc1a exciter(&bus, data); + exciter.getSignals().template assignSignalNode(&efd_node); + exciter.getSignals().template attachSignalNode(&omega_node); + exciter.getSignals().template attachSignalNode(&vs_node); + + bus.allocate(); + bus.initialize(); + exciter.allocate(); + efd_node.init(1.2); + + TestStatus success = true; + success *= (exciter.verify() == 0); + success *= (exciter.initialize() == 0); + success *= (exciter.evaluateResidual() == 0); + success *= isEqual(exciter.getResidual()[idx(Internal::VHV)], static_cast(0.0), kTol); + + return success; + } + + bool timeConstantClampSemantics() + { + auto data = makeDefaultData(); + data.parameters[Params::Tr] = 0.0; + data.parameters[Params::Tb] = 0.0; + data.parameters[Params::Tc] = 0.1; + data.parameters[Params::Tf1] = 0.0; + + Fixture fixture(data); + if (!fixture.allocateAndInitialize(1.2)) + { + return false; + } + + fixture.exciter.tagDifferentiable(); + + bool success = fixture.exciter.tag()[idx(Internal::VC)]; + success = success && fixture.exciter.tag()[idx(Internal::VF)]; + success = success && fixture.exciter.tag()[idx(Internal::XLL)]; + + fixture.exciter.y()[idx(Internal::VC)] += 0.1; + fixture.exciter.evaluateResidual(); + success = success + && fixture.exciter.getResidual()[idx(Internal::VC)] + < static_cast(0.0); + + fixture.exciter.y()[idx(Internal::VC)] -= 0.1; + fixture.exciter.y()[idx(Internal::VF)] += 0.1; + fixture.exciter.evaluateResidual(); + success = success + && fixture.exciter.getResidual()[idx(Internal::VF)] + < static_cast(0.0); + + fixture.exciter.y()[idx(Internal::VF)] -= 0.1; + fixture.exciter.y()[idx(Internal::XLL)] += 0.1; + fixture.exciter.evaluateResidual(); + success = success + && fixture.exciter.getResidual()[idx(Internal::XLL)] + < static_cast(0.0); + + return success; + } + + bool uelRoutingSelector() + { + Fixture hv_gate(makeDefaultData()); + if (!hv_gate.allocateAndInitialize(1.2)) + { + return false; + } + + hv_gate.vuel_value = hv_gate.exciter.y()[idx(Internal::VLL)] + 0.1; + hv_gate.exciter.evaluateResidual(); + bool success = hv_gate.exciter.getResidual()[idx(Internal::VHV)] > static_cast(0.0); + success = success && isEqual(hv_gate.exciter.getResidual()[idx(Internal::EV)], static_cast(0.0), kTol); + + auto sum_data = makeDefaultData(); + sum_data.parameters[Params::UEL] = static_cast(2); + Fixture sum_junction(sum_data); + sum_junction.vuel_value = 0.0; + if (!sum_junction.allocateAndInitialize(1.2)) + { + return false; + } + + sum_junction.vuel_value = 0.1; + sum_junction.exciter.evaluateResidual(); + success = success && sum_junction.exciter.getResidual()[idx(Internal::EV)] > static_cast(0.0); + success = success && isEqual(sum_junction.exciter.getResidual()[idx(Internal::VHV)], static_cast(0.0), kTol); + + return success; + } + + bool exciterFeedbackLimiter() + { + auto limited_data = makeDefaultData(); + limited_data.parameters[Params::Ke] = -0.2; + limited_data.parameters[Params::Se1] = 0.0; + limited_data.parameters[Params::Se2] = 0.0; + limited_data.parameters[Params::exclim] = 1.0; + Fixture limited(limited_data); + if (!limited.allocateAndInitialize(1.2)) + { + return false; + } + + limited.exciter.y()[idx(Internal::EFDP)] = 1.0; + limited.exciter.y()[idx(Internal::SE)] = 0.0; + limited.exciter.y()[idx(Internal::VFE)] = 0.0; + limited.exciter.evaluateResidual(); + bool success = std::abs(limited.exciter.getResidual()[idx(Internal::VFE)]) < kTol; + + auto unlimited_data = limited_data; + unlimited_data.parameters[Params::exclim] = 0.0; + Fixture unlimited(unlimited_data); + if (!unlimited.allocateAndInitialize(1.2)) + { + return false; + } + + unlimited.exciter.y()[idx(Internal::EFDP)] = 1.0; + unlimited.exciter.y()[idx(Internal::SE)] = 0.0; + unlimited.exciter.y()[idx(Internal::VFE)] = 0.0; + unlimited.exciter.evaluateResidual(); + success = success && unlimited.exciter.getResidual()[idx(Internal::VFE)] < static_cast(0.0); + + return success; + } + +#ifdef GRIDKIT_ENABLE_ENZYME + std::vector + dependencyTrackingJacobian(const DataT& data) + { + using Variable = DependencyTracking::Variable; + + PhasorDynamics::Bus bus(Variable{3.0}, Variable{4.0}); + PhasorDynamics::SignalNode efd_node; + PhasorDynamics::SignalNode omega_node; + PhasorDynamics::SignalNode vs_node; + PhasorDynamics::SignalNode vuel_node; + + Variable efd_value{0.0}; + Variable omega_value{0.0}; + Variable vs_value{0.0}; + Variable vuel_value{0.0}; + + IdxT efd_index = INVALID_INDEX; + IdxT omega_index = 13; + IdxT vs_index = 14; + IdxT vuel_index = 15; + + efd_node.set(&efd_value, &efd_index); + omega_node.set(&omega_value, &omega_index); + vs_node.set(&vs_value, &vs_index); + vuel_node.set(&vuel_value, &vuel_index); + + PhasorDynamics::Exciter::Esdc1a exciter(&bus, data); + exciter.getSignals().template assignSignalNode(&efd_node); + exciter.getSignals().template attachSignalNode(&omega_node); + exciter.getSignals().template attachSignalNode(&vs_node); + exciter.getSignals().template attachSignalNode(&vuel_node); + + bus.allocate(); + exciter.allocate(); + bus.initialize(); + efd_node.init(Variable{1.2}); + exciter.initialize(); + + for (size_t i = 0; i < exciter.size(); ++i) + { + exciter.y()[i].setVariableNumber(i); + } + for (size_t i = 0; i < bus.size(); ++i) + { + bus.y()[i].setVariableNumber(i + exciter.size()); + } + omega_value.setVariableNumber(13); + vs_value.setVariableNumber(14); + vuel_value.setVariableNumber(15); + + bus.evaluateResidual(); + exciter.evaluateResidual(); + auto residual_y = exciter.getResidual(); + + omega_value = 0.0; + vs_value = 0.0; + vuel_value = 0.0; + bus.initialize(); + efd_node.init(Variable{1.2}); + exciter.initialize(); + + for (size_t i = 0; i < exciter.size(); ++i) + { + exciter.yp()[i].setVariableNumber(i); + } + + bus.evaluateResidual(); + exciter.evaluateResidual(); + auto residual_yp = exciter.getResidual(); + + std::vector dependencies(residual_y.size()); + for (IdxT i = 0; i < residual_y.size(); ++i) + { + auto dependency_y = residual_y[static_cast(i)].getDependencies(); + auto dependency_yp = residual_yp[static_cast(i)].getDependencies(); + + for (const auto& pair_y : dependency_y) + { + auto index_y = pair_y.first; + auto value_y = pair_y.second; + auto it_yp = dependency_yp.find(index_y); + if (it_yp != dependency_yp.end()) + { + dependencies[static_cast(i)].insert(std::make_pair(index_y, value_y + it_yp->second)); + } + else + { + dependencies[static_cast(i)].insert(pair_y); + } + } + + for (const auto& pair_yp : dependency_yp) + { + if (dependency_y.find(pair_yp.first) == dependency_y.end()) + { + dependencies[static_cast(i)].insert(pair_yp); + } + } + } + + return dependencies; + } + + std::vector + enzymeJacobian(const DataT& data) + { + PhasorDynamics::Bus bus(3.0, 4.0); + PhasorDynamics::SignalNode efd_node; + PhasorDynamics::SignalNode omega_node; + PhasorDynamics::SignalNode vs_node; + PhasorDynamics::SignalNode vuel_node; + + ScalarT efd_value{0.0}; + ScalarT omega_value{0.0}; + ScalarT vs_value{0.0}; + ScalarT vuel_value{0.0}; + + IdxT efd_index = INVALID_INDEX; + IdxT omega_index = 13; + IdxT vs_index = 14; + IdxT vuel_index = 15; + + efd_node.set(&efd_value, &efd_index); + omega_node.set(&omega_value, &omega_index); + vs_node.set(&vs_value, &vs_index); + vuel_node.set(&vuel_value, &vuel_index); + + PhasorDynamics::Exciter::Esdc1a exciter(&bus, data); + exciter.getSignals().template assignSignalNode(&efd_node); + exciter.getSignals().template attachSignalNode(&omega_node); + exciter.getSignals().template attachSignalNode(&vs_node); + exciter.getSignals().template attachSignalNode(&vuel_node); + + bus.allocate(); + exciter.allocate(); + bus.initialize(); + efd_node.init(1.2); + exciter.initialize(); + exciter.updateTime(0.0, 1.0); + + for (size_t i = 0; i < bus.size(); ++i) + { + bus.setVariableIndex(i, static_cast(i + exciter.size())); + bus.setResidualIndex(i, static_cast(i + exciter.size())); + } + + bus.evaluateResidual(); + exciter.evaluateResidual(); + exciter.evaluateJacobian(); + exciter.constructCsr(); + + auto* model_jacobian = exciter.getCsrJacobian(); + + return MapFromCsr(model_jacobian); + } +#endif + }; + } // namespace Testing +} // namespace GridKit diff --git a/tests/UnitTests/PhasorDynamics/runExciterEsdc1aTests.cpp b/tests/UnitTests/PhasorDynamics/runExciterEsdc1aTests.cpp new file mode 100644 index 000000000..04c40b955 --- /dev/null +++ b/tests/UnitTests/PhasorDynamics/runExciterEsdc1aTests.cpp @@ -0,0 +1,18 @@ +#include "ExciterEsdc1aTests.hpp" + +int main() +{ + GridKit::Testing::TestingResults result; + + GridKit::Testing::ExciterEsdc1aTests test; + + result += test.constructor(); + result += test.zeroInitialResidual(); + result += test.blockDiagramSemantics(); + result += test.parameterValidation(); +#ifdef GRIDKIT_ENABLE_ENZYME + result += test.jacobianStructureAndValues(); +#endif + + return result.summary(); +}